MOM_wave_speed.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Routines for calculating baroclinic wave speeds
7
9use mom_error_handler, only : mom_error, fatal, warning
10use mom_file_parser, only : log_version
11use mom_grid, only : ocean_grid_type
12use mom_interface_heights, only : thickness_to_dz
17use mom_eos, only : calculate_density_derivs, calculate_specific_vol_derivs
18
19implicit none ; private
20
21#include <MOM_memory.h>
22
24
25! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
26! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
27! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
28! vary with the Boussinesq approximation, the Boussinesq variant is given first.
29
30!> Control structure for MOM_wave_speed
31type, public :: wave_speed_cs ; private
32 logical :: initialized = .false. !< True if this control structure has been initialized.
33 logical :: use_ebt_mode = .false. !< If true, calculate the equivalent barotropic wave speed instead
34 !! of the first baroclinic wave speed.
35 !! This parameter controls the default behavior of wave_speed() which
36 !! can be overridden by optional arguments.
37 logical :: better_cg1_est = .false. !< If true, use an improved estimate of the first mode
38 !! internal wave speed.
39 real :: mono_n2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as
40 !! monotonic for the purposes of calculating the equivalent barotropic
41 !! wave speed [nondim]. This parameter controls the default behavior of
42 !! wave_speed() which can be overridden by optional arguments.
43 real :: mono_n2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of
44 !! calculating the equivalent barotropic wave speed [H ~> m or kg m-2].
45 !! If this parameter is negative, this limiting does not occur.
46 !! This parameter controls the default behavior of wave_speed() which
47 !! can be overridden by optional arguments.
48 real :: min_speed2 = 0. !< The minimum mode 1 internal wave speed squared [L2 T-2 ~> m2 s-2]
49 real :: wave_speed_tol = 0.001 !< The fractional tolerance with which to solve for the wave
50 !! speeds [nondim]
51 real :: c1_thresh = -1.0 !< A minimal value of the first mode internal wave speed
52 !! below which all higher mode speeds are not calculated but
53 !! are simply reported as 0 [L T-1 ~> m s-1]. A non-negative
54 !! value must be specified via a call to wave_speed_init for
55 !! the subroutine wave_speeds to be used (but not wave_speed).
56 type(remapping_cs) :: remap_2018_cs !< Used for vertical remapping when calculating equivalent barotropic
57 !! mode structure for answer dates below 20190101.
58 type(remapping_cs) :: remap_cs !< Used for vertical remapping when calculating equivalent barotropic
59 !! mode structure.
60 integer :: remap_answer_date = 99991231 !< The vintage of the order of arithmetic and expressions to use
61 !! for remapping. Values below 20190101 recover the remapping
62 !! answers from 2018, while higher values use more robust
63 !! forms of the same remapping expressions.
64 type(diag_ctrl), pointer :: diag !< Diagnostics control structure
65end type wave_speed_cs
66
67contains
68
69!> Calculates the wave speed of the first baroclinic mode.
70subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N2_column_fraction, &
71 mono_N2_depth, modal_structure)
72 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
73 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
74 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
75 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
76 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
77 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
78 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cg1 !< First mode internal wave speed [L T-1 ~> m s-1]
79 type(wave_speed_cs), intent(in) :: cs !< Wave speed control struct
80 integer, optional, intent(in) :: halo_size !< Width of halo within which to
81 !! calculate wave speeds
82 logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
83 !! barotropic mode instead of the first baroclinic mode.
84 real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction
85 !! of water column over which N2 is limited as monotonic
86 !! for the purposes of calculating vertical modal structure [nondim].
87 real, optional, intent(in) :: mono_n2_depth !< A depth below which N2 is limited as
88 !! monotonic for the purposes of calculating vertical
89 !! modal structure [H ~> m or kg m-2].
90 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
91 optional, intent(out) :: modal_structure !< Normalized model structure [nondim]
92
93 ! Local variables
94 real, dimension(SZK_(GV)+1) :: &
95 drho_dt, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
96 drho_ds, & ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
97 dspv_dt, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
98 dspv_ds, & ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
99 pres, & ! Interface pressure [R L2 T-2 ~> Pa]
100 t_int, & ! Temperature interpolated to interfaces [C ~> degC]
101 s_int, & ! Salinity interpolated to interfaces [S ~> ppt]
102 h_top, & ! The distance of each filtered interface from the ocean surface [H ~> m or kg m-2]
103 h_bot, & ! The distance of each filtered interface from the bottom [H ~> m or kg m-2]
104 gprime ! The reduced gravity across each interface [L2 H-1 T-2 ~> m s-2 or m4 s-2 kg-1].
105 real, dimension(SZK_(GV)) :: &
106 igl, igu ! The inverse of the reduced gravity across an interface times
107 ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].
108 real, dimension(SZK_(GV),SZI_(G)) :: &
109 hf, & ! Layer thicknesses after very thin layers are combined [H ~> m or kg m-2]
110 tf, & ! Layer temperatures after very thin layers are combined [C ~> degC]
111 sf, & ! Layer salinities after very thin layers are combined [S ~> ppt]
112 rf ! Layer densities after very thin layers are combined [R ~> kg m-3]
113 real, dimension(SZK_(GV)) :: &
114 hc, & ! A column of layer thicknesses after convective instabilities are removed [H ~> m or kg m-2]
115 tc, & ! A column of layer temperatures after convective instabilities are removed [C ~> degC]
116 sc, & ! A column of layer salinities after convective instabilities are removed [S ~> ppt]
117 rc ! A column of layer densities after convective instabilities are removed [R ~> kg m-3]
118 real :: i_htot ! The inverse of the total filtered thicknesses [H-1 ~> m-1 or m2 kg-1]
119 real :: det, ddet ! Determinant of the eigen system and its derivative with lam. Because the
120 ! units of the eigenvalue change with the number of layers and because of the
121 ! dynamic rescaling that is used to keep det in a numerically representable range,
122 ! the units of of det are hard to interpret, but det/ddet is always in units
123 ! of [T2 L-2 ~> s2 m-2]
124 real :: lam ! The eigenvalue [T2 L-2 ~> s2 m-2]
125 real :: dlam ! The change in estimates of the eigenvalue [T2 L-2 ~> s2 m-2]
126 real :: lam0 ! The first guess of the eigenvalue [T2 L-2 ~> s2 m-2]
127 real :: h_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
128 real, dimension(SZI_(G)) :: &
129 htot, hmin, & ! Thicknesses [H ~> m or kg m-2]
130 h_here, & ! A thickness [H ~> m or kg m-2]
131 hxt_here, & ! A layer integrated temperature [C H ~> degC m or degC kg m-2]
132 hxs_here, & ! A layer integrated salinity [S H ~> ppt m or ppt kg m-2]
133 hxr_here ! A layer integrated density [R H ~> kg m-2 or kg2 m-5]
134 real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
135 real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2]
136 real :: cg1_est ! An initial estimate of the squared first mode speed [L2 T-2 ~> m2 s-2]
137 real :: i_hnew ! The inverse of a new layer thickness [H-1 ~> m-1 or m2 kg-1]
138 real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R H ~> kg m-2 or kg2 m-5]
139 real :: dspvxh_sum ! The sum of specific volume differences across interfaces times
140 ! thicknesses [H R-1 ~> m4 kg-1 or m], negative for stable stratification.
141 real :: g_rho0 ! G_Earth/Rho0 [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
142 real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant and
143 ! its derivative with lam between rows of the Thomas algorithm solver [L2 s2 T-2 m-2 ~> nondim].
144 ! The exact value should not matter for the final result if it is an even power of 2.
145 real :: tol_hfrac ! Layers that together are smaller than this fraction of
146 ! the total water column can be merged for efficiency [nondim].
147 real :: min_h_frac ! tol_Hfrac divided by the total number of layers [nondim].
148 real :: tol_solve ! The fractional tolerance with which to solve for the wave speeds [nondim]
149 real :: tol_merge ! The fractional change in estimated wave speed that is allowed
150 ! when deciding to merge layers in the calculation [nondim]
151 real :: rescale ! A rescaling factor to control the magnitude of the determinant [nondim]
152 real :: i_rescale ! The reciprocal of the rescaling factor to control the magnitude of the determinant [nondim]
153 integer :: kf(szi_(g)) ! The number of active layers after filtering.
154 integer, parameter :: max_itt = 10
155 real :: lam_it(max_itt) ! The guess at the eignevalue with each iteration [T2 L-2 ~> s2 m-2]
156 real :: det_it(max_itt), ddet_it(max_itt) ! The determinant of the matrix and its derivative with lam
157 ! with each iteration. Because of all of the dynamic rescaling of the determinant
158 ! between rows, its units are not easily interpretable, but the ratio of det/ddet
159 ! always has units of [T2 L-2 ~> s2 m-2]
160 logical :: use_eos ! If true, density or specific volume is calculated from T & S using an equation of state.
161 logical :: nonbous ! If true, do not make the Boussinesq approximation.
162 logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed.
163 logical :: merge ! If true, merge the current layer with the one above.
164 integer :: kc ! The number of layers in the column after merging
165 integer :: i, j, k, k2, itt, is, ie, js, je, nz, halo
166 real :: hw ! The mean of the adjacent layer thicknesses [H ~> m or kg m-2]
167 real :: sum_hc ! The sum of the layer thicknesses [H ~> m or kg m-2]
168 real :: gp ! A limited local copy of gprime [L2 H-1 T-2 ~> m s-2 or m4 s-2 kg-1]
169 real :: n2min ! A minimum buoyancy frequency, including a slope rescaling factor [L2 H-2 T-2 ~> s-2 or m6 kg-2 s-2]
170 logical :: below_mono_n2_frac ! True if an interface is below the fractional depth where N2 should not increase.
171 logical :: below_mono_n2_depth ! True if an interface is below the absolute depth where N2 should not increase.
172 logical :: l_use_ebt_mode, calc_modal_structure
173 real :: l_mono_n2_column_fraction ! A local value of mono_N2_column_fraction [nondim]
174 real :: l_mono_n2_depth ! A local value of mono_N2_column_depth [H ~> m or kg m-2]
175 real :: mode_struct(szk_(gv)) ! The mode structure [nondim], but it is also temporarily
176 ! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6.
177 real :: ms_min, ms_max ! The minimum and maximum mode structure values returned from tdma6 [L2 T-2 ~> m2 s-2]
178 real :: ms_sq ! The sum of the square of the values returned from tdma6 [L4 T-4 ~> m4 s-4]
179
180 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke ; halo = 0
181
182 if (.not. cs%initialized) call mom_error(fatal, "MOM_wave_speed / wave_speed: "// &
183 "Module must be initialized before it is used.")
184
185 if (present(halo_size)) then
186 halo = halo_size
187 is = g%isc - halo ; ie = g%iec + halo ; js = g%jsc - halo ; je = g%jec + halo
188 endif
189
190 l_use_ebt_mode = cs%use_ebt_mode
191 if (present(use_ebt_mode)) l_use_ebt_mode = use_ebt_mode
192 l_mono_n2_column_fraction = cs%mono_N2_column_fraction
193 if (present(mono_n2_column_fraction)) l_mono_n2_column_fraction = mono_n2_column_fraction
194 l_mono_n2_depth = cs%mono_N2_depth
195 if (present(mono_n2_depth)) l_mono_n2_depth = mono_n2_depth
196 calc_modal_structure = l_use_ebt_mode
197 if (present(modal_structure)) calc_modal_structure = .true.
198 if (calc_modal_structure) then
199 do k=1,nz ; do j=js,je ; do i=is,ie
200 modal_structure(i,j,k) = 0.0
201 enddo ; enddo ; enddo
202 endif
203
204 nonbous = .not.(gv%Boussinesq .or. gv%semi_Boussinesq)
205 h_to_pres = gv%H_to_RZ * gv%g_Earth
206 ! Note that g_Rho0 = H_to_pres / GV%Rho0**2
207 if (.not.nonbous) g_rho0 = gv%g_Earth*gv%H_to_Z / gv%Rho0
208 use_eos = associated(tv%eqn_of_state)
209
210 better_est = cs%better_cg1_est
211
212 if (better_est) then
213 tol_solve = cs%wave_speed_tol
214 tol_hfrac = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
215 else
216 tol_solve = 0.001 ; tol_hfrac = 0.0001 ; tol_merge = 0.001
217 endif
218
219 ! The rescaling below can control the growth of the determinant provided that
220 ! (tol_merge*cg1_min2/c2_scale > I_rescale). For default values, this suggests a stable lower
221 ! bound on min_speed of sqrt(nz/(tol_solve*rescale)) or 3e2/1024**2 = 2.9e-4 m/s for 90 layers.
222 ! The upper bound on the rate of increase in the determinant is g'H/c2_scale < rescale or in the
223 ! worst possible oceanic case of g'H < 0.5*10m/s2*1e4m = 5.e4 m2/s2 < 1024**2*c2_scale, suggesting
224 ! that c2_scale can safely be set to 1/(16*1024**2), which would decrease the stable floor on
225 ! min_speed to ~6.9e-8 m/s for 90 layers or 2.33e-7 m/s for 1000 layers.
226 cg1_min2 = cs%min_speed2
227 rescale = 1024.0**4 ; i_rescale = 1.0/rescale
228 c2_scale = us%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.
229
230 min_h_frac = tol_hfrac / real(nz)
231 !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,tv,use_EOS,nonBous, &
232 !$OMP CS,min_h_frac,calc_modal_structure,l_use_ebt_mode, &
233 !$OMP modal_structure,l_mono_N2_column_fraction,l_mono_N2_depth, &
234 !$OMP H_to_pres,cg1,g_Rho0,rescale,I_rescale,cg1_min2, &
235 !$OMP better_est,tol_solve,tol_merge,c2_scale)
236 do j=js,je
237 ! First merge very thin layers with the one above (or below if they are
238 ! at the top). This also transposes the row order so that columns can
239 ! be worked upon one at a time.
240 do i=is,ie ; htot(i) = 0.0 ; enddo
241 do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
242
243 do i=is,ie
244 hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
245 hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
246 enddo
247 if (use_eos) then
248 do k=1,nz ; do i=is,ie
249 if ((h_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then
250 hf(kf(i),i) = h_here(i)
251 tf(kf(i),i) = hxt_here(i) / h_here(i)
252 sf(kf(i),i) = hxs_here(i) / h_here(i)
253 kf(i) = kf(i) + 1
254
255 ! Start a new layer
256 h_here(i) = h(i,j,k)
257 hxt_here(i) = h(i,j,k) * tv%T(i,j,k)
258 hxs_here(i) = h(i,j,k) * tv%S(i,j,k)
259 else
260 h_here(i) = h_here(i) + h(i,j,k)
261 hxt_here(i) = hxt_here(i) + h(i,j,k) * tv%T(i,j,k)
262 hxs_here(i) = hxs_here(i) + h(i,j,k) * tv%S(i,j,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)
267 tf(kf(i),i) = hxt_here(i) / h_here(i)
268 sf(kf(i),i) = hxs_here(i) / h_here(i)
269 endif ; enddo
270 else ! .not. (use_EOS)
271 do k=1,nz ; do i=is,ie
272 if ((h_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then
273 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
274 kf(i) = kf(i) + 1
275
276 ! Start a new layer
277 h_here(i) = h(i,j,k)
278 hxr_here(i) = h(i,j,k)*gv%Rlay(k)
279 else
280 h_here(i) = h_here(i) + h(i,j,k)
281 hxr_here(i) = hxr_here(i) + h(i,j,k)*gv%Rlay(k)
282 endif
283 enddo ; enddo
284 do i=is,ie ; if (h_here(i) > 0.0) then
285 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
286 endif ; enddo
287 endif
288
289 ! From this point, we can work on individual columns without causing memory to have page faults.
290 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
291 if (use_eos) then
292 pres(1) = 0.0 ; h_top(1) = 0.0
293 do k=2,kf(i)
294 pres(k) = pres(k-1) + h_to_pres*hf(k-1,i)
295 t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
296 s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
297 h_top(k) = h_top(k-1) + hf(k-1,i)
298 enddo
299 if (nonbous) then
300 call calculate_specific_vol_derivs(t_int, s_int, pres, dspv_dt, dspv_ds, &
301 tv%eqn_of_state, (/2,kf(i)/) )
302 else
303 call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, &
304 tv%eqn_of_state, (/2,kf(i)/) )
305 endif
306
307 ! Sum the reduced gravities to find out how small a density difference is negligibly small.
308 drxh_sum = 0.0 ; dspvxh_sum = 0.0
309 if (better_est) then
310 ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for
311 ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.
312 ! For a uniform stratification and a huge number of layers uniformly distributed in
313 ! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.
314 if (h_top(kf(i)) > 0.0) then
315 i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
316 h_bot(kf(i)+1) = 0.0
317 if (nonbous) then
318 do k=kf(i),2,-1
319 h_bot(k) = h_bot(k+1) + hf(k,i)
320 dspvxh_sum = dspvxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
321 min(0.0, dspv_dt(k)*(tf(k,i)-tf(k-1,i)) + dspv_ds(k)*(sf(k,i)-sf(k-1,i)))
322 enddo
323 else
324 do k=kf(i),2,-1
325 h_bot(k) = h_bot(k+1) + hf(k,i)
326 drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
327 max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
328 enddo
329 endif
330 endif
331 else
332 ! This estimate is problematic in that it goes like 1/nz for a large number of layers,
333 ! but it is an overestimate (as desired) for a small number of layers, by at a factor
334 ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.
335 if (nonbous) then
336 do k=2,kf(i)
337 dspvxh_sum = dspvxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
338 min(0.0, dspv_dt(k)*(tf(k,i)-tf(k-1,i)) + dspv_ds(k)*(sf(k,i)-sf(k-1,i)))
339 enddo
340 else
341 do k=2,kf(i)
342 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
343 max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
344 enddo
345 endif
346 endif
347 else ! .not. (use_EOS)
348 drxh_sum = 0.0 ; dspvxh_sum = 0.0
349 if (better_est) then
350 h_top(1) = 0.0
351 do k=2,kf(i) ; h_top(k) = h_top(k-1) + hf(k-1,i) ; enddo
352 if (h_top(kf(i)) > 0.0) then
353 i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
354 h_bot(kf(i)+1) = 0.0
355 if (nonbous) then
356 do k=kf(i),2,-1
357 h_bot(k) = h_bot(k+1) + hf(k,i)
358 dspvxh_sum = dspvxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
359 min(0.0, (rf(k-1,i)-rf(k,i)) / (rf(k,i)*rf(k-1,i)))
360 enddo
361 else
362 do k=kf(i),2,-1
363 h_bot(k) = h_bot(k+1) + hf(k,i)
364 drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * max(0.0,rf(k,i)-rf(k-1,i))
365 enddo
366 endif
367 endif
368 else
369 if (nonbous) then
370 do k=2,kf(i)
371 dspvxh_sum = dspvxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
372 min(0.0, (rf(k-1,i)-rf(k,i)) / (rf(k,i)*rf(k-1,i)))
373 enddo
374 else
375 do k=2,kf(i)
376 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
377 enddo
378 endif
379 endif
380 endif ! use_EOS
381
382 if (nonbous) then
383 ! Note that dSpVxh_sum is negative for stable stratification.
384 cg1_est = h_to_pres * abs(dspvxh_sum)
385 else
386 cg1_est = g_rho0 * drxh_sum
387 endif
388
389 ! Find gprime across each internal interface, taking care of convective instabilities by
390 ! merging layers. If the estimated wave speed is too small, simply return zero.
391 if (cg1_est <= cg1_min2) then
392 cg1(i,j) = 0.0
393 if (present(modal_structure)) modal_structure(i,j,:) = 0.
394 else
395 ! Merge layers to eliminate convective instabilities or exceedingly
396 ! small reduced gravities. Merging layers reduces the estimated wave speed by
397 ! (rho(2)-rho(1))*h(1)*h(2) / H_tot.
398 if (use_eos) then
399 kc = 1
400 hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
401 do k=2,kf(i)
402 if (better_est .and. nonbous) then
403 merge = ((dspv_dt(k)*(tc(kc)-tf(k,i)) + dspv_ds(k)*(sc(kc)-sf(k,i))) * &
404 ((hc(kc) * hf(k,i))*i_htot) < abs(2.0 * tol_merge * dspvxh_sum))
405 elseif (better_est) then
406 merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
407 ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
408 elseif (nonbous) then
409 merge = ((dspv_dt(k)*(tc(kc)-tf(k,i)) + dspv_ds(k)*(sc(kc)-sf(k,i))) * &
410 (hc(kc) + hf(k,i)) < abs(2.0 * tol_merge * dspvxh_sum))
411 else
412 merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
413 (hc(kc) + hf(k,i)) < 2.0 * tol_merge*drxh_sum)
414 endif
415 if (merge) then
416 ! Merge this layer with the one above and backtrack.
417 i_hnew = 1.0 / (hc(kc) + hf(k,i))
418 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
419 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
420 hc(kc) = (hc(kc) + hf(k,i))
421 ! Backtrack to remove any convective instabilities above... Note
422 ! that the tolerance is a factor of two larger, to avoid limit how
423 ! far back we go.
424 do k2=kc,2,-1
425 if (better_est .and. nonbous) then
426 merge = ( (dspv_dt(k2)*(tc(k2-1)-tc(k2)) + dspv_ds(k2)*(sc(k2-1)-sc(k2))) * &
427 ((hc(k2) * hc(k2-1))*i_htot) < abs(tol_merge * dspvxh_sum) )
428 elseif (better_est) then
429 merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
430 ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
431 elseif (nonbous) then
432 merge = ( (dspv_dt(k2)*(tc(k2-1)-tc(k2)) + dspv_ds(k2)*(sc(k2-1)-sc(k2))) * &
433 (hc(k2) + hc(k2-1)) < abs(tol_merge * dspvxh_sum) )
434 else
435 merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
436 (hc(k2) + hc(k2-1)) < tol_merge*drxh_sum)
437 endif
438 if (merge) then
439 ! Merge the two bottommost layers. At this point kc = k2.
440 i_hnew = 1.0 / (hc(kc) + hc(kc-1))
441 tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
442 sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
443 hc(kc-1) = (hc(kc) + hc(kc-1))
444 kc = kc - 1
445 else ; exit ; endif
446 enddo
447 else
448 ! Add a new layer to the column.
449 kc = kc + 1
450 if (nonbous) then
451 dspv_ds(kc) = dspv_ds(k) ; dspv_dt(kc) = dspv_dt(k)
452 else
453 drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
454 endif
455 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
456 endif
457 enddo
458 ! At this point there are kc layers and the gprimes should be positive.
459 if (nonbous) then
460 do k=2,kc
461 gprime(k) = h_to_pres * (dspv_dt(k)*(tc(k-1)-tc(k)) + dspv_ds(k)*(sc(k-1)-sc(k)))
462 enddo
463 else
464 do k=2,kc
465 gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + drho_ds(k)*(sc(k)-sc(k-1)))
466 enddo
467 endif
468 else ! .not. (use_EOS)
469 ! Do the same with density directly...
470 kc = 1
471 hc(1) = hf(1,i) ; rc(1) = rf(1,i)
472 do k=2,kf(i)
473 if (nonbous .and. better_est) then
474 merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < &
475 (rc(kc)*rf(k,i)) * abs(2.0 * tol_merge * dspvxh_sum))
476 elseif (nonbous) then
477 merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < &
478 (rc(kc)*rf(k,i)) * abs(2.0 * tol_merge * dspvxh_sum))
479 elseif (better_est) then
480 merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < 2.0*tol_merge*drxh_sum)
481 else
482 merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol_merge*drxh_sum)
483 endif
484 if (merge) then
485 ! Merge this layer with the one above and backtrack.
486 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
487 hc(kc) = (hc(kc) + hf(k,i))
488 ! Backtrack to remove any convective instabilities above... Note
489 ! that the tolerance is a factor of two larger, to avoid limit how
490 ! far back we go.
491 do k2=kc,2,-1
492 if (nonbous .and. better_est) then
493 merge = ((rc(k2) - rc(k2-1)) * ((hc(kc) * hf(k,i))*i_htot) < &
494 (rc(k2-1)*rc(k2)) * abs(2.0 * tol_merge * dspvxh_sum))
495 elseif (nonbous) then
496 merge = ((rc(k2) - rc(k2-1)) * (hc(kc) + hf(k,i)) < &
497 (rc(k2-1)*rc(k2)) * abs(2.0 * tol_merge * dspvxh_sum))
498 elseif (better_est) then
499 merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
500 else
501 merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol_merge*drxh_sum)
502 endif
503 if (merge) then
504 ! Merge the two bottommost layers. At this point kc = k2.
505 rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
506 hc(kc-1) = (hc(kc) + hc(kc-1))
507 kc = kc - 1
508 else ; exit ; endif
509 enddo
510 else
511 ! Add a new layer to the column.
512 kc = kc + 1
513 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
514 endif
515 enddo
516 ! At this point there are kc layers and the gprimes should be positive.
517 if (nonbous) then
518 do k=2,kc
519 gprime(k) = h_to_pres * (rc(k) - rc(k-1)) / (rc(k) * rc(k-1))
520 enddo
521 else
522 do k=2,kc
523 gprime(k) = g_rho0 * (rc(k)-rc(k-1))
524 enddo
525 endif
526 endif ! use_EOS
527
528 ! Sum the contributions from all of the interfaces to give an over-estimate
529 ! of the first-mode wave speed. Also populate Igl and Igu which are the
530 ! non-leading diagonals of the tridiagonal matrix.
531 if (kc >= 2) then
532 speed2_tot = 0.0
533 if (better_est) then
534 h_top(1) = 0.0 ; h_bot(kc+1) = 0.0
535 do k=2,kc+1 ; h_top(k) = h_top(k-1) + hc(k-1) ; enddo
536 do k=kc,2,-1 ; h_bot(k) = h_bot(k+1) + hc(k) ; enddo
537 i_htot = 0.0 ; if (h_top(kc+1) > 0.0) i_htot = 1.0 / h_top(kc+1)
538 endif
539
540 if (l_use_ebt_mode) then
541 igu(1) = 0. ! Neumann condition for pressure modes
542 sum_hc = hc(1)
543 n2min = gprime(2)/hc(1)
544
545 below_mono_n2_frac = .false.
546 below_mono_n2_depth = .false.
547 do k=2,kc
548 hw = 0.5*(hc(k-1)+hc(k))
549 gp = gprime(k)
550
551 if (l_mono_n2_column_fraction>0. .or. l_mono_n2_depth>=0.) then
552 ! Determine whether N2 estimates should not be allowed to increase with depth.
553 if (l_mono_n2_column_fraction>0.) then
554 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
555 below_mono_n2_frac = &
556 (max(g%meanSL(i,j) + g%bathyT(i,j), 0.0) - gv%H_to_Z * sum_hc < &
557 l_mono_n2_column_fraction * max(g%meanSL(i,j) + g%bathyT(i,j), 0.0))
558 else
559 below_mono_n2_frac = (htot(i) - sum_hc < l_mono_n2_column_fraction*htot(i))
560 endif
561 endif
562 if (l_mono_n2_depth >= 0.) below_mono_n2_depth = (sum_hc > l_mono_n2_depth)
563
564 if ( (gp > n2min*hw) .and. (below_mono_n2_frac .or. below_mono_n2_depth) ) then
565 ! Filters out regions where N2 increases with depth, but only in a lower fraction
566 ! of the water column or below a certain depth.
567 gp = n2min * hw
568 else
569 n2min = gp / hw
570 endif
571 endif
572
573 igu(k) = 1.0/(gp*hc(k))
574 igl(k-1) = 1.0/(gp*hc(k-1))
575 sum_hc = sum_hc + hc(k)
576
577 if (better_est) then
578 ! Estimate that the ebt_mode is sqrt(2) times the speed of the flat bottom modes.
579 speed2_tot = speed2_tot + 2.0 * gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
580 else ! The ebt_mode wave should be faster than the flat-bottom mode, so 0.707 should be > 1?
581 speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))*0.707
582 endif
583 enddo
584 !Igl(kc) = 0. ! Neumann condition for pressure modes
585 igl(kc) = 2.*igu(kc) ! Dirichlet condition for pressure modes
586 else ! .not. l_use_ebt_mode
587 do k=2,kc
588 igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
589 if (better_est) then
590 speed2_tot = speed2_tot + gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
591 else
592 speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
593 endif
594 enddo
595 endif
596
597 if (calc_modal_structure) then
598 mode_struct(:) = 0.
599 mode_struct(1:kc) = 1. ! Uniform flow, first guess
600 endif
601
602 ! Under estimate the first eigenvalue (overestimate the speed) to start with.
603 if (calc_modal_structure) then
604 lam0 = 0.5 / speed2_tot ; lam = lam0
605 else
606 lam0 = 1.0 / speed2_tot ; lam = lam0
607 endif
608 ! Find the determinant and its derivative with lam.
609 do itt=1,max_itt
610 lam_it(itt) = lam
611 if (l_use_ebt_mode) then
612 ! This initialization of det,ddet imply Neumann boundary conditions for horizontal
613 ! velocity or pressure modes, so that first 3 rows of the matrix are
614 ! / b(1)-lam igl(1) 0 0 0 ... \
615 ! | igu(2) b(2)-lam igl(2) 0 0 ... |
616 ! | 0 igu(3) b(3)-lam igl(3) 0 ... |
617 ! The last two rows of the pressure equation matrix are
618 ! | ... 0 igu(kc-1) b(kc-1)-lam igl(kc-1) |
619 ! \ ... 0 0 igu(kc) b(kc)-lam /
620 call tridiag_det(igu, igl, 1, kc, lam, det, ddet, row_scale=c2_scale)
621 else
622 ! This initialization of det,ddet imply Dirichlet boundary conditions for vertical
623 ! velocity modes, so that first 3 rows of the matrix are
624 ! / b(2)-lam igl(2) 0 0 0 ... |
625 ! | igu(3) b(3)-lam igl(3) 0 0 ... |
626 ! | 0 igu(4) b(4)-lam igl(4) 0 ... |
627 ! The last three rows of the w equation matrix are
628 ! | ... 0 igu(kc-2) b(kc-2)-lam igl(kc-2) 0 |
629 ! | ... 0 0 igu(kc-1) b(kc-1)-lam igl(kc-1) |
630 ! \ ... 0 0 0 igu(kc) b(kc)-lam /
631 call tridiag_det(igu, igl, 2, kc, lam, det, ddet, row_scale=c2_scale)
632 endif
633 ! Use Newton's method iteration to find a new estimate of lam.
634 det_it(itt) = det ; ddet_it(itt) = ddet
635
636 if ((ddet >= 0.0) .or. (-det > -0.5*lam*ddet)) then
637 ! lam was not an under-estimate, as intended, so Newton's method
638 ! may not be reliable; lam must be reduced, but not by more
639 ! than half.
640 lam = 0.5 * lam
641 dlam = -lam
642 else ! Newton's method is OK.
643 dlam = - det / ddet
644 lam = lam + dlam
645 endif
646
647 if (calc_modal_structure) then
648 call tdma6(kc, igu, igl, lam, mode_struct)
649 ! Note that tdma6 changes the units of mode_struct to [L2 T-2 ~> m2 s-2]
650 ms_min = mode_struct(1)
651 ms_max = mode_struct(1)
652 ms_sq = mode_struct(1)**2
653 do k = 2,kc
654 ms_min = min(ms_min, mode_struct(k))
655 ms_max = max(ms_max, mode_struct(k))
656 ms_sq = ms_sq + mode_struct(k)**2
657 enddo
658 if (ms_min<0. .and. ms_max>0.) then ! Any zero crossings => lam is too high
659 lam = 0.5 * ( lam - dlam )
660 dlam = -lam
661 mode_struct(1:kc) = abs(mode_struct(1:kc)) / sqrt( ms_sq )
662 else
663 mode_struct(1:kc) = mode_struct(1:kc) / sqrt( ms_sq )
664 endif
665 ! After the nondimensionalization above, mode_struct is once again [nondim]
666 endif
667
668 if (abs(dlam) < tol_solve*lam) exit
669 enddo
670
671 cg1(i,j) = 0.0
672 if (lam > 0.0) cg1(i,j) = 1.0 / sqrt(lam)
673
674 if (present(modal_structure)) then
675 if (mode_struct(1)/=0.) then ! Normalize
676 mode_struct(1:kc) = mode_struct(1:kc) / mode_struct(1)
677 else
678 mode_struct(1:kc)=0.
679 endif
680
681 if (cs%remap_answer_date < 20190101) then
682 call remapping_core_h(cs%remap_2018_CS, kc, hc(:), mode_struct, &
683 nz, h(i,j,:), modal_structure(i,j,:))
684 else
685 call remapping_core_h(cs%remap_CS, kc, hc(:), mode_struct, &
686 nz, h(i,j,:), modal_structure(i,j,:))
687 endif
688 endif
689 else
690 cg1(i,j) = 0.0
691 if (present(modal_structure)) modal_structure(i,j,:) = 0.
692 endif
693 endif ! cg1 /= 0.0
694 else
695 cg1(i,j) = 0.0 ! This is a land point.
696 if (present(modal_structure)) modal_structure(i,j,:) = 0.
697 endif ; enddo ! i-loop
698 enddo ! j-loop
699
700end subroutine wave_speed
701
702!> Solve a non-symmetric tridiagonal problem with the sum of the upper and lower diagonals minus a
703!! scalar contribution as the leading diagonal.
704!! This uses the Thomas algorithm rather than the Hallberg algorithm since the matrix is not symmetric.
705subroutine tdma6(n, a, c, lam, y)
706 integer, intent(in) :: n !< Number of rows of matrix
707 real, dimension(:), intent(in) :: a !< Lower diagonal [T2 L-2 ~> s2 m-2]
708 real, dimension(:), intent(in) :: c !< Upper diagonal [T2 L-2 ~> s2 m-2]
709 real, intent(in) :: lam !< Scalar subtracted from leading diagonal [T2 L-2 ~> s2 m-2]
710 real, dimension(:), intent(inout) :: y !< RHS on entry [A ~> a], result on exit [A L2 T-2 ~> a m2 s-2]
711
712 ! Local variables
713 real :: lambda ! A temporary variable in [T2 L-2 ~> s2 m-2]
714 real :: beta(n) ! A temporary variable in [T2 L-2 ~> s2 m-2]
715 real :: I_beta(n) ! A temporary variable in [L2 T-2 ~> m2 s-2]
716 real :: yy(n) ! A temporary variable with the same units as y on entry [A ~> a]
717 integer :: k, m
718
719 lambda = lam
720 beta(1) = (a(1)+c(1)) - lambda
721 if (beta(1)==0.) then ! lam was chosen too perfectly
722 ! Change lambda and redo this first row
723 lambda = (1. + 1.e-5) * lambda
724 beta(1) = (a(1)+c(1)) - lambda
725 endif
726 i_beta(1) = 1. / beta(1)
727 yy(1) = y(1)
728 do k = 2, n
729 beta(k) = ( (a(k)+c(k)) - lambda ) - a(k) * c(k-1) * i_beta(k-1)
730 ! Perhaps the following 0 needs to become a tolerance to handle underflow?
731 if (beta(k)==0.) then ! lam was chosen too perfectly
732 ! Change lambda and redo everything up to row k
733 lambda = (1. + 1.e-5) * lambda
734 i_beta(1) = 1. / ( (a(1)+c(1)) - lambda )
735 do m = 2, k
736 i_beta(m) = 1. / ( ( (a(m)+c(m)) - lambda ) - a(m) * c(m-1) * i_beta(m-1) )
737 yy(m) = y(m) + a(m) * yy(m-1) * i_beta(m-1)
738 enddo
739 else
740 i_beta(k) = 1. / beta(k)
741 endif
742 yy(k) = y(k) + a(k) * yy(k-1) * i_beta(k-1)
743 enddo
744 ! The units of y change by a factor of [L2 T-2 ~> m2 s-2] in the following lines.
745 y(n) = yy(n) * i_beta(n)
746 do k = n-1, 1, -1
747 y(k) = ( yy(k) + c(k) * y(k+1) ) * i_beta(k)
748 enddo
749
750end subroutine tdma6
751
752!> Calculates the wave speeds for the first few barolinic modes.
753subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_struct_max, u_struct_bot, Nb, int_w2, &
754 int_U2, int_N2w2, halo_size)
755 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
756 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
757 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
758 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
759 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
760 integer, intent(in) :: nmodes !< Number of modes
761 type(wave_speed_cs), intent(in) :: cs !< Wave speed control struct
762 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1,nmodes),intent(out) :: w_struct !< Wave vertical velocity profile [nondim]
763 real, dimension(SZI_(G),SZJ_(G),SZK_(GV),nmodes),intent(out) :: u_struct !< Wave horizontal velocity profile
764 !! [Z-1 ~> m-1]
765 real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: cn !< Waves speeds [L T-1 ~> m s-1]
766 real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: u_struct_max !< Maximum of wave horizontal velocity
767 !! profile [Z-1 ~> m-1]
768 real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: u_struct_bot !< Bottom value of wave horizontal
769 !! velocity profile [Z-1 ~> m-1]
770 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: nb !< Bottom value of buoyancy freqency
771 !! [T-1 ~> s-1]
772 real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_w2 !< depth-integrated vertical velocity
773 !! profile squared [H ~> m or kg m-2]
774 real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_u2 !< depth-integrated horizontal velocity
775 !! profile squared [H Z-2 ~> m-1 or kg m-4]
776 real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_n2w2 !< depth-integrated buoyancy frequency
777 !! times vertical velocity profile
778 !! squared [H T-2 ~> m s-2 or kg m-2 s-2]
779 integer, optional, intent(in) :: halo_size !< Width of halo within which to
780 !! calculate wave speeds
781
782 ! Local variables
783 real, dimension(SZK_(GV)+1) :: &
784 drho_dt, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
785 drho_ds, & ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
786 dspv_dt, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
787 dspv_ds, & ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
788 pres, & ! Interface pressure [R L2 T-2 ~> Pa]
789 t_int, & ! Temperature interpolated to interfaces [C ~> degC]
790 s_int, & ! Salinity interpolated to interfaces [S ~> ppt]
791 h_top, & ! The distance of each filtered interface from the ocean surface [H ~> m or kg m-2]
792 h_bot, & ! The distance of each filtered interface from the bottom [H ~> m or kg m-2]
793 gprime, & ! The reduced gravity across each interface [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
794 n2 ! The buoyancy freqency squared [T-2 ~> s-2]
795 real, dimension(SZK_(GV),SZI_(G)) :: &
796 hf, & ! Layer thicknesses after very thin layers are combined [H ~> m or kg m-2]
797 dzf, & ! Layer vertical extents after very thin layers are combined [Z ~> m]
798 tf, & ! Layer temperatures after very thin layers are combined [C ~> degC]
799 sf, & ! Layer salinities after very thin layers are combined [S ~> ppt]
800 rf ! Layer densities after very thin layers are combined [R ~> kg m-3]
801 real, dimension(SZI_(G),SZK_(GV)) :: &
802 dz_2d ! Height change across layers [Z ~> m]
803 real, dimension(SZK_(GV)) :: &
804 igl, igu, & ! The inverse of the reduced gravity across an interface times
805 ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].
806 hc, & ! A column of layer thicknesses after convective instabilities are removed [H ~> m or kg m-2]
807 dzc, & ! A column of layer vertical extents after convective instabilities are removed [Z ~> m]
808 tc, & ! A column of layer temperatures after convective instabilities are removed [C ~> degC]
809 sc, & ! A column of layer salinities after convective instabilities are removed [S ~> ppt]
810 rc ! A column of layer densities after convective instabilities are removed [R ~> kg m-3]
811 real :: i_htot ! The inverse of the total filtered thicknesses [H-1 ~> m-1 or m2 kg-1]
812 real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant and its
813 ! derivative with lam between rows of the Thomas algorithm solver [L2 s2 T-2 m-2 ~> nondim].
814 ! The exact value should not matter for the final result if it is an even power of 2.
815 real :: det, ddet ! Determinant of the eigen system and its derivative with lam. Because the
816 ! units of the eigenvalue change with the number of layers and because of the
817 ! dynamic rescaling that is used to keep det in a numerically representable range,
818 ! the units of of det are hard to interpret, but det/ddet is always in units
819 ! of [T2 L-2 ~> s2 m-2]
820 real :: lam_1 ! approximate mode-1 eigenvalue [T2 L-2 ~> s2 m-2]
821 real :: lam_n ! approximate mode-n eigenvalue [T2 L-2 ~> s2 m-2]
822 real :: dlam ! The change in estimates of the eigenvalue [T2 L-2 ~> s2 m-2]
823 real :: lammin ! minimum lam value for root searching range [T2 L-2 ~> s2 m-2]
824 real :: lammax ! maximum lam value for root searching range [T2 L-2 ~> s2 m-2]
825 real :: laminc ! width of moving window for root searching [T2 L-2 ~> s2 m-2]
826 real :: det_l, ddet_l ! determinant of the eigensystem and its derivative with lam at the lower
827 ! end of the range of values bracketing a particular root, in dynamically
828 ! rescaled units that may differ from the other det variables, but such
829 ! that the units of det_l/ddet_l are [T2 L-2 ~> s2 m-2]
830 real :: det_r, ddet_r ! determinant and its derivative with lam at the lower end of the
831 ! bracket in arbitrarily rescaled units, but such that the units of
832 ! det_r/ddet_r are [T2 L-2 ~> s2 m-2]
833 real :: det_sub, ddet_sub ! determinant and its derivative with lam at a subinterval endpoint that
834 ! is a candidate for a new bracket endpoint in arbitrarily rescaled units,
835 ! but such that the units of det_sub/ddet_sub are [T2 L-2 ~> s2 m-2]
836 real :: xl, xr ! lam guesses at left and right of window [T2 L-2 ~> s2 m-2]
837 real :: xl_sub ! lam guess at left of subinterval window [T2 L-2 ~> s2 m-2]
838 real, dimension(nmodes) :: &
839 xbl, xbr ! lam guesses bracketing a zero-crossing (root) [T2 L-2 ~> s2 m-2]
840 integer :: numint ! number of widows (intervals) in root searching range
841 integer :: nrootsfound ! number of extra roots found (not including 1st root)
842 real :: h_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
843 real, dimension(SZI_(G)) :: &
844 htot, hmin, & ! Thicknesses [H ~> m or kg m-2]
845 h_here, & ! A layer thickness [H ~> m or kg m-2]
846 dz_here, & ! A layer vertical extent [Z ~> m]
847 hxt_here, & ! A layer integrated temperature [C H ~> degC m or degC kg m-2]
848 hxs_here, & ! A layer integrated salinity [S H ~> ppt m or ppt kg m-2]
849 hxr_here ! A layer integrated density [R H ~> kg m-2 or kg2 m-5]
850 real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
851 real :: speed2_min ! minimum mode speed (squared) to consider in root searching [L2 T-2 ~> m2 s-2]
852 real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2]
853 real :: cg1_est ! An initial estimate of the squared first mode speed [L2 T-2 ~> m2 s-2]
854 real, parameter :: reduct_factor = 0.5 ! A factor used in setting speed2_min [nondim]
855 real :: i_hnew ! The inverse of a new layer thickness [H-1 ~> m-1 or m2 kg-1]
856 real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R H ~> kg m-2 or kg2 m-5]
857 real :: dspvxh_sum ! The sum of specific volume differences across interfaces times
858 ! thicknesses [H R-1 ~> m4 kg-1 or m], negative for stable stratification.
859 real :: g_rho0 ! G_Earth/Rho0 [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
860 real :: tol_hfrac ! Layers that together are smaller than this fraction of
861 ! the total water column can be merged for efficiency [nondim].
862 real :: min_h_frac ! tol_Hfrac divided by the total number of layers [nondim].
863 real :: tol_solve ! The fractional tolerance with which to solve for the wave speeds [nondim].
864 real :: tol_merge ! The fractional change in estimated wave speed that is allowed
865 ! when deciding to merge layers in the calculation [nondim]
866 integer :: kf(szi_(g)) ! The number of active layers after filtering.
867 integer, parameter :: max_itt = 30
868 logical :: use_eos ! If true, density or specific volume is calculated from T & S using the equation of state.
869 logical :: nonbous ! If true, do not make the Boussinesq approximation.
870 logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed.
871 logical :: merge ! If true, merge the current layer with the one above.
872 integer :: nsub ! number of subintervals used for root finding
873 integer, parameter :: sub_it_max = 4
874 ! maximum number of times to subdivide interval
875 ! for root finding (# intervals = 2**sub_it_max)
876 logical :: sub_rootfound ! if true, subdivision has located root
877 integer :: kc ! The number of layers in the column after merging
878 integer :: sub, sub_it
879 integer :: i, j, k, k2, itt, is, ie, js, je, nz, iint, m, halo
880 real, dimension(SZK_(GV)+1) :: modal_structure !< Normalized model structure [nondim]
881 real, dimension(SZK_(GV)) :: modal_structure_fder !< Normalized model structure [Z-1 ~> m-1]
882 real :: mode_struct(szk_(gv)+1) ! The mode structure [nondim], but it is also temporarily
883 ! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6.
884 real :: mode_struct_fder(szk_(gv)) ! The mode structure 1st derivative [Z-1 ~> m-1], but it is also temporarily
885 ! in units of [L2 Z-1 T-2 ~> m s-2] after it is modified inside of tdma6.
886 real :: mode_struct_sq(szk_(gv)+1) ! The square of mode structure [nondim]
887 real :: mode_struct_fder_sq(szk_(gv)) ! The square of mode structure 1st derivative [Z-2 ~> m-2]
888
889 real :: w2avg ! A total for renormalization [H L4 T-4 ~> m5 s-4 or kg m2 s-4]
890 real, parameter :: a_int = 0.5 ! Integral total for normalization [nondim]
891 real :: renorm ! Normalization factor [T2 L-2 ~> s2 m-2]
892
893 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke ; halo = 0
894
895 if (.not. cs%initialized) call mom_error(fatal, "MOM_wave_speed / wave_speeds: "// &
896 "Module must be initialized before it is used.")
897
898 if (present(halo_size)) then
899 halo = halo_size
900 is = g%isc - halo ; ie = g%iec + halo ; js = g%jsc - halo ; je = g%jec + halo
901 endif
902
903 nonbous = .not.(gv%Boussinesq .or. gv%semi_Boussinesq)
904 h_to_pres = gv%H_to_RZ * gv%g_Earth
905 if (.not.nonbous) g_rho0 = gv%g_Earth * gv%H_to_Z / gv%Rho0
906 use_eos = associated(tv%eqn_of_state)
907
908 if (cs%c1_thresh < 0.0) &
909 call mom_error(fatal, "INTERNAL_WAVE_CG1_THRESH must be set to a non-negative "//&
910 "value via wave_speed_init for wave_speeds to be used.")
911 c2_scale = us%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.
912
913 better_est = cs%better_cg1_est
914 if (better_est) then
915 tol_solve = cs%wave_speed_tol
916 tol_hfrac = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
917 else
918 tol_solve = 0.001 ; tol_hfrac = 0.0001 ; tol_merge = 0.001
919 endif
920 cg1_min2 = cs%min_speed2
921
922 ! Zero out all local values. Values over land or for columns that are too weakly stratified
923 ! are not changed from this zero value.
924 cn(:,:,:) = 0.0
925 u_struct_max(:,:,:) = 0.0
926 u_struct_bot(:,:,:) = 0.0
927 nb(:,:) = 0.0
928 int_w2(:,:,:) = 0.0
929 int_n2w2(:,:,:) = 0.0
930 int_u2(:,:,:) = 0.0
931 u_struct(:,:,:,:) = 0.0
932 w_struct(:,:,:,:) = 0.0
933
934 min_h_frac = tol_hfrac / real(nz)
935 !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,CS,use_EOS,nonBous, &
936 !$OMP min_h_frac,H_to_pres,tv,cn,g_Rho0,nmodes,cg1_min2, &
937 !$OMP better_est,tol_solve,tol_merge,c2_scale)
938 do j=js,je
939 ! First merge very thin layers with the one above (or below if they are
940 ! at the top). This also transposes the row order so that columns can
941 ! be worked upon one at a time.
942 do i=is,ie ; htot(i) = 0.0 ; enddo
943 do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
944
945 call thickness_to_dz(h, tv, dz_2d, j, g, gv, halo_size=halo)
946
947 do i=is,ie
948 hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0 ; dz_here(i) = 0.0
949 hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
950 enddo
951 if (use_eos) then
952 do k=1,nz ; do i=is,ie
953 if ((h_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then
954 hf(kf(i),i) = h_here(i)
955 dzf(kf(i),i) = dz_here(i)
956 tf(kf(i),i) = hxt_here(i) / h_here(i)
957 sf(kf(i),i) = hxs_here(i) / h_here(i)
958 kf(i) = kf(i) + 1
959
960 ! Start a new layer
961 h_here(i) = h(i,j,k)
962 dz_here(i) = dz_2d(i,k)
963 hxt_here(i) = h(i,j,k)*tv%T(i,j,k)
964 hxs_here(i) = h(i,j,k)*tv%S(i,j,k)
965 else
966 h_here(i) = h_here(i) + h(i,j,k)
967 dz_here(i) = dz_here(i) + dz_2d(i,k)
968 hxt_here(i) = hxt_here(i) + h(i,j,k)*tv%T(i,j,k)
969 hxs_here(i) = hxs_here(i) + h(i,j,k)*tv%S(i,j,k)
970 endif
971 enddo ; enddo
972 do i=is,ie ; if (h_here(i) > 0.0) then
973 hf(kf(i),i) = h_here(i)
974 dzf(kf(i),i) = dz_here(i)
975 tf(kf(i),i) = hxt_here(i) / h_here(i)
976 sf(kf(i),i) = hxs_here(i) / h_here(i)
977 endif ; enddo
978 else ! .not. (use_EOS)
979 do k=1,nz ; do i=is,ie
980 if ((h_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then
981 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
982 dzf(kf(i),i) = dz_here(i)
983 kf(i) = kf(i) + 1
984
985 ! Start a new layer
986 h_here(i) = h(i,j,k)
987 dz_here(i) = dz_2d(i,k)
988 hxr_here(i) = h(i,j,k)*gv%Rlay(k)
989 else
990 h_here(i) = h_here(i) + h(i,j,k)
991 dz_here(i) = dz_here(i) + dz_2d(i,k)
992 hxr_here(i) = hxr_here(i) + h(i,j,k)*gv%Rlay(k)
993 endif
994 enddo ; enddo
995 do i=is,ie ; if (h_here(i) > 0.0) then
996 hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
997 dzf(kf(i),i) = dz_here(i)
998 endif ; enddo
999 endif
1000
1001 ! From this point, we can work on individual columns without causing memory to have page faults.
1002 do i=is,ie
1003 if (g%mask2dT(i,j) > 0.0) then
1004 if (use_eos) then
1005 pres(1) = 0.0 ; h_top(1) = 0.0
1006 do k=2,kf(i)
1007 pres(k) = pres(k-1) + h_to_pres*hf(k-1,i)
1008 t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
1009 s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
1010 h_top(k) = h_top(k-1) + hf(k-1,i)
1011 enddo
1012 if (nonbous) then
1013 call calculate_specific_vol_derivs(t_int, s_int, pres, dspv_dt, dspv_ds, &
1014 tv%eqn_of_state, (/2,kf(i)/) )
1015 else
1016 call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, &
1017 tv%eqn_of_state, (/2,kf(i)/) )
1018 endif
1019
1020 ! Sum the reduced gravities to find out how small a density difference is negligibly small.
1021 drxh_sum = 0.0 ; dspvxh_sum = 0.0
1022 if (better_est) then
1023 ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for
1024 ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.
1025 ! For a uniform stratification and a huge number of layers uniformly distributed in
1026 ! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.
1027 if (h_top(kf(i)) > 0.0) then
1028 i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
1029 h_bot(kf(i)+1) = 0.0
1030 if (nonbous) then
1031 do k=kf(i),2,-1
1032 h_bot(k) = h_bot(k+1) + hf(k,i)
1033 dspvxh_sum = dspvxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
1034 min(0.0, dspv_dt(k)*(tf(k,i)-tf(k-1,i)) + dspv_ds(k)*(sf(k,i)-sf(k-1,i)))
1035 enddo
1036 else
1037 do k=kf(i),2,-1
1038 h_bot(k) = h_bot(k+1) + hf(k,i)
1039 drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
1040 max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
1041 enddo
1042 endif
1043 endif
1044 else
1045 ! This estimate is problematic in that it goes like 1/nz for a large number of layers,
1046 ! but it is an overestimate (as desired) for a small number of layers, by at a factor
1047 ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.
1048 if (nonbous) then
1049 do k=2,kf(i)
1050 dspvxh_sum = dspvxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
1051 min(0.0, dspv_dt(k)*(tf(k,i)-tf(k-1,i)) + dspv_ds(k)*(sf(k,i)-sf(k-1,i)))
1052 enddo
1053 else
1054 do k=2,kf(i)
1055 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
1056 max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
1057 enddo
1058 endif
1059 endif
1060 else ! Not use_EOS
1061 drxh_sum = 0.0 ; dspvxh_sum = 0.0
1062 if (better_est) then
1063 h_top(1) = 0.0
1064 do k=2,kf(i) ; h_top(k) = h_top(k-1) + hf(k-1,i) ; enddo
1065 if (h_top(kf(i)) > 0.0) then
1066 i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
1067 h_bot(kf(i)+1) = 0.0
1068 if (nonbous) then
1069 do k=kf(i),2,-1
1070 h_bot(k) = h_bot(k+1) + hf(k,i)
1071 dspvxh_sum = dspvxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
1072 min(0.0, (rf(k-1,i)-rf(k,i)) / (rf(k,i)*rf(k-1,i)))
1073 enddo
1074 else
1075 do k=kf(i),2,-1
1076 h_bot(k) = h_bot(k+1) + hf(k,i)
1077 drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * max(0.0,rf(k,i)-rf(k-1,i))
1078 enddo
1079 endif
1080 endif
1081 else
1082 do k=2,kf(i)
1083 drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
1084 enddo
1085 endif
1086 endif
1087
1088 if (nonbous) then
1089 ! Note that dSpVxh_sum is negative for stable stratification.
1090 cg1_est = h_to_pres * abs(dspvxh_sum)
1091 else
1092 cg1_est = g_rho0 * drxh_sum
1093 endif
1094
1095 ! Find gprime across each internal interface, taking care of convective
1096 ! instabilities by merging layers.
1097 if (cg1_est > cg1_min2) then
1098 ! Merge layers to eliminate convective instabilities or exceedingly
1099 ! small reduced gravities. Merging layers reduces the estimated wave speed by
1100 ! (rho(2)-rho(1))*h(1)*h(2) / H_tot.
1101 if (use_eos) then
1102 kc = 1
1103 hc(1) = hf(1,i) ; dzc(1) = dzf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
1104 do k=2,kf(i)
1105 if (better_est .and. nonbous) then
1106 merge = ((dspv_dt(k)*(tc(kc)-tf(k,i)) + dspv_ds(k)*(sc(kc)-sf(k,i))) * &
1107 ((hc(kc) * hf(k,i))*i_htot) < abs(2.0 * tol_merge * dspvxh_sum))
1108 elseif (better_est) then
1109 merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
1110 ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
1111 elseif (nonbous) then
1112 merge = ((dspv_dt(k)*(tc(kc)-tf(k,i)) + dspv_ds(k)*(sc(kc)-sf(k,i))) * &
1113 (hc(kc) + hf(k,i)) < abs(2.0 * tol_merge * dspvxh_sum))
1114 else
1115 merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
1116 (hc(kc) + hf(k,i)) < 2.0 * tol_merge*drxh_sum)
1117 endif
1118 if (merge) then
1119 ! Merge this layer with the one above and backtrack.
1120 i_hnew = 1.0 / (hc(kc) + hf(k,i))
1121 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
1122 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
1123 hc(kc) = hc(kc) + hf(k,i)
1124 dzc(kc) = dzc(kc) + dzf(k,i)
1125 ! Backtrack to remove any convective instabilities above... Note
1126 ! that the tolerance is a factor of two larger, to avoid limit how
1127 ! far back we go.
1128 do k2=kc,2,-1
1129 if (better_est .and. nonbous) then
1130 merge = ( (dspv_dt(k2)*(tc(k2-1)-tc(k2)) + dspv_ds(k2)*(sc(k2-1)-sc(k2))) * &
1131 ((hc(k2) * hc(k2-1))*i_htot) < abs(tol_merge * dspvxh_sum) )
1132 elseif (better_est) then
1133 merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
1134 ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
1135 elseif (nonbous) then
1136 merge = ( (dspv_dt(k2)*(tc(k2-1)-tc(k2)) + dspv_ds(k2)*(sc(k2-1)-sc(k2))) * &
1137 (hc(k2) + hc(k2-1)) < abs(tol_merge * dspvxh_sum) )
1138 else
1139 merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
1140 (hc(k2) + hc(k2-1)) < tol_merge*drxh_sum)
1141 endif
1142 if (merge) then
1143 ! Merge the two bottommost layers. At this point kc = k2.
1144 i_hnew = 1.0 / (hc(kc) + hc(kc-1))
1145 tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
1146 sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
1147 hc(kc-1) = hc(kc) + hc(kc-1)
1148 dzc(kc-1) = dzc(kc) + dzc(kc-1)
1149 kc = kc - 1
1150 else ; exit ; endif
1151 enddo
1152 else
1153 ! Add a new layer to the column.
1154 kc = kc + 1
1155 if (nonbous) then
1156 dspv_ds(kc) = dspv_ds(k) ; dspv_dt(kc) = dspv_dt(k)
1157 else
1158 drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
1159 endif
1160 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i) ; dzc(kc) = dzf(k,i)
1161 endif
1162 enddo
1163 ! At this point there are kc layers and the gprimes should be positive.
1164 if (nonbous) then
1165 do k=2,kc
1166 gprime(k) = h_to_pres * (dspv_dt(k)*(tc(k-1)-tc(k)) + dspv_ds(k)*(sc(k-1)-sc(k)))
1167 enddo
1168 else
1169 do k=2,kc
1170 gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + drho_ds(k)*(sc(k)-sc(k-1)))
1171 enddo
1172 endif
1173 else ! .not. (use_EOS)
1174 ! Do the same with density directly...
1175 kc = 1
1176 hc(1) = hf(1,i) ; dzc(1) = dzf(1,i) ; rc(1) = rf(1,i)
1177 do k=2,kf(i)
1178 if (nonbous .and. better_est) then
1179 merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < &
1180 (rc(kc)*rf(k,i)) * abs(2.0 * tol_merge * dspvxh_sum))
1181 elseif (nonbous) then
1182 merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < &
1183 (rc(kc)*rf(k,i)) * abs(2.0 * tol_merge * dspvxh_sum))
1184 elseif (better_est) then
1185 merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < 2.0*tol_merge*drxh_sum)
1186 else
1187 merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol_merge*drxh_sum)
1188 endif
1189 if (merge) then
1190 ! Merge this layer with the one above and backtrack.
1191 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
1192 hc(kc) = hc(kc) + hf(k,i)
1193 dzc(kc) = dzc(kc) + dzf(k,i)
1194 ! Backtrack to remove any convective instabilities above... Note
1195 ! that the tolerance is a factor of two larger, to avoid limit how
1196 ! far back we go.
1197 do k2=kc,2,-1
1198 if (better_est) then
1199 merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
1200 else
1201 merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol_merge*drxh_sum)
1202 endif
1203 if (merge) then
1204 ! Merge the two bottommost layers. At this point kc = k2.
1205 rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
1206 hc(kc-1) = hc(kc) + hc(kc-1)
1207 dzc(kc-1) = dzc(kc) + dzc(kc-1)
1208 kc = kc - 1
1209 else ; exit ; endif
1210 enddo
1211 else
1212 ! Add a new layer to the column.
1213 kc = kc + 1
1214 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i) ; dzc(kc) = dzf(k,i)
1215 endif
1216 enddo
1217 ! At this point there are kc layers and the gprimes should be positive.
1218 if (nonbous) then
1219 do k=2,kc
1220 gprime(k) = h_to_pres * (rc(k) - rc(k-1)) / (rc(k) * rc(k-1))
1221 enddo
1222 else
1223 do k=2,kc
1224 gprime(k) = g_rho0 * (rc(k)-rc(k-1))
1225 enddo
1226 endif
1227 endif ! use_EOS
1228
1229 !-----------------NOW FIND WAVE SPEEDS---------------------------------------
1230 ! ig = i + G%idg_offset ; jg = j + G%jdg_offset
1231 ! Sum the contributions from all of the interfaces to give an over-estimate
1232 ! of the first-mode wave speed. Also populate Igl and Igu which are the
1233 ! non-leading diagonals of the tridiagonal matrix.
1234 if (kc >= 2) then
1235 ! initialize speed2_tot
1236 speed2_tot = 0.0
1237 if (better_est) then
1238 h_top(1) = 0.0 ; h_bot(kc+1) = 0.0
1239 do k=2,kc+1 ; h_top(k) = h_top(k-1) + hc(k-1) ; enddo
1240 do k=kc,2,-1 ; h_bot(k) = h_bot(k+1) + hc(k) ; enddo
1241 i_htot = 0.0 ; if (h_top(kc+1) > 0.0) i_htot = 1.0 / h_top(kc+1)
1242 endif
1243
1244 ! Calculate Igu, Igl, depth, and N2 at each interior interface
1245 ! [excludes surface (K=1) and bottom (K=kc+1)]
1246 igl(:) = 0.
1247 igu(:) = 0.
1248 n2(:) = 0.
1249
1250 do k=2,kc
1251 igl(k) = 1.0 / (gprime(k)*hc(k)) ; igu(k) = 1.0 / (gprime(k)*hc(k-1))
1252 if (nonbous) then
1253 n2(k) = 2.0*us%L_to_Z**2*gprime(k) * (hc(k) + hc(k-1)) / & ! Units are [T-2 ~> s-2]
1254 (dzc(k) + dzc(k-1))**2
1255 else
1256 n2(k) = 2.0*us%L_to_Z**2*gv%Z_to_H*gprime(k) / (dzc(k) + dzc(k-1)) ! Units are [T-2 ~> s-2]
1257 endif
1258 if (better_est) then
1259 speed2_tot = speed2_tot + gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
1260 else
1261 speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
1262 endif
1263 enddo
1264
1265 ! Set stratification for surface and bottom (setting equal to nearest interface for now)
1266 n2(1) = n2(2) ; n2(kc+1) = n2(kc)
1267 ! set bottom stratification
1268 nb(i,j) = sqrt(n2(kc+1))
1269
1270 ! Under estimate the first eigenvalue (overestimate the speed) to start with.
1271 lam_1 = 1.0 / speed2_tot
1272
1273 ! init and first guess for mode structure
1274 mode_struct(:) = 0.
1275 mode_struct_fder(:) = 0.
1276 mode_struct(2:kc) = 1. ! Uniform flow, first guess
1277 modal_structure(:) = 0.
1278 modal_structure_fder(:) = 0.
1279
1280 ! Find the first eigen value
1281 do itt=1,max_itt
1282 ! calculate the determinant of (A-lam_1*I)
1283 call tridiag_det(igu, igl, 2, kc, lam_1, det, ddet, row_scale=c2_scale)
1284
1285 ! If possible, use Newton's method iteration to find a new estimate of lam_1
1286 !det = det_it(itt) ; ddet = ddet_it(itt)
1287 if ((ddet >= 0.0) .or. (-det > -0.5*lam_1*ddet)) then
1288 ! lam_1 was not an under-estimate, as intended, so Newton's method
1289 ! may not be reliable; lam_1 must be reduced, but not by more than half.
1290 lam_1 = 0.5 * lam_1
1291 dlam = -lam_1
1292 else ! Newton's method is OK.
1293 dlam = - det / ddet
1294 lam_1 = lam_1 + dlam
1295 endif
1296
1297 call tdma6(kc-1, igu(2:kc), igl(2:kc), lam_1, mode_struct(2:kc))
1298 ! Note that tdma6 changes the units of mode_struct to [L2 T-2 ~> m2 s-2]
1299 ! apply BC
1300 mode_struct(1) = 0.
1301 mode_struct(kc+1) = 0.
1302
1303 ! renormalization of the integral of the profile
1304 w2avg = 0.0
1305 do k=1,kc
1306 w2avg = w2avg + 0.5*(mode_struct(k)**2+mode_struct(k+1)**2)*hc(k) ! [H L4 T-4 ~> m5 s-4 or kg m2 s-4]
1307 enddo
1308 renorm = sqrt(htot(i)*a_int/w2avg) ! [T2 L-2 ~> s2 m-2]
1309 do k=1,kc+1 ; mode_struct(k) = renorm * mode_struct(k) ; enddo
1310 ! after renorm, mode_struct is again [nondim]
1311 if (abs(dlam) < tol_solve*lam_1) exit
1312 enddo
1313
1314 if (lam_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam_1)
1315
1316 ! sign of wave structure is irrelevant, flip to positive if needed
1317 if (mode_struct(2)<0.) then
1318 mode_struct(2:kc) = -1. * mode_struct(2:kc)
1319 endif
1320
1321 ! vertical derivative of w at interfaces lives on the layer points
1322 do k=1,kc
1323 mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / dzc(k)
1324 enddo
1325
1326 ! boundary condition for derivative is no-gradient
1327 do k=kc+1,nz
1328 mode_struct_fder(k) = mode_struct_fder(kc)
1329 enddo
1330
1331 ! now save maximum value and bottom value
1332 u_struct_bot(i,j,1) = mode_struct_fder(kc)
1333 u_struct_max(i,j,1) = maxval(abs(mode_struct_fder(1:kc)))
1334
1335 ! Calculate terms for vertically integrated energy equation
1336 do k=1,kc
1337 mode_struct_fder_sq(k) = mode_struct_fder(k)**2
1338 enddo
1339 do k=1,kc+1
1340 mode_struct_sq(k) = mode_struct(k)**2
1341 enddo
1342
1343 ! sum over layers for quantities defined on layer
1344 do k=1,kc
1345 int_u2(i,j,1) = int_u2(i,j,1) + mode_struct_fder_sq(k) * hc(k)
1346 enddo
1347
1348 ! vertical integration with Trapezoidal rule for values at interfaces
1349 do k=1,kc
1350 int_w2(i,j,1) = int_w2(i,j,1) + 0.5*(mode_struct_sq(k)+mode_struct_sq(k+1)) * hc(k)
1351 int_n2w2(i,j,1) = int_n2w2(i,j,1) + 0.5*(mode_struct_sq(k)*n2(k) + &
1352 mode_struct_sq(k+1)*n2(k+1)) * hc(k)
1353 enddo
1354
1355 ! for w (diag) interpolate onto all interfaces
1356 call interpolate_column(kc, hc(1:kc), mode_struct(1:kc+1), &
1357 nz, h(i,j,:), modal_structure(:), .false.)
1358
1359 ! for u (remap) onto all layers
1360 call remapping_core_h(cs%remap_CS, kc, hc(1:kc), mode_struct_fder(1:kc), &
1361 nz, h(i,j,:), modal_structure_fder(:))
1362
1363 ! write the wave structure
1364 do k=1,nz+1
1365 w_struct(i,j,k,1) = modal_structure(k)
1366 enddo
1367
1368 do k=1,nz
1369 u_struct(i,j,k,1) = modal_structure_fder(k)
1370 enddo
1371
1372 ! Find other eigen values if c1 is of significant magnitude, > cn_thresh
1373 nrootsfound = 0 ! number of extra roots found (not including 1st root)
1374 if ((nmodes > 1) .and. (kc >= nmodes+1) .and. (cn(i,j,1) > cs%c1_thresh)) then
1375 ! Set the range to look for the other desired eigen values
1376 ! set min value just greater than the 1st root (found above)
1377 lammin = lam_1*(1.0 + tol_solve)
1378 ! set max value based on a low guess at wavespeed for highest mode
1379 speed2_min = (reduct_factor*cn(i,j,1)/real(nmodes))**2
1380 lammax = 1.0 / speed2_min
1381 ! set width of interval (not sure about this - BDM)
1382 laminc = 0.5*lam_1
1383 ! set number of intervals within search range
1384 numint = nint((lammax - lammin)/laminc)
1385
1386 ! Find intervals containing zero-crossings (roots) of the determinant
1387 ! that are beyond the first root
1388
1389 ! find det_l of first interval (det at left endpoint)
1390 call tridiag_det(igu, igl, 2, kc, lammin, det_l, ddet_l, row_scale=c2_scale)
1391 ! move interval window looking for zero-crossings************************
1392 do iint=1,numint
1393 xr = lammin + laminc * iint
1394 xl = xr - laminc
1395 call tridiag_det(igu, igl, 2, kc, xr, det_r, ddet_r, row_scale=c2_scale)
1396 if (det_l*det_r < 0.0) then ! if function changes sign
1397 if (det_l*ddet_l < 0.0) then ! if function at left is headed to zero
1398 nrootsfound = nrootsfound + 1
1399 xbl(nrootsfound) = xl
1400 xbr(nrootsfound) = xr
1401 else
1402 ! function changes sign but has a local max/min in interval,
1403 ! try subdividing interval as many times as necessary (or sub_it_max).
1404 ! loop that increases number of subintervals:
1405 !call MOM_error(WARNING, "determinant changes sign "// &
1406 ! "but has a local max/min in interval; "//&
1407 ! "reduce increment in lam.")
1408 ! begin subdivision loop -------------------------------------------
1409 sub_rootfound = .false. ! initialize
1410 do sub_it=1,sub_it_max
1411 nsub = 2**sub_it ! number of subintervals; nsub=2,4,8,...
1412 ! loop over each subinterval:
1413 do sub=1,nsub-1,2 ! only check odds; sub = 1; 1,3; 1,3,5,7; ...
1414 xl_sub = xl + laminc/(nsub)*sub
1415 call tridiag_det(igu, igl, 2, kc, xl_sub, det_sub, ddet_sub, &
1416 row_scale=c2_scale)
1417 if (det_sub*det_r < 0.0) then ! if function changes sign
1418 if (det_sub*ddet_sub < 0.0) then ! if function at left is headed to zero
1419 sub_rootfound = .true.
1420 nrootsfound = nrootsfound + 1
1421 xbl(nrootsfound) = xl_sub
1422 xbr(nrootsfound) = xr
1423 exit ! exit sub loop
1424 endif ! headed toward zero
1425 endif ! sign change
1426 enddo ! sub-loop
1427 if (sub_rootfound) exit ! root has been found, exit sub_it loop
1428 ! Otherwise, function changes sign but has a local max/min in one of the
1429 ! sub intervals, try subdividing again unless sub_it_max has been reached.
1430 if (sub_it == sub_it_max) then
1431 call mom_error(warning, "wave_speed: root not found "// &
1432 "after sub_it_max subdivisions of original "// &
1433 "interval.")
1434 endif ! sub_it == sub_it_max
1435 enddo ! sub_it-loop-------------------------------------------------
1436 endif ! det_l*ddet_l < 0.0
1437 endif ! det_l*det_r < 0.0
1438 ! exit iint-loop if all desired roots have been found
1439 if (nrootsfound >= nmodes-1) then
1440 ! exit if all additional roots found
1441 exit
1442 elseif (iint == numint) then
1443 ! oops, lamMax not large enough - could add code to increase (BDM)
1444 ! set unfound modes to zero for now (BDM)
1445 ! cn(i,j,nrootsfound+2:nmodes) = 0.0
1446 else
1447 ! else shift interval and keep looking until nmodes or numint is reached
1448 det_l = det_r
1449 ddet_l = ddet_r
1450 endif
1451 enddo ! iint-loop
1452
1453 ! Use Newton's method to find the roots within the identified windows
1454 do m=1,nrootsfound ! loop over the root-containing widows (excluding 1st mode)
1455 lam_n = xbl(m) ! first guess is left edge of window
1456
1457 ! init and first guess for mode structure
1458 mode_struct(:) = 0.
1459 mode_struct_fder(:) = 0.
1460 mode_struct(2:kc) = 1. ! Uniform flow, first guess
1461 modal_structure(:) = 0.
1462 modal_structure_fder(:) = 0.
1463
1464 do itt=1,max_itt
1465 ! calculate the determinant of (A-lam_n*I)
1466 call tridiag_det(igu, igl, 2, kc, lam_n, det, ddet, row_scale=c2_scale)
1467 ! Use Newton's method to find a new estimate of lam_n
1468 dlam = - det / ddet
1469 lam_n = lam_n + dlam
1470
1471 call tdma6(kc-1, igu(2:kc), igl(2:kc), lam_n, mode_struct(2:kc))
1472 ! Note that tdma6 changes the units of mode_struct to [L2 T-2 ~> m2 s-2]
1473 ! apply BC
1474 mode_struct(1) = 0.
1475 mode_struct(kc+1) = 0.
1476
1477 ! renormalization of the integral of the profile
1478 w2avg = 0.0
1479 do k=1,kc
1480 w2avg = w2avg + 0.5*(mode_struct(k)**2+mode_struct(k+1)**2)*hc(k)
1481 enddo
1482 renorm = sqrt(htot(i)*a_int/w2avg)
1483 do k=1,kc+1 ; mode_struct(k) = renorm * mode_struct(k) ; enddo
1484
1485 if (abs(dlam) < tol_solve*lam_1) exit
1486 enddo ! itt-loop
1487
1488 ! calculate nth mode speed
1489 if (lam_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam_n)
1490
1491 ! sign is irrelevant, flip to positive if needed
1492 if (mode_struct(2)<0.) then
1493 mode_struct(2:kc) = -1. * mode_struct(2:kc)
1494 endif
1495
1496 ! derivative of vertical profile (i.e. dw/dz) is evaluated at the layer point
1497 do k=1,kc
1498 mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / dzc(k)
1499 enddo
1500
1501 ! boundary condition for 1st derivative is no-gradient
1502 do k=kc+1,nz
1503 mode_struct_fder(k) = mode_struct_fder(kc)
1504 enddo
1505
1506 ! now save maximum value and bottom value
1507 u_struct_bot(i,j,m) = mode_struct_fder(kc)
1508 u_struct_max(i,j,m) = maxval(abs(mode_struct_fder(1:kc)))
1509
1510 ! Calculate terms for vertically integrated energy equation
1511 do k=1,kc
1512 mode_struct_fder_sq(k) = mode_struct_fder(k)**2
1513 enddo
1514 do k=1,kc+1
1515 mode_struct_sq(k) = mode_struct(k)**2
1516 enddo
1517
1518 ! sum over layers for integral of quantities defined at layer points
1519 do k=1,kc
1520 int_u2(i,j,m) = int_u2(i,j,m) + mode_struct_fder_sq(k) * hc(k)
1521 enddo
1522
1523 ! vertical integration with Trapezoidal rule for quantities on interfaces
1524 do k=1,kc
1525 int_w2(i,j,m) = int_w2(i,j,m) + 0.5*(mode_struct_sq(k)+mode_struct_sq(k+1)) * hc(k)
1526 int_n2w2(i,j,m) = int_n2w2(i,j,m) + 0.5*(mode_struct_sq(k)*n2(k) + &
1527 mode_struct_sq(k+1)*n2(k+1)) * hc(k)
1528 enddo
1529
1530 ! for w (diag) interpolate onto all interfaces
1531 call interpolate_column(kc, hc(1:kc), mode_struct(1:kc+1), &
1532 nz, h(i,j,:), modal_structure(:), .false.)
1533
1534 ! for u (remap) onto all layers
1535 call remapping_core_h(cs%remap_CS, kc, hc(1:kc), mode_struct_fder(1:kc), &
1536 nz, h(i,j,:), modal_structure_fder(:))
1537
1538 ! write the wave structure
1539 ! note that m=1 solves for 2nd mode,...
1540 do k=1,nz+1
1541 w_struct(i,j,k,m+1) = modal_structure(k)
1542 enddo
1543
1544 do k=1,nz
1545 u_struct(i,j,k,m+1) = modal_structure_fder(k)
1546 enddo
1547
1548 enddo ! n-loop
1549 endif ! if nmodes>1 .and. kc>nmodes .and. c1>c1_thresh
1550 endif ! if more than 2 layers
1551 endif ! if drxh_sum < 0
1552 endif ! if not land
1553 enddo ! i-loop
1554 enddo ! j-loop
1555
1556end subroutine wave_speeds
1557
1558!> Calculate the determinant of a tridiagonal matrix with diagonals a,b-lam,c and its derivative
1559!! with lam, where lam is constant across rows. Only the ratio of det to its derivative and their
1560!! signs are typically used, so internal rescaling by consistent factors are used to avoid
1561!! over- or underflow.
1562subroutine tridiag_det(a, c, ks, ke, lam, det, ddet, row_scale)
1563 real, dimension(:), intent(in) :: a !< Lower diagonal of matrix (first entry unused) [T2 L-2 ~> s2 m-2]
1564 real, dimension(:), intent(in) :: c !< Upper diagonal of matrix (last entry unused) [T2 L-2 ~> s2 m-2]
1565 integer, intent(in) :: ks !< Starting index to use in determinant
1566 integer, intent(in) :: ke !< Ending index to use in determinant
1567 real, intent(in) :: lam !< Value subtracted from b [T2 L-2 ~> s2 m-2]
1568 real, intent(out):: det !< Determinant of the matrix in dynamically rescaled units that
1569 !! depend on the number of rows and the cumulative magnitude of
1570 !! det and are therefore difficult to interpret, but the units
1571 !! of det/ddet are always in [T2 L-2 ~> s2 m-2]
1572 real, intent(out):: ddet !< Derivative of determinant with lam in units that are dynamically
1573 !! rescaled along with those of det, such that the units of
1574 !! det/ddet are always in [T2 L-2 ~> s2 m-2]
1575 real, intent(in) :: row_scale !< A scaling factor of the rows of the matrix to
1576 !! limit the growth of the determinant [L2 s2 T-2 m-2 ~> 1]
1577 ! Local variables
1578 real :: detKm1, detKm2 ! Cumulative value of the determinant for the previous two layers in units
1579 ! that vary with the number of layers that have been worked on [various]
1580 real :: ddetKm1, ddetKm2 ! Derivative of the cumulative determinant with lam for the previous two
1581 ! layers [various], but the units of detKm1/ddetKm1 are [T2 L-2 ~> s2 m-2]
1582 real, parameter :: rescale = 1024.0**4 ! max value of determinant allowed before rescaling [nondim]
1583 real :: I_rescale ! inverse of rescale [nondim]
1584 integer :: k ! row (layer interface) index
1585
1586 i_rescale = 1.0 / rescale
1587
1588 detkm1 = 1.0 ; ddetkm1 = 0.0
1589 det = (a(ks)+c(ks)) - lam ; ddet = -1.0
1590 do k=ks+1,ke
1591 ! Shift variables and rescale rows to avoid over- or underflow.
1592 detkm2 = row_scale*detkm1 ; ddetkm2 = row_scale*ddetkm1
1593 detkm1 = row_scale*det ; ddetkm1 = row_scale*ddet
1594
1595 det = ((a(k)+c(k))-lam)*detkm1 - (a(k)*c(k-1))*detkm2
1596 ddet = ((a(k)+c(k))-lam)*ddetkm1 - (a(k)*c(k-1))*ddetkm2 - detkm1
1597
1598 ! Rescale det & ddet if det is getting too large or too small.
1599 if (abs(det) > rescale) then
1600 det = i_rescale*det ; detkm1 = i_rescale*detkm1
1601 ddet = i_rescale*ddet ; ddetkm1 = i_rescale*ddetkm1
1602 elseif (abs(det) < i_rescale) then
1603 det = rescale*det ; detkm1 = rescale*detkm1
1604 ddet = rescale*ddet ; ddetkm1 = rescale*ddetkm1
1605 endif
1606 enddo
1607
1608end subroutine tridiag_det
1609
1610!> Initialize control structure for MOM_wave_speed
1611subroutine wave_speed_init(CS, GV, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, &
1612 remap_answer_date, better_speed_est, om4_remap_via_sub_cells, &
1613 min_speed, wave_speed_tol, c1_thresh)
1614 type(wave_speed_cs), intent(inout) :: cs !< Wave speed control struct
1615 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1616 logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1617 !! barotropic mode instead of the first baroclinic mode.
1618 real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction of water column over
1619 !! which N2 is limited as monotonic for the purposes of
1620 !! calculating the vertical modal structure [nondim].
1621 real, optional, intent(in) :: mono_n2_depth !< The depth below which N2 is limited
1622 !! as monotonic for the purposes of calculating the
1623 !! vertical modal structure [H ~> m or kg m-2].
1624 logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
1625 !! that recover the remapping answers from 2018. Otherwise
1626 !! use more robust but mathematically equivalent expressions.
1627 integer, optional, intent(in) :: remap_answer_date !< The vintage of the order of arithmetic and expressions
1628 !! to use for remapping. Values below 20190101 recover the remapping
1629 !! answers from 2018, while higher values use more robust
1630 !! forms of the same remapping expressions.
1631 logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
1632 !! mode speed as the starting point for iterations.
1633 logical, optional, intent(in) :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells
1634 !! for calculating the EBT structure
1635 real, optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
1636 !! below which 0 is returned [L T-1 ~> m s-1].
1637 real, optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
1638 !! wave speeds [nondim]
1639 real, optional, intent(in) :: c1_thresh !< A minimal value of the first mode internal wave speed
1640 !! below which all higher mode speeds are not calculated but are
1641 !! simply reported as 0 [L T-1 ~> m s-1]. A non-negative value
1642 !! must be specified for wave_speeds to be used (but not wave_speed).
1643
1644 ! This include declares and sets the variable "version".
1645# include "version_variable.h"
1646 character(len=40) :: mdl = "MOM_wave_speed" ! This module's name.
1647
1648 cs%initialized = .true.
1649
1650 ! Write all relevant parameters to the model log.
1651 call log_version(mdl, version)
1652
1653 call wave_speed_set_param(cs, use_ebt_mode=use_ebt_mode, mono_n2_column_fraction=mono_n2_column_fraction, &
1654 mono_n2_depth=mono_n2_depth, better_speed_est=better_speed_est, &
1655 min_speed=min_speed, wave_speed_tol=wave_speed_tol, &
1656 remap_answers_2018=remap_answers_2018, remap_answer_date=remap_answer_date, &
1657 c1_thresh=c1_thresh)
1658
1659 ! The following remapping is only used for wave_speed with pre-2019 answers.
1660 if (cs%remap_answer_date < 20190101) &
1661 call initialize_remapping(cs%remap_2018_CS, 'PLM', boundary_extrapolation=.false., &
1662 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
1663 answer_date=cs%remap_answer_date, &
1664 h_neglect=1.0e-30*gv%m_to_H, h_neglect_edge=1.0e-10*gv%m_to_H)
1665
1666 ! This is used in wave_speeds in all cases, and in wave_speed with newer answers.
1667 call initialize_remapping(cs%remap_CS, 'PLM', boundary_extrapolation=.false., &
1668 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
1669 answer_date=cs%remap_answer_date, &
1670 h_neglect=gv%H_subroundoff, h_neglect_edge=gv%H_subroundoff)
1671
1672end subroutine wave_speed_init
1673
1674!> Sets internal parameters for MOM_wave_speed
1675subroutine wave_speed_set_param(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, &
1676 remap_answer_date, better_speed_est, min_speed, wave_speed_tol, c1_thresh)
1677 type(wave_speed_cs), intent(inout) :: cs
1678 !< Control structure for MOM_wave_speed
1679 logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1680 !! barotropic mode instead of the first baroclinic mode.
1681 real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction of water column over
1682 !! which N2 is limited as monotonic for the purposes of
1683 !! calculating the vertical modal structure [nondim].
1684 real, optional, intent(in) :: mono_n2_depth !< The depth below which N2 is limited
1685 !! as monotonic for the purposes of calculating the
1686 !! vertical modal structure [H ~> m or kg m-2].
1687 logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
1688 !! that recover the remapping answers from 2018. Otherwise
1689 !! use more robust but mathematically equivalent expressions.
1690 integer, optional, intent(in) :: remap_answer_date !< The vintage of the order of arithmetic and expressions
1691 !! to use for remapping. Values below 20190101 recover the remapping
1692 !! answers from 2018, while higher values use more robust
1693 !! forms of the same remapping expressions.
1694 logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
1695 !! mode speed as the starting point for iterations.
1696 real, optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
1697 !! below which 0 is returned [L T-1 ~> m s-1].
1698 real, optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
1699 !! wave speeds [nondim]
1700 real, optional, intent(in) :: c1_thresh !< A minimal value of the first mode internal wave speed
1701 !! below which all higher mode speeds are not calculated but are
1702 !! simply reported as 0 [L T-1 ~> m s-1]. A non-negative value
1703 !! must be specified for wave_speeds to be used (but not wave_speed).
1704
1705 if (present(use_ebt_mode)) cs%use_ebt_mode = use_ebt_mode
1706 if (present(mono_n2_column_fraction)) cs%mono_N2_column_fraction = mono_n2_column_fraction
1707 if (present(mono_n2_depth)) cs%mono_N2_depth = mono_n2_depth
1708 if (present(remap_answers_2018)) then
1709 if (remap_answers_2018) then
1710 cs%remap_answer_date = 20181231
1711 else
1712 cs%remap_answer_date = 20190101
1713 endif
1714 endif
1715 if (present(remap_answer_date)) cs%remap_answer_date = remap_answer_date
1716 if (present(better_speed_est)) cs%better_cg1_est = better_speed_est
1717 if (present(min_speed)) cs%min_speed2 = min_speed**2
1718 if (present(wave_speed_tol)) cs%wave_speed_tol = wave_speed_tol
1719 if (present(c1_thresh)) cs%c1_thresh = c1_thresh
1720
1721end subroutine wave_speed_set_param
1722
1723!> \namespace mom_wave_speed
1724
1725!!
1726!! Subroutine wave_speed() solves for the first baroclinic mode wave speed. (It could
1727!! solve for all the wave speeds, but the iterative approach taken here means
1728!! that this is not particularly efficient.)
1729!!
1730!! If `e(k)` is the perturbation interface height, this means solving for the
1731!! smallest eigenvalue (`lam` = 1/c^2) of the system
1732!!
1733!! \verbatim
1734!! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0
1735!! \endverbatim
1736!!
1737!! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving
1738!!
1739!! \verbatim
1740!! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0
1741!! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0
1742!! \endverbatim
1743!!
1744!! Here
1745!! \verbatim
1746!! Igl(k) = 1.0/(gprime(K)*h(k)) ; Igu(k) = 1.0/(gprime(K)*h(k-1))
1747!! \endverbatim
1748!!
1749!! Alternately, these same eigenvalues can be found from the second smallest
1750!! eigenvalue of the Montgomery potential (M(k)) calculation:
1751!!
1752!! \verbatim
1753!! -Igl(k)*M(k-1) + (Igl(k)+Igu(k+1)-lam)*M(k) - Igu(k+1)*M(k+1) = 0.0
1754!! \endverbatim
1755!!
1756!! with rigid lid and flat bottom boundary conditions
1757!!
1758!! \verbatim
1759!! (Igu(2)-lam)*M(1) - Igu(2)*M(2) = 0.0
1760!! -Igl(nz)*M(nz-1) + (Igl(nz)-lam)*M(nz) = 0.0
1761!! \endverbatim
1762!!
1763!! Note that the barotropic mode has been eliminated from the rigid lid
1764!! interface height equations, hence the matrix is one row smaller. Without
1765!! the rigid lid, the top boundary condition is simpler to implement with
1766!! the M equations.
1767
1768end module mom_wave_speed