MOM6
MOM_wave_structure.F90
Go to the documentation of this file.
1 !> Vertical structure functions for first baroclinic mode wave speed
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 ! By Benjamin Mater & Robert Hallberg, 2015
7 
8 ! The subroutine in this module calculates the vertical structure
9 ! functions of the first baroclinic mode internal wave speed.
10 ! Calculation of interface values is the same as done in
11 ! MOM_wave_speed by Hallberg, 2008.
12 
13 use mom_debugging, only : isnan => is_nan
15 use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
17 use mom_error_handler, only : mom_error, fatal, warning
19 use mom_grid, only : ocean_grid_type
23 
24 implicit none ; private
25 
26 #include <MOM_memory.h>
27 
29 
30 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34 
35 !> The control structure for the MOM_wave_structure module
36 type, public :: wave_structure_cs ; !private
37  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
38  !! regulate the timing of diagnostic output.
39  real, allocatable, dimension(:,:,:) :: w_strct
40  !< Vertical structure of vertical velocity (normalized) [m s-1].
41  real, allocatable, dimension(:,:,:) :: u_strct
42  !< Vertical structure of horizontal velocity (normalized) [m s-1].
43  real, allocatable, dimension(:,:,:) :: w_profile
44  !< Vertical profile of w_hat(z), where
45  !! w(x,y,z,t) = w_hat(z)*exp(i(kx+ly-freq*t)) is the full time-
46  !! varying vertical velocity with w_hat(z) = W0*w_strct(z) [Z T-1 ~> m s-1].
47  real, allocatable, dimension(:,:,:) :: uavg_profile
48  !< Vertical profile of the magnitude of horizontal velocity,
49  !! (u^2+v^2)^0.5, averaged over a period [L T-1 ~> m s-1].
50  real, allocatable, dimension(:,:,:) :: z_depths
51  !< Depths of layer interfaces [m].
52  real, allocatable, dimension(:,:,:) :: n2
53  !< Squared buoyancy frequency at each interface [s-2].
54  integer, allocatable, dimension(:,:):: num_intfaces
55  !< Number of layer interfaces (including surface and bottom)
56  real :: int_tide_source_x !< X Location of generation site
57  !! for internal tide for testing (BDM)
58  real :: int_tide_source_y !< Y Location of generation site
59  !! for internal tide for testing (BDM)
60 
61 end type wave_structure_cs
62 
63 contains
64 
65 !> This subroutine determines the internal wave velocity structure for any mode.
66 !!
67 !! This subroutine solves for the eigen vector [vertical structure, e(k)] associated with
68 !! the first baroclinic mode speed [i.e., smallest eigen value (lam = 1/c^2)] of the
69 !! system d2e/dz2 = -(N2/cn2)e, or (A-lam*I)e = 0, where A = -(1/N2)(d2/dz2), lam = 1/c^2,
70 !! and I is the identity matrix. 2nd order discretization in the vertical lets this system
71 !! be represented as
72 !!
73 !! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0
74 !!
75 !! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving
76 !!
77 !! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0
78 !! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0
79 !!
80 !! where, upon noting N2 = reduced gravity/layer thickness, we get
81 !! Igl(k) = 1.0/(gprime(k)*H(k)) ; Igu(k) = 1.0/(gprime(k)*H(k-1))
82 !!
83 !! The eigen value for this system is approximated using "wave_speed." This subroutine uses
84 !! these eigen values (mode speeds) to estimate the corresponding eigen vectors (velocity
85 !! structure) using the "inverse iteration with shift" method. The algorithm is
86 !!
87 !! Pick a starting vector reasonably close to mode structure and with unit magnitude, b_guess
88 !! For n=1,2,3,...
89 !! Solve (A-lam*I)e = e_guess for e
90 !! Set e_guess=e/|e| and repeat, with each iteration refining the estimate of e
91 subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halos)
92  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
93  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
94  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
95  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
96  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
97  !! thermodynamic variables.
98  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: cn !< The (non-rotational) mode internal
99  !! gravity wave speed [L T-1 ~> m s-1].
100  integer, intent(in) :: modenum !< Mode number
101  real, intent(in) :: freq !< Intrinsic wave frequency [T-1 ~> s-1].
102  type(wave_structure_cs), pointer :: cs !< The control structure returned by a
103  !! previous call to wave_structure_init.
104  real, dimension(SZI_(G),SZJ_(G)), &
105  optional, intent(in) :: en !< Internal wave energy density [R Z3 T-2 ~> J m-2]
106  logical, optional, intent(in) :: full_halos !< If true, do the calculation
107  !! over the entire computational domain.
108  ! Local variables
109  real, dimension(SZK_(G)+1) :: &
110  drho_dt, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
111  drho_ds, & ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
112  pres, & ! Interface pressure [Pa]
113  t_int, & ! Temperature interpolated to interfaces [degC]
114  s_int, & ! Salinity interpolated to interfaces [ppt]
115  gprime ! The reduced gravity across each interface [m2 Z-1 s-2 ~> m s-2].
116  real, dimension(SZK_(G)) :: &
117  igl, igu ! The inverse of the reduced gravity across an interface times
118  ! the thickness of the layer below (Igl) or above (Igu) it [s2 m-2].
119  real, dimension(SZK_(G),SZI_(G)) :: &
120  hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m]
121  tf, & ! Layer temperatures after very thin layers are combined [degC]
122  sf, & ! Layer salinities after very thin layers are combined [ppt]
123  rf ! Layer densities after very thin layers are combined [R ~> kg m-3]
124  real, dimension(SZK_(G)) :: &
125  hc, & ! A column of layer thicknesses after convective istabilities are removed [Z ~> m]
126  tc, & ! A column of layer temperatures after convective istabilities are removed [degC]
127  sc, & ! A column of layer salinites after convective istabilities are removed [ppt]
128  rc, & ! A column of layer densities after convective istabilities are removed [R ~> kg m-3]
129  det, ddet
130  real, dimension(SZI_(G),SZJ_(G)) :: &
131  htot ! The vertical sum of the thicknesses [Z ~> m]
132  real :: lam
133  real :: min_h_frac
134  real :: h_to_pres
135  real, dimension(SZI_(G)) :: &
136  hmin, & ! Thicknesses [Z ~> m]
137  h_here, & ! A thickness [Z ~> m]
138  hxt_here, & ! A layer integrated temperature [degC Z ~> degC m]
139  hxs_here, & ! A layer integrated salinity [ppt Z ~> ppt m]
140  hxr_here ! A layer integrated density [R Z ~> kg m-2]
141  real :: speed2_tot
142  real :: i_hnew ! The inverse of a new layer thickness [Z-1 ~> m-1]
143  real :: drxh_sum ! The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-2]
144  real, parameter :: tol1 = 0.0001, tol2 = 0.001
145  real, pointer, dimension(:,:,:) :: t => null(), s => null()
146  real :: g_rho0 ! G_Earth/Rho0 in [m2 s-2 Z-1 R-1 ~> m4 s-2 kg-1].
147  ! real :: rescale, I_rescale
148  integer :: kf(szi_(g))
149  integer, parameter :: max_itt = 1 ! number of times to iterate in solving for eigenvector
150  real :: cg_subro ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1]
151  real, parameter :: a_int = 0.5 ! value of normalized integral: \int(w_strct^2)dz = a_int
152  real :: i_a_int ! inverse of a_int
153  real :: f2 ! squared Coriolis frequency [T-2 ~> s-2]
154  real :: kmag2 ! magnitude of horizontal wave number squared
155  logical :: use_eos ! If true, density is calculated from T & S using an
156  ! equation of state.
157  real, dimension(SZK_(G)+1) :: w_strct, u_strct, w_profile, uavg_profile, z_int, n2
158  ! local representations of variables in CS; note,
159  ! not all rows will be filled if layers get merged!
160  real, dimension(SZK_(G)+1) :: w_strct2, u_strct2
161  ! squared values
162  real, dimension(SZK_(G)) :: dz ! thicknesses of merged layers (same as Hc I hope)
163  ! real, dimension(SZK_(G)+1) :: dWdz_profile ! profile of dW/dz
164  real :: w2avg ! average of squared vertical velocity structure funtion
165  real :: int_dwdz2
166  real :: int_w2
167  real :: int_n2w2
168  real :: ke_term ! terms in vertically averaged energy equation
169  real :: pe_term ! terms in vertically averaged energy equation
170  real :: w0 ! A vertical velocity magnitude [Z T-1 ~> m s-1]
171  real :: gp_unscaled ! A version of gprime rescaled to [m s-2].
172  real, dimension(SZK_(G)-1) :: lam_z ! product of eigen value and gprime(k); one value for each
173  ! interface (excluding surface and bottom)
174  real, dimension(SZK_(G)-1) :: a_diag, b_diag, c_diag
175  ! diagonals of tridiagonal matrix; one value for each
176  ! interface (excluding surface and bottom)
177  real, dimension(SZK_(G)-1) :: e_guess ! guess at eigen vector with unit amplitde (for TDMA)
178  real, dimension(SZK_(G)-1) :: e_itt ! improved guess at eigen vector (from TDMA)
179  real :: pi
180  integer :: kc
181  integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop
182 
183  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
184  i_a_int = 1/a_int
185 
186  !if (present(CS)) then
187  if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_structure: "// &
188  "Module must be initialized before it is used.")
189  !endif
190 
191  if (present(full_halos)) then ; if (full_halos) then
192  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
193  endif ; endif
194 
195  pi = (4.0*atan(1.0))
196 
197  s => tv%S ; t => tv%T
198  g_rho0 = us%L_T_to_m_s**2 * gv%g_Earth / gv%Rho0
199  cg_subro = 1e-100*us%m_s_to_L_T ! The hard-coded value here might need to increase.
200  use_eos = associated(tv%eqn_of_state)
201 
202  h_to_pres = gv%Z_to_H*gv%H_to_Pa
203  ! rescale = 1024.0**4 ; I_rescale = 1.0/rescale
204 
205  min_h_frac = tol1 / real(nz)
206 
207  do j=js,je
208  ! First merge very thin layers with the one above (or below if they are
209  ! at the top). This also transposes the row order so that columns can
210  ! be worked upon one at a time.
211  do i=is,ie ; htot(i,j) = 0.0 ; enddo
212  do k=1,nz ; do i=is,ie ; htot(i,j) = htot(i,j) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo
213 
214  do i=is,ie
215  hmin(i) = htot(i,j)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
216  hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
217  enddo
218  if (use_eos) then
219  do k=1,nz ; do i=is,ie
220  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
221  hf(kf(i),i) = h_here(i)
222  tf(kf(i),i) = hxt_here(i) / h_here(i)
223  sf(kf(i),i) = hxs_here(i) / h_here(i)
224  kf(i) = kf(i) + 1
225 
226  ! Start a new layer
227  h_here(i) = h(i,j,k)*gv%H_to_Z
228  hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
229  hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
230  else
231  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
232  hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
233  hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
234  endif
235  enddo ; enddo
236  do i=is,ie ; if (h_here(i) > 0.0) then
237  hf(kf(i),i) = h_here(i)
238  tf(kf(i),i) = hxt_here(i) / h_here(i)
239  sf(kf(i),i) = hxs_here(i) / h_here(i)
240  endif ; enddo
241  else
242  do k=1,nz ; do i=is,ie
243  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
244  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
245  kf(i) = kf(i) + 1
246 
247  ! Start a new layer
248  h_here(i) = h(i,j,k)*gv%H_to_Z
249  hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
250  else
251  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
252  hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
253  endif
254  enddo ; enddo
255  do i=is,ie ; if (h_here(i) > 0.0) then
256  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
257  endif ; enddo
258  endif ! use_EOS?
259 
260  ! From this point, we can work on individual columns without causing memory
261  ! to have page faults.
262  do i=is,ie ; if (cn(i,j)>0.0)then
263  !----for debugging, remove later----
264  ig = i + g%idg_offset ; jg = j + g%jdg_offset
265  !if (ig == CS%int_tide_source_x .and. jg == CS%int_tide_source_y) then
266  !-----------------------------------
267  if (g%mask2dT(i,j) > 0.5) then
268 
269  lam = 1/(us%L_T_to_m_s**2 * cn(i,j)**2)
270 
271  ! Calculate drxh_sum
272  if (use_eos) then
273  pres(1) = 0.0
274  do k=2,kf(i)
275  pres(k) = pres(k-1) + h_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  enddo
279  call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, 2, &
280  kf(i)-1, tv%eqn_of_state, scale=us%kg_m3_to_R)
281 
282  ! Sum the reduced gravities to find out how small a density difference
283  ! is negligibly small.
284  drxh_sum = 0.0
285  do k=2,kf(i)
286  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
287  max(0.0,drho_dt(k)*(tf(k,i)-tf(k-1,i)) + &
288  drho_ds(k)*(sf(k,i)-sf(k-1,i)))
289  enddo
290  else
291  drxh_sum = 0.0
292  do k=2,kf(i)
293  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
294  max(0.0,rf(k,i)-rf(k-1,i))
295  enddo
296  endif ! use_EOS?
297 
298  ! Find gprime across each internal interface, taking care of convective
299  ! instabilities by merging layers.
300  if (drxh_sum >= 0.0) then
301  ! Merge layers to eliminate convective instabilities or exceedingly
302  ! small reduced gravities.
303  if (use_eos) then
304  kc = 1
305  hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
306  do k=2,kf(i)
307  if ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
308  (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum) then
309  ! Merge this layer with the one above and backtrack.
310  i_hnew = 1.0 / (hc(kc) + hf(k,i))
311  tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
312  sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
313  hc(kc) = (hc(kc) + hf(k,i))
314  ! Backtrack to remove any convective instabilities above... Note
315  ! that the tolerance is a factor of two larger, to avoid limit how
316  ! far back we go.
317  do k2=kc,2,-1
318  if ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
319  (hc(k2) + hc(k2-1)) < tol2*drxh_sum) then
320  ! Merge the two bottommost layers. At this point kc = k2.
321  i_hnew = 1.0 / (hc(kc) + hc(kc-1))
322  tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
323  sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
324  hc(kc-1) = (hc(kc) + hc(kc-1))
325  kc = kc - 1
326  else ; exit ; endif
327  enddo
328  else
329  ! Add a new layer to the column.
330  kc = kc + 1
331  drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
332  tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
333  endif
334  enddo
335  ! At this point there are kc layers and the gprimes should be positive.
336  do k=2,kc ! Revisit this if non-Boussinesq.
337  gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + &
338  drho_ds(k)*(sc(k)-sc(k-1)))
339  enddo
340  else ! .not.use_EOS
341  ! Do the same with density directly...
342  kc = 1
343  hc(1) = hf(1,i) ; rc(1) = rf(1,i)
344  do k=2,kf(i)
345  if ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum) then
346  ! Merge this layer with the one above and backtrack.
347  rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
348  hc(kc) = (hc(kc) + hf(k,i))
349  ! Backtrack to remove any convective instabilities above... Note
350  ! that the tolerance is a factor of two larger, to avoid limit how
351  ! far back we go.
352  do k2=kc,2,-1
353  if ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum) then
354  ! Merge the two bottommost layers. At this point kc = k2.
355  rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
356  hc(kc-1) = (hc(kc) + hc(kc-1))
357  kc = kc - 1
358  else ; exit ; endif
359  enddo
360  else
361  ! Add a new layer to the column.
362  kc = kc + 1
363  rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
364  endif
365  enddo
366  ! At this point there are kc layers and the gprimes should be positive.
367  do k=2,kc ! Revisit this if non-Boussinesq.
368  gprime(k) = g_rho0 * (rc(k)-rc(k-1))
369  enddo
370  endif ! use_EOS?
371 
372  !-----------------NOW FIND WAVE STRUCTURE-------------------------------------
373  ! Construct and solve tridiagonal system for the interior interfaces
374  ! Note that kc = number of layers,
375  ! kc+1 = nzm = number of interfaces,
376  ! kc-1 = number of interior interfaces (excluding surface and bottom)
377  ! Also, note that "K" refers to an interface, while "k" refers to the layer below.
378  ! Need at least 3 layers (2 internal interfaces) to generate a matrix, also
379  ! need number of layers to be greater than the mode number
380  if (kc >= modenum + 1) then
381  ! Set depth at surface
382  z_int(1) = 0.0
383  ! Calculate Igu, Igl, depth, and N2 at each interior interface
384  ! [excludes surface (K=1) and bottom (K=kc+1)]
385  do k=2,kc
386  igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
387  z_int(k) = z_int(k-1) + hc(k-1)
388  n2(k) = us%m_to_Z**2*gprime(k)/(0.5*(hc(k)+hc(k-1)))
389  enddo
390  ! Set stratification for surface and bottom (setting equal to nearest interface for now)
391  n2(1) = n2(2) ; n2(kc+1) = n2(kc)
392  ! Calcualte depth at bottom
393  z_int(kc+1) = z_int(kc)+hc(kc)
394  ! check that thicknesses sum to total depth
395  if (abs(z_int(kc+1)-htot(i,j)) > 1.e-14*htot(i,j)) then
396  call mom_error(fatal, "wave_structure: mismatch in total depths")
397  endif
398 
399  ! Note that many of the calcluation from here on revert to using vertical
400  ! distances in m, not Z.
401 
402  ! Populate interior rows of tridiagonal matrix; must multiply through by
403  ! gprime to get tridiagonal matrix to the symmetrical form:
404  ! [-1/H(k-1)]e(k-1) + [1/H(k-1)+1/H(k)-lam_z]e(k) + [-1/H(k)]e(k+1) = 0,
405  ! where lam_z = lam*gprime is now a function of depth.
406  ! Frist, populate interior rows
407  do k=3,kc-1
408  row = k-1 ! indexing for TD matrix rows
409  gp_unscaled = us%m_to_Z*gprime(k)
410  lam_z(row) = lam*gp_unscaled
411  a_diag(row) = gp_unscaled*(-igu(k))
412  b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
413  c_diag(row) = gp_unscaled*(-igl(k))
414  if (isnan(lam_z(row)))then ; print *, "Wave_structure: lam_z(row) is NAN" ; endif
415  if (isnan(a_diag(row)))then ; print *, "Wave_structure: a(k) is NAN" ; endif
416  if (isnan(b_diag(row)))then ; print *, "Wave_structure: b(k) is NAN" ; endif
417  if (isnan(c_diag(row)))then ; print *, "Wave_structure: c(k) is NAN" ; endif
418  enddo
419  ! Populate top row of tridiagonal matrix
420  k=2 ; row = k-1 ;
421  gp_unscaled = us%m_to_Z*gprime(k)
422  lam_z(row) = lam*gp_unscaled
423  a_diag(row) = 0.0
424  b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
425  c_diag(row) = gp_unscaled*(-igl(k))
426  ! Populate bottom row of tridiagonal matrix
427  k=kc ; row = k-1
428  gp_unscaled = us%m_to_Z*gprime(k)
429  lam_z(row) = lam*gp_unscaled
430  a_diag(row) = gp_unscaled*(-igu(k))
431  b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
432  c_diag(row) = 0.0
433 
434  ! Guess a vector shape to start with (excludes surface and bottom)
435  e_guess(1:kc-1) = sin((z_int(2:kc)/htot(i,j)) *pi)
436  e_guess(1:kc-1) = e_guess(1:kc-1)/sqrt(sum(e_guess(1:kc-1)**2))
437 
438  ! Perform inverse iteration with tri-diag solver
439  do itt=1,max_itt
440  call tridiag_solver(a_diag(1:kc-1),b_diag(1:kc-1),c_diag(1:kc-1), &
441  -lam_z(1:kc-1),e_guess(1:kc-1),"TDMA_H",e_itt)
442  e_guess(1:kc-1) = e_itt(1:kc-1) / sqrt(sum(e_itt(1:kc-1)**2))
443  enddo ! itt-loop
444  w_strct(2:kc) = e_guess(1:kc-1)
445  w_strct(1) = 0.0 ! rigid lid at surface
446  w_strct(kc+1) = 0.0 ! zero-flux at bottom
447 
448  ! Check to see if solver worked
449  ig_stop = 0 ; jg_stop = 0
450  if (isnan(sum(w_strct(1:kc+1))))then
451  print *, "Wave_structure: w_strct has a NAN at ig=", ig, ", jg=", jg
452  if (i<g%isc .or. i>g%iec .or. j<g%jsc .or. j>g%jec)then
453  print *, "This is occuring at a halo point."
454  endif
455  ig_stop = ig ; jg_stop = jg
456  endif
457 
458  ! Normalize vertical structure function of w such that
459  ! \int(w_strct)^2dz = a_int (a_int could be any value, e.g., 0.5)
460  nzm = kc+1 ! number of layer interfaces after merging
461  !(including surface and bottom)
462  w2avg = 0.0
463  do k=1,nzm-1
464  dz(k) = us%Z_to_m*hc(k)
465  w2avg = w2avg + 0.5*(w_strct(k)**2+w_strct(k+1)**2)*dz(k)
466  enddo
467  !### Some mathematical cancellations could occur in the next two lines.
468  w2avg = w2avg / htot(i,j)
469  w_strct(:) = w_strct(:) / sqrt(htot(i,j)*w2avg*i_a_int)
470 
471  ! Calculate vertical structure function of u (i.e. dw/dz)
472  do k=2,nzm-1
473  u_strct(k) = 0.5*((w_strct(k-1) - w_strct(k) )/dz(k-1) + &
474  (w_strct(k) - w_strct(k+1))/dz(k))
475  enddo
476  u_strct(1) = (w_strct(1) - w_strct(2) )/dz(1)
477  u_strct(nzm) = (w_strct(nzm-1)- w_strct(nzm))/dz(nzm-1)
478 
479  ! Calculate wavenumber magnitude
480  f2 = g%CoriolisBu(i,j)**2
481  !f2 = 0.25*US%s_to_T**2 *((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
482  ! (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
483  kmag2 = us%m_to_L**2 * (freq**2 - f2) / (cn(i,j)**2 + cg_subro**2)
484 
485  ! Calculate terms in vertically integrated energy equation
486  int_dwdz2 = 0.0 ; int_w2 = 0.0 ; int_n2w2 = 0.0
487  u_strct2(:) = u_strct(1:nzm)**2
488  w_strct2(:) = w_strct(1:nzm)**2
489  ! vertical integration with Trapezoidal rule
490  do k=1,nzm-1
491  int_dwdz2 = int_dwdz2 + 0.5*(u_strct2(k)+u_strct2(k+1)) * us%m_to_Z*dz(k)
492  int_w2 = int_w2 + 0.5*(w_strct2(k)+w_strct2(k+1)) * us%m_to_Z*dz(k)
493  int_n2w2 = int_n2w2 + 0.5*(w_strct2(k)*n2(k)+w_strct2(k+1)*n2(k+1)) * us%m_to_Z*dz(k)
494  enddo
495 
496  ! Back-calculate amplitude from energy equation
497  if (present(en) .and. (freq**2*kmag2 > 0.0)) then
498  ! Units here are [R
499  ke_term = 0.25*gv%Rho0*( ((freq**2 + f2) / (freq**2*kmag2))*int_dwdz2 + int_w2 )
500  pe_term = 0.25*gv%Rho0*( int_n2w2 / (us%s_to_T*freq)**2 )
501  if (en(i,j) >= 0.0) then
502  w0 = sqrt( en(i,j) / (ke_term + pe_term) )
503  else
504  call mom_error(warning, "wave_structure: En < 0.0; setting to W0 to 0.0")
505  print *, "En(i,j)=", en(i,j), " at ig=", ig, ", jg=", jg
506  w0 = 0.0
507  endif
508  ! Calculate actual vertical velocity profile and derivative
509  w_profile(:) = w0*w_strct(:)
510  ! dWdz_profile(:) = W0*u_strct(:)
511  ! Calculate average magnitude of actual horizontal velocity over a period
512  uavg_profile(:) = us%Z_to_L*abs(w0*u_strct(:)) * sqrt((freq**2 + f2) / (2.0*freq**2*kmag2))
513  else
514  w_profile(:) = 0.0
515  ! dWdz_profile(:) = 0.0
516  uavg_profile(:) = 0.0
517  endif
518 
519  ! Store values in control structure
520  cs%w_strct(i,j,1:nzm) = w_strct(:)
521  cs%u_strct(i,j,1:nzm) = u_strct(:)
522  cs%W_profile(i,j,1:nzm) = w_profile(:)
523  cs%Uavg_profile(i,j,1:nzm)= uavg_profile(:)
524  cs%z_depths(i,j,1:nzm) = us%Z_to_m*z_int(:)
525  cs%N2(i,j,1:nzm) = n2(:)
526  cs%num_intfaces(i,j) = nzm
527  else
528  ! If not enough layers, default to zero
529  nzm = kc+1
530  cs%w_strct(i,j,1:nzm) = 0.0
531  cs%u_strct(i,j,1:nzm) = 0.0
532  cs%W_profile(i,j,1:nzm) = 0.0
533  cs%Uavg_profile(i,j,1:nzm)= 0.0
534  cs%z_depths(i,j,1:nzm) = 0.0 ! could use actual values
535  cs%N2(i,j,1:nzm) = 0.0 ! could use with actual values
536  cs%num_intfaces(i,j) = nzm
537  endif ! kc >= 3 and kc > ModeNum + 1?
538  endif ! drxh_sum >= 0?
539  !else ! if at test point - delete later
540  ! return ! if at test point - delete later
541  !endif ! if at test point - delete later
542  endif ! mask2dT > 0.5?
543  else
544  ! if cn=0.0, default to zero
545  nzm = nz+1! could use actual values
546  cs%w_strct(i,j,1:nzm) = 0.0
547  cs%u_strct(i,j,1:nzm) = 0.0
548  cs%W_profile(i,j,1:nzm) = 0.0
549  cs%Uavg_profile(i,j,1:nzm)= 0.0
550  cs%z_depths(i,j,1:nzm) = 0.0 ! could use actual values
551  cs%N2(i,j,1:nzm) = 0.0 ! could use with actual values
552  cs%num_intfaces(i,j) = nzm
553  endif ; enddo ! if cn>0.0? ; i-loop
554  enddo ! j-loop
555 
556 end subroutine wave_structure
557 
558 !> Solves a tri-diagonal system Ax=y using either the standard
559 !! Thomas algorithm (TDMA_T) or its more stable variant that invokes the
560 !! "Hallberg substitution" (TDMA_H).
561 subroutine tridiag_solver(a, b, c, h, y, method, x)
562  real, dimension(:), intent(in) :: a !< lower diagonal with first entry equal to zero.
563  real, dimension(:), intent(in) :: b !< middle diagonal.
564  real, dimension(:), intent(in) :: c !< upper diagonal with last entry equal to zero.
565  real, dimension(:), intent(in) :: h !< vector of values that have already been added to b; used
566  !! for systems of the form (e.g. average layer thickness in vertical diffusion case):
567  !! [ -alpha(k-1/2) ] * e(k-1) +
568  !! [ alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k) +
569  !! [ -alpha(k+1/2) ] * e(k+1) = y(k)
570  !! where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)],
571  !! and c(k)=[-alpha(k+1/2)]. Only used with TDMA_H method.
572  real, dimension(:), intent(in) :: y !< vector of known values on right hand side.
573  character(len=*), intent(in) :: method !< A string describing the algorithm to use
574  real, dimension(:), intent(out) :: x !< vector of unknown values to solve for.
575  ! Local variables
576  integer :: nrow ! number of rows in A matrix
577 ! real, allocatable, dimension(:,:) :: A_check ! for solution checking
578 ! real, allocatable, dimension(:) :: y_check ! for solution checking
579  real, allocatable, dimension(:) :: c_prime, y_prime, q, alpha
580  ! intermediate values for solvers
581  real :: Q_prime, beta ! intermediate values for solver
582  integer :: k ! row (e.g. interface) index
583  integer :: i,j
584 
585  nrow = size(y)
586  allocate(c_prime(nrow))
587  allocate(y_prime(nrow))
588  allocate(q(nrow))
589  allocate(alpha(nrow))
590 ! allocate(A_check(nrow,nrow))
591 ! allocate(y_check(nrow))
592 
593  if (method == 'TDMA_T') then
594  ! Standard Thomas algoritim (4th variant).
595  ! Note: Requires A to be non-singular for accuracy/stability
596  c_prime(:) = 0.0 ; y_prime(:) = 0.0
597  c_prime(1) = c(1)/b(1) ; y_prime(1) = y(1)/b(1)
598 
599  ! Forward sweep
600  do k=2,nrow-1
601  c_prime(k) = c(k)/(b(k)-a(k)*c_prime(k-1))
602  enddo
603  !print *, 'c_prime=', c_prime(1:nrow)
604  do k=2,nrow
605  y_prime(k) = (y(k)-a(k)*y_prime(k-1))/(b(k)-a(k)*c_prime(k-1))
606  enddo
607  !print *, 'y_prime=', y_prime(1:nrow)
608  x(nrow) = y_prime(nrow)
609 
610  ! Backward sweep
611  do k=nrow-1,1,-1
612  x(k) = y_prime(k)-c_prime(k)*x(k+1)
613  enddo
614  !print *, 'x=',x(1:nrow)
615 
616  ! Check results - delete later
617  !do j=1,nrow ; do i=1,nrow
618  ! if (i==j)then ; A_check(i,j) = b(i)
619  ! elseif (i==j+1)then ; A_check(i,j) = a(i)
620  ! elseif (i==j-1)then ; A_check(i,j) = c(i)
621  ! endif
622  !enddo ; enddo
623  !print *, 'A(2,1),A(2,2),A(1,2)=', A_check(2,1), A_check(2,2), A_check(1,2)
624  !y_check = matmul(A_check,x)
625  !if (all(y_check /= y))then
626  ! print *, "tridiag_solver: Uh oh, something's not right!"
627  ! print *, "y=", y
628  ! print *, "y_check=", y_check
629  !endif
630 
631  elseif (method == 'TDMA_H') then
632  ! Thomas algoritim (4th variant) w/ Hallberg substitution.
633  ! For a layered system where k is at interfaces, alpha{k+1/2} refers to
634  ! some property (e.g. inverse thickness for mode-structure problem) of the
635  ! layer below and alpha{k-1/2} refers to the layer above.
636  ! Here, alpha(k)=alpha{k+1/2} and alpha(k-1)=alpha{k-1/2}.
637  ! Strictly speaking, this formulation requires A to be a non-singular,
638  ! symmetric, diagonally dominant matrix, with h>0.
639  ! Need to add a check for these conditions.
640  do k=1,nrow-1
641  if (abs(a(k+1)-c(k)) > 1.e-10*(abs(a(k+1))+abs(c(k)))) then
642  call mom_error(warning, "tridiag_solver: matrix not symmetric; need symmetry when invoking TDMA_H")
643  endif
644  enddo
645  alpha = -c
646  ! Alpha of the bottom-most layer is not necessarily zero. Therefore,
647  ! back out the value from the provided b(nrow and h(nrow) values
648  alpha(nrow) = b(nrow)-h(nrow)-alpha(nrow-1)
649  ! Prime other variables
650  beta = 1/b(1)
651  y_prime(:) = 0.0 ; q(:) = 0.0
652  y_prime(1) = beta*y(1) ; q(1) = beta*alpha(1)
653  q_prime = 1-q(1)
654 
655  ! Forward sweep
656  do k=2,nrow-1
657  beta = 1/(h(k)+alpha(k-1)*q_prime+alpha(k))
658  if (isnan(beta))then ; print *, "Tridiag_solver: beta is NAN" ; endif
659  q(k) = beta*alpha(k)
660  y_prime(k) = beta*(y(k)+alpha(k-1)*y_prime(k-1))
661  q_prime = beta*(h(k)+alpha(k-1)*q_prime)
662  enddo
663  if ((h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow)) == 0.0)then
664  call mom_error(fatal, "Tridiag_solver: this system is not stable.") ! ; overriding beta(nrow)
665  ! This has hard-coded dimensions: beta = 1/(1e-15) ! place holder for unstable systems - delete later
666  else
667  beta = 1/(h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow))
668  endif
669  y_prime(nrow) = beta*(y(nrow)+alpha(nrow-1)*y_prime(nrow-1))
670  x(nrow) = y_prime(nrow)
671  ! Backward sweep
672  do k=nrow-1,1,-1
673  x(k) = y_prime(k)+q(k)*x(k+1)
674  enddo
675  !print *, 'yprime=',y_prime(1:nrow)
676  !print *, 'x=',x(1:nrow)
677  endif
678 
679  deallocate(c_prime,y_prime,q,alpha)
680 ! deallocate(A_check,y_check)
681 
682 end subroutine tridiag_solver
683 
684 !> Allocate memory associated with the wave structure module and read parameters.
685 subroutine wave_structure_init(Time, G, param_file, diag, CS)
686  type(time_type), intent(in) :: time !< The current model time.
687  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
688  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
689  !! parameters.
690  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
691  !! diagnostic output.
692  type(wave_structure_cs), pointer :: cs !< A pointer that is set to point to the
693  !! control structure for this module.
694 ! This include declares and sets the variable "version".
695 #include "version_variable.h"
696  character(len=40) :: mdl = "MOM_wave_structure" ! This module's name.
697  integer :: isd, ied, jsd, jed, nz
698 
699  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
700 
701  if (associated(cs)) then
702  call mom_error(warning, "wave_structure_init called with an "// &
703  "associated control structure.")
704  return
705  else ; allocate(cs) ; endif
706 
707  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
708  "X Location of generation site for internal tide", default=1.)
709  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
710  "Y Location of generation site for internal tide", default=1.)
711 
712  cs%diag => diag
713 
714  ! Allocate memory for variable in control structure; note,
715  ! not all rows will be filled if layers get merged!
716  allocate(cs%w_strct(isd:ied,jsd:jed,nz+1))
717  allocate(cs%u_strct(isd:ied,jsd:jed,nz+1))
718  allocate(cs%W_profile(isd:ied,jsd:jed,nz+1))
719  allocate(cs%Uavg_profile(isd:ied,jsd:jed,nz+1))
720  allocate(cs%z_depths(isd:ied,jsd:jed,nz+1))
721  allocate(cs%N2(isd:ied,jsd:jed,nz+1))
722  allocate(cs%num_intfaces(isd:ied,jsd:jed))
723 
724  ! Write all relevant parameters to the model log.
725  call log_version(param_file, mdl, version, "")
726 
727 end subroutine wave_structure_init
728 
729 end module mom_wave_structure
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Call this subroutine to determine whether the averaging is currently enabled. .true....
Definition: MOM_diag_mediator.F90:1850
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:78
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_wave_structure
Vertical structure functions for first baroclinic mode wave speed.
Definition: MOM_wave_structure.F90:2
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_wave_structure::wave_structure
subroutine, public wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halos)
This subroutine determines the internal wave velocity structure for any mode.
Definition: MOM_wave_structure.F90:92
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_structure::wave_structure_cs
The control structure for the MOM_wave_structure module.
Definition: MOM_wave_structure.F90:36
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:70
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_structure::wave_structure_init
subroutine, public wave_structure_init(Time, G, param_file, diag, CS)
Allocate memory associated with the wave structure module and read parameters.
Definition: MOM_wave_structure.F90:686
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
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:239
mom_wave_structure::tridiag_solver
subroutine tridiag_solver(a, b, c, h, y, method, x)
Solves a tri-diagonal system Ax=y using either the standard Thomas algorithm (TDMA_T) or its more sta...
Definition: MOM_wave_structure.F90:562