MOM6
Data Types | Functions/Subroutines
mom_wave_structure Module Reference

Detailed Description

Vertical structure functions for first baroclinic mode wave speed.

Data Types

type  wave_structure_cs
 The control structure for the MOM_wave_structure module. More...
 

Functions/Subroutines

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. More...
 
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 stable variant that invokes the "Hallberg substitution" (TDMA_H). More...
 
subroutine, public wave_structure_init (Time, G, param_file, diag, CS)
 Allocate memory associated with the wave structure module and read parameters. More...
 

Function/Subroutine Documentation

◆ tridiag_solver()

subroutine mom_wave_structure::tridiag_solver ( real, dimension(:), intent(in)  a,
real, dimension(:), intent(in)  b,
real, dimension(:), intent(in)  c,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  y,
character(len=*), intent(in)  method,
real, dimension(:), intent(out)  x 
)
private

Solves a tri-diagonal system Ax=y using either the standard Thomas algorithm (TDMA_T) or its more stable variant that invokes the "Hallberg substitution" (TDMA_H).

Parameters
[in]alower diagonal with first entry equal to zero.
[in]bmiddle diagonal.
[in]cupper diagonal with last entry equal to zero.
[in]hvector of values that have already been added to b; used for systems of the form (e.g. average layer thickness in vertical diffusion case): [ -alpha(k-1/2) ] * e(k-1) + [ alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k) + [ -alpha(k+1/2) ] * e(k+1) = y(k) where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)], and c(k)=[-alpha(k+1/2)]. Only used with TDMA_H method.
[in]yvector of known values on right hand side.
[in]methodA string describing the algorithm to use
[out]xvector of unknown values to solve for.

Definition at line 562 of file MOM_wave_structure.F90.

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 

References mom_error_handler::mom_error().

Referenced by wave_structure().

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

◆ wave_structure()

subroutine, public mom_wave_structure::wave_structure ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  cn,
integer, intent(in)  ModeNum,
real, intent(in)  freq,
type(wave_structure_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional  En,
logical, intent(in), optional  full_halos 
)

This subroutine determines the internal wave velocity structure for any mode.

This subroutine solves for the eigen vector [vertical structure, e(k)] associated with the first baroclinic mode speed [i.e., smallest eigen value (lam = 1/c^2)] of the system d2e/dz2 = -(N2/cn2)e, or (A-lam*I)e = 0, where A = -(1/N2)(d2/dz2), lam = 1/c^2, and I is the identity matrix. 2nd order discretization in the vertical lets this system be represented as

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

with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving

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

where, upon noting N2 = reduced gravity/layer thickness, we get Igl(k) = 1.0/(gprime(k)*H(k)) ; Igu(k) = 1.0/(gprime(k)*H(k-1))

The eigen value for this system is approximated using "wave_speed." This subroutine uses these eigen values (mode speeds) to estimate the corresponding eigen vectors (velocity structure) using the "inverse iteration with shift" method. The algorithm is

Pick a starting vector reasonably close to mode structure and with unit magnitude, b_guess For n=1,2,3,... Solve (A-lam*I)e = e_guess for e Set e_guess=e/|e| and repeat, with each iteration refining the estimate of e

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvA structure pointing to various thermodynamic variables.
[in]cnThe (non-rotational) mode internal gravity wave speed [L T-1 ~> m s-1].
[in]modenumMode number
[in]freqIntrinsic wave frequency [T-1 ~> s-1].
csThe control structure returned by a previous call to wave_structure_init.
[in]enInternal wave energy density [R Z3 T-2 ~> J m-2]
[in]full_halosIf true, do the calculation over the entire computational domain.

Definition at line 92 of file MOM_wave_structure.F90.

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 

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

Referenced by mom_internal_tides::propagate_int_tide().

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

◆ wave_structure_init()

subroutine, public mom_wave_structure::wave_structure_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(in), target  diag,
type(wave_structure_cs), pointer  CS 
)

Allocate memory associated with the wave structure module and read parameters.

Parameters
[in]timeThe current model time.
[in]gThe ocean's grid structure.
[in]param_fileA structure to parse for run-time parameters.
[in]diagA structure that is used to regulate diagnostic output.
csA pointer that is set to point to the control structure for this module.

Definition at line 686 of file MOM_wave_structure.F90.

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 

References mom_error_handler::mom_error().

Referenced by mom_internal_tides::internal_tides_init().

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