MOM6
Data Types | Functions/Subroutines
mom_wave_structure Module Reference

Data Types

type  wave_structure_cs
 

Functions/Subroutines

subroutine, public wave_structure (h, tv, G, GV, 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)
 This subroutine solves a tri-diagonal system Ax=y using either the standard Thomas algorithim (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)
 

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

This subroutine solves a tri-diagonal system Ax=y using either the standard Thomas algorithim (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.
[out]xvector of unknown values to solve for.

Definition at line 619 of file MOM_wave_structure.F90.

References mom_error_handler::mom_error().

Referenced by wave_structure().

619  real, dimension(:), intent(in) :: a !< lower diagonal with first entry equal to zero.
620  real, dimension(:), intent(in) :: b !< middle diagonal.
621  real, dimension(:), intent(in) :: c !< upper diagonal with last entry equal to zero.
622  real, dimension(:), intent(in) :: h !< vector of values that have already been added to b; used
623  !! for systems of the form (e.g. average layer thickness in vertical diffusion case):
624  !! [ -alpha(k-1/2) ] * e(k-1) +
625  !! [ alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k) +
626  !! [ -alpha(k+1/2) ] * e(k+1) = y(k)
627  !! where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)],
628  !! and c(k)=[-alpha(k+1/2)]. Only used with TDMA_H method.
629  real, dimension(:), intent(in) :: y !< vector of known values on right hand side.
630  character(len=*), intent(in) :: method
631  real, dimension(:), intent(out) :: x !< vector of unknown values to solve for.
632 
633 ! This subroutine solves a tri-diagonal system Ax=y using either the standard
634 ! Thomas algorithim (TDMA_T) or its more stable variant that invokes the
635 ! "Hallberg substitution" (TDMA_H).
636 !
637 ! Arguments:
638 ! (in) a - lower diagonal with first entry equal to zero
639 ! (in) b - middle diagonal
640 ! (in) c - upper diagonal with last entry equal to zero
641 ! (in) h - vector of values that have already been added to b; used for
642 ! systems of the form (e.g. average layer thickness in vertical diffusion case):
643 ! [ -alpha(k-1/2) ] * e(k-1) +
644 ! [ alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k) +
645 ! [ -alpha(k+1/2) ] * e(k+1) = y(k)
646 ! where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)],
647 ! and c(k)=[-alpha(k+1/2)]. Only used with TDMA_H method.
648 ! (in) y - vector of known values on right hand side
649 ! (out) x - vector of unknown values to solve for
650 
651  integer :: nrow ! number of rows in A matrix
652  real, allocatable, dimension(:,:) :: a_check ! for solution checking
653  real, allocatable, dimension(:) :: y_check ! for solution checking
654  real, allocatable, dimension(:) :: c_prime, y_prime, q, alpha
655  ! intermediate values for solvers
656  real :: q_prime, beta ! intermediate values for solver
657  integer :: k ! row (e.g. interface) index
658  integer :: i,j
659 
660  nrow = size(y)
661  allocate(c_prime(nrow))
662  allocate(y_prime(nrow))
663  allocate(q(nrow))
664  allocate(alpha(nrow))
665  allocate(a_check(nrow,nrow))
666  allocate(y_check(nrow))
667 
668  if (method == 'TDMA_T') then
669  ! Standard Thomas algoritim (4th variant).
670  ! Note: Requires A to be non-singular for accuracy/stability
671  c_prime(:) = 0.0 ; y_prime(:) = 0.0
672  c_prime(1) = c(1)/b(1) ; y_prime(1) = y(1)/b(1)
673 
674  ! Forward sweep
675  do k=2,nrow-1
676  c_prime(k) = c(k)/(b(k)-a(k)*c_prime(k-1))
677  enddo
678  !print *, 'c_prime=', c_prime(1:nrow)
679  do k=2,nrow
680  y_prime(k) = (y(k)-a(k)*y_prime(k-1))/(b(k)-a(k)*c_prime(k-1))
681  enddo
682  !print *, 'y_prime=', y_prime(1:nrow)
683  x(nrow) = y_prime(nrow)
684 
685  ! Backward sweep
686  do k=nrow-1,1,-1
687  x(k) = y_prime(k)-c_prime(k)*x(k+1)
688  enddo
689  !print *, 'x=',x(1:nrow)
690 
691  ! Check results - delete later
692  !do j=1,nrow ; do i=1,nrow
693  ! if(i==j)then ; A_check(i,j) = b(i)
694  ! elseif(i==j+1)then ; A_check(i,j) = a(i)
695  ! elseif(i==j-1)then ; A_check(i,j) = c(i)
696  ! endif
697  !enddo ; enddo
698  !print *, 'A(2,1),A(2,2),A(1,2)=', A_check(2,1), A_check(2,2), A_check(1,2)
699  !y_check = matmul(A_check,x)
700  !if(all(y_check .ne. y))then
701  ! print *, "tridiag_solver: Uh oh, something's not right!"
702  ! print *, "y=", y
703  ! print *, "y_check=", y_check
704  !endif
705 
706  elseif (method == 'TDMA_H') then
707  ! Thomas algoritim (4th variant) w/ Hallberg substitution.
708  ! For a layered system where k is at interfaces, alpha{k+1/2} refers to
709  ! some property (e.g. inverse thickness for mode-structure problem) of the
710  ! layer below and alpha{k-1/2} refers to the layer above.
711  ! Here, alpha(k)=alpha{k+1/2} and alpha(k-1)=alpha{k-1/2}.
712  ! Strictly speaking, this formulation requires A to be a non-singular,
713  ! symmetric, diagonally dominant matrix, with h>0.
714  ! Need to add a check for these conditions.
715  do k=1,nrow-1
716  if (abs(a(k+1)-c(k)) > 1.e-10) then
717  call mom_error(warning, "tridiag_solver: matrix not symmetric; need symmetry when invoking TDMA_H")
718  endif
719  enddo
720  alpha = -c
721  ! Alpha of the bottom-most layer is not necessarily zero. Therefore,
722  ! back out the value from the provided b(nrow and h(nrow) values
723  alpha(nrow) = b(nrow)-h(nrow)-alpha(nrow-1)
724  ! Prime other variables
725  beta = 1/b(1)
726  y_prime(:) = 0.0 ; q(:) = 0.0
727  y_prime(1) = beta*y(1) ; q(1) = beta*alpha(1)
728  q_prime = 1-q(1)
729 
730  ! Forward sweep
731  do k=2,nrow-1
732  beta = 1/(h(k)+alpha(k-1)*q_prime+alpha(k))
733  if(isnan(beta))then ; print *, "Tridiag_solver: beta is NAN" ; endif
734  q(k) = beta*alpha(k)
735  y_prime(k) = beta*(y(k)+alpha(k-1)*y_prime(k-1))
736  q_prime = beta*(h(k)+alpha(k-1)*q_prime)
737  enddo
738  if((h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow)) == 0.0)then
739  call mom_error(warning, "Tridiag_solver: this system is not stable; overriding beta(nrow).")
740  beta = 1/(1e-15) ! place holder for unstable systems - delete later
741  else
742  beta = 1/(h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow))
743  endif
744  y_prime(nrow) = beta*(y(nrow)+alpha(nrow-1)*y_prime(nrow-1))
745  x(nrow) = y_prime(nrow)
746  ! Backward sweep
747  do k=nrow-1,1,-1
748  x(k) = y_prime(k)+q(k)*x(k+1)
749  enddo
750  !print *, 'yprime=',y_prime(1:nrow)
751  !print *, 'x=',x(1:nrow)
752  endif
753 
754  deallocate(c_prime,y_prime,q,alpha,a_check,y_check)
755 
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,
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.

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

Definition at line 94 of file MOM_wave_structure.F90.

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

Referenced by mom_internal_tides::propagate_int_tide().

94  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
95  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid
96  !! structure.
97  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H
98  !! (usually m or kg m-2).
99  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
100  !! thermodynamic variables.
101  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: cn !< The (non-rotational) mode
102  !! internal gravity wave speed,
103  !! in m s-1.
104  integer, intent(in) :: modenum !< Mode number
105  real, intent(in) :: freq !< Intrinsic wave frequency, in s-1.
106  type(wave_structure_cs), pointer :: cs !< The control structure returned
107  !! by a previous call to
108  !! wave_structure_init.
109  real, dimension(SZI_(G),SZJ_(G)), &
110  optional, intent(in) :: en !< Internal wave energy density,
111  !! in Jm-2.
112  logical,optional, intent(in) :: full_halos !< If true, do the calculation
113  !! over the entire computational
114  !! domain.
115 
116 ! This subroutine determines the internal wave velocity structure for any mode.
117 ! Arguments: h - Layer thickness, in m or kg m-2.
118 ! (in) tv - A structure containing the thermobaric variables.
119 ! (in) G - The ocean's grid structure.
120 ! (in) GV - The ocean's vertical grid structure.
121 ! (in) cn - The (non-rotational) mode internal gravity wave speed, in m s-1.
122 ! (in) ModeNum - mode number
123 ! (in) freq - intrinsic wave frequency, in s-1
124 ! (in) CS - The control structure returned by a previous call to
125 ! wave_structure_init.
126 ! (in,opt) En - Internal wave energy density, in Jm-2
127 ! (in,opt) full_halos - If true, do the calculation over the entire
128 ! computational domain.
129 !
130 ! This subroutine solves for the eigen vector [vertical structure, e(k)] associated with
131 ! the first baroclinic mode speed [i.e., smallest eigen value (lam = 1/c^2)] of the
132 ! system d2e/dz2 = -(N2/cn2)e, or (A-lam*I)e = 0, where A = -(1/N2)(d2/dz2), lam = 1/c^2,
133 ! and I is the identity matrix. 2nd order discretization in the vertical lets this system
134 ! be represented as
135 !
136 ! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0
137 !
138 ! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving
139 !
140 ! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0
141 ! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0
142 !
143 ! where, upon noting N2 = reduced gravity/layer thickness, we get
144 ! Igl(k) = 1.0/(gprime(k)*H(k)) ; Igu(k) = 1.0/(gprime(k)*H(k-1))
145 !
146 ! The eigen value for this system is approximated using "wave_speed." This subroutine uses
147 ! these eigen values (mode speeds) to estimate the corresponding eigen vectors (velocity
148 ! structure) using the "inverse iteration with shift" method. The algorithm is
149 !
150 ! Pick a starting vector reasonably close to mode structure and with unit magnitude, b_guess
151 ! For n=1,2,3,...
152 ! Solve (A-lam*I)e = e_guess for e
153 ! Set e_guess=e/|e| and repeat, with each iteration refining the estimate of e
154 
155  real, dimension(SZK_(G)+1) :: &
156  drho_dt, drho_ds, &
157  pres, t_int, s_int, &
158  gprime ! The reduced gravity across each interface, in m s-2.
159  real, dimension(SZK_(G)) :: &
160  igl, igu ! The inverse of the reduced gravity across an interface times
161  ! the thickness of the layer below (Igl) or above (Igu) it,
162  ! in units of s2 m-2.
163  real, dimension(SZK_(G),SZI_(G)) :: &
164  hf, tf, sf, rf
165  real, dimension(SZK_(G)) :: &
166  hc, tc, sc, rc, &
167  det, ddet
168  real, dimension(SZI_(G),SZJ_(G)) :: &
169  htot
170  real :: lam
171  real :: min_h_frac
172  real :: h_to_pres
173  real :: h_to_m ! Local copy of a unit conversion factor.
174  real, dimension(SZI_(G)) :: &
175  hmin, & ! Thicknesses in m.
176  h_here, hxt_here, hxs_here, hxr_here
177  real :: speed2_tot
178  real :: i_hnew, drxh_sum
179  real, parameter :: tol1 = 0.0001, tol2 = 0.001
180  real, pointer, dimension(:,:,:) :: t, s
181  real :: g_rho0 ! G_Earth/Rho0 in m4 s-2 kg-1.
182  real :: rescale, i_rescale
183  integer :: kf(szi_(g))
184  integer, parameter :: max_itt = 1 ! number of times to iterate in solving for eigenvector
185  real, parameter :: cg_subro = 1e-100 ! a very small number
186  real, parameter :: a_int = 0.5 ! value of normalized integral: \int(w_strct^2)dz = a_int
187  real :: i_a_int ! inverse of a_int
188  real :: f2 ! squared Coriolis frequency
189  real :: kmag2 ! magnitude of horizontal wave number squared
190  logical :: use_eos ! If true, density is calculated from T & S using an
191  ! equation of state.
192  real, dimension(SZK_(G)+1) :: w_strct, u_strct, w_profile, uavg_profile, z_int, n2
193  ! local representations of variables in CS; note,
194  ! not all rows will be filled if layers get merged!
195  real, dimension(SZK_(G)+1) :: w_strct2, u_strct2
196  ! squared values
197  real, dimension(SZK_(G)) :: dz ! thicknesses of merged layers (same as Hc I hope)
198  real, dimension(SZK_(G)+1) :: dwdz_profile ! profile of dW/dz
199  real :: w2avg ! average of squared vertical velocity structure funtion
200  real :: int_dwdz2, int_w2, int_n2w2, ke_term, pe_term, w0
201  ! terms in vertically averaged energy equation
202  real, dimension(SZK_(G)-1) :: lam_z ! product of eigen value and gprime(k); one value for each
203  ! interface (excluding surface and bottom)
204  real, dimension(SZK_(G)-1) :: a_diag, b_diag, c_diag
205  ! diagonals of tridiagonal matrix; one value for each
206  ! interface (excluding surface and bottom)
207  real, dimension(SZK_(G)-1) :: e_guess ! guess at eigen vector with unit amplitde (for TDMA)
208  real, dimension(SZK_(G)-1) :: e_itt ! improved guess at eigen vector (from TDMA)
209  real :: pi
210  integer :: kc
211  integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop
212 
213  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
214  i_a_int = 1/a_int
215 
216  !if (present(CS)) then
217  if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_structure: "// &
218  "Module must be initialized before it is used.")
219  !endif
220 
221  if (present(full_halos)) then ; if (full_halos) then
222  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
223  endif ; endif
224 
225  pi = (4.0*atan(1.0))
226 
227  s => tv%S ; t => tv%T
228  g_rho0 = gv%g_Earth/gv%Rho0
229  use_eos = associated(tv%eqn_of_state)
230 
231  h_to_pres = gv%g_Earth * gv%Rho0
232  h_to_m = gv%H_to_m
233  rescale = 1024.0**4 ; i_rescale = 1.0/rescale
234 
235  min_h_frac = tol1 / real(nz)
236 
237  do j=js,je
238  ! First merge very thin layers with the one above (or below if they are
239  ! at the top). This also transposes the row order so that columns can
240  ! be worked upon one at a time.
241  do i=is,ie ; htot(i,j) = 0.0 ; enddo
242  do k=1,nz ; do i=is,ie ; htot(i,j) = htot(i,j) + h(i,j,k)*h_to_m ; enddo ; enddo
243 
244  do i=is,ie
245  hmin(i) = htot(i,j)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
246  hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
247  enddo
248  if (use_eos) then
249  do k=1,nz ; do i=is,ie
250  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*h_to_m > hmin(i))) then
251  hf(kf(i),i) = h_here(i)
252  tf(kf(i),i) = hxt_here(i) / h_here(i)
253  sf(kf(i),i) = hxs_here(i) / h_here(i)
254  kf(i) = kf(i) + 1
255 
256  ! Start a new layer
257  h_here(i) = h(i,j,k)*h_to_m
258  hxt_here(i) = (h(i,j,k)*h_to_m)*t(i,j,k)
259  hxs_here(i) = (h(i,j,k)*h_to_m)*s(i,j,k)
260  else
261  h_here(i) = h_here(i) + h(i,j,k)*h_to_m
262  hxt_here(i) = hxt_here(i) + (h(i,j,k)*h_to_m)*t(i,j,k)
263  hxs_here(i) = hxs_here(i) + (h(i,j,k)*h_to_m)*s(i,j,k)
264  endif
265  enddo ; enddo
266  do i=is,ie ; if (h_here(i) > 0.0) then
267  hf(kf(i),i) = h_here(i)
268  tf(kf(i),i) = hxt_here(i) / h_here(i)
269  sf(kf(i),i) = hxs_here(i) / h_here(i)
270  endif ; enddo
271  else
272  do k=1,nz ; do i=is,ie
273  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*h_to_m > hmin(i))) then
274  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
275  kf(i) = kf(i) + 1
276 
277  ! Start a new layer
278  h_here(i) = h(i,j,k)*h_to_m
279  hxr_here(i) = (h(i,j,k)*h_to_m)*gv%Rlay(k)
280  else
281  h_here(i) = h_here(i) + h(i,j,k)*h_to_m
282  hxr_here(i) = hxr_here(i) + (h(i,j,k)*h_to_m)*gv%Rlay(k)
283  endif
284  enddo ; enddo
285  do i=is,ie ; if (h_here(i) > 0.0) then
286  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
287  endif ; enddo
288  endif ! use_EOS?
289 
290  ! From this point, we can work on individual columns without causing memory
291  ! to have page faults.
292  do i=is,ie ; if(cn(i,j)>0.0)then
293  !----for debugging, remove later----
294  ig = i + g%idg_offset ; jg = j + g%jdg_offset
295  !if(ig .eq. CS%int_tide_source_x .and. jg .eq. CS%int_tide_source_y) then
296  !-----------------------------------
297  if (g%mask2dT(i,j) > 0.5) then
298 
299  lam = 1/(cn(i,j)**2)
300 
301  ! Calculate drxh_sum
302  if (use_eos) then
303  pres(1) = 0.0
304  do k=2,kf(i)
305  pres(k) = pres(k-1) + h_to_pres*hf(k-1,i)
306  t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
307  s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
308  enddo
309  call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, 2, &
310  kf(i)-1, tv%eqn_of_state)
311 
312  ! Sum the reduced gravities to find out how small a density difference
313  ! is negligibly small.
314  drxh_sum = 0.0
315  do k=2,kf(i)
316  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
317  max(0.0,drho_dt(k)*(tf(k,i)-tf(k-1,i)) + &
318  drho_ds(k)*(sf(k,i)-sf(k-1,i)))
319  enddo
320  else
321  drxh_sum = 0.0
322  do k=2,kf(i)
323  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
324  max(0.0,rf(k,i)-rf(k-1,i))
325  enddo
326  endif ! use_EOS?
327 
328  ! Find gprime across each internal interface, taking care of convective
329  ! instabilities by merging layers.
330  if (drxh_sum >= 0.0) then
331  ! Merge layers to eliminate convective instabilities or exceedingly
332  ! small reduced gravities.
333  if (use_eos) then
334  kc = 1
335  hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
336  do k=2,kf(i)
337  if ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
338  (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum) then
339  ! Merge this layer with the one above and backtrack.
340  i_hnew = 1.0 / (hc(kc) + hf(k,i))
341  tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
342  sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
343  hc(kc) = (hc(kc) + hf(k,i))
344  ! Backtrack to remove any convective instabilities above... Note
345  ! that the tolerance is a factor of two larger, to avoid limit how
346  ! far back we go.
347  do k2=kc,2,-1
348  if ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
349  (hc(k2) + hc(k2-1)) < tol2*drxh_sum) then
350  ! Merge the two bottommost layers. At this point kc = k2.
351  i_hnew = 1.0 / (hc(kc) + hc(kc-1))
352  tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
353  sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
354  hc(kc-1) = (hc(kc) + hc(kc-1))
355  kc = kc - 1
356  else ; exit ; endif
357  enddo
358  else
359  ! Add a new layer to the column.
360  kc = kc + 1
361  drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
362  tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
363  endif
364  enddo
365  ! At this point there are kc layers and the gprimes should be positive.
366  do k=2,kc ! Revisit this if non-Boussinesq.
367  gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + &
368  drho_ds(k)*(sc(k)-sc(k-1)))
369  enddo
370  else ! .not.use_EOS
371  ! Do the same with density directly...
372  kc = 1
373  hc(1) = hf(1,i) ; rc(1) = rf(1,i)
374  do k=2,kf(i)
375  if ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum) then
376  ! Merge this layer with the one above and backtrack.
377  rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
378  hc(kc) = (hc(kc) + hf(k,i))
379  ! Backtrack to remove any convective instabilities above... Note
380  ! that the tolerance is a factor of two larger, to avoid limit how
381  ! far back we go.
382  do k2=kc,2,-1
383  if ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum) then
384  ! Merge the two bottommost layers. At this point kc = k2.
385  rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
386  hc(kc-1) = (hc(kc) + hc(kc-1))
387  kc = kc - 1
388  else ; exit ; endif
389  enddo
390  else
391  ! Add a new layer to the column.
392  kc = kc + 1
393  rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
394  endif
395  enddo
396  ! At this point there are kc layers and the gprimes should be positive.
397  do k=2,kc ! Revisit this if non-Boussinesq.
398  gprime(k) = g_rho0 * (rc(k)-rc(k-1))
399  enddo
400  endif ! use_EOS?
401 
402  !-----------------NOW FIND WAVE STRUCTURE-------------------------------------
403  ! Construct and solve tridiagonal system for the interior interfaces
404  ! Note that kc = number of layers,
405  ! kc+1 = nzm = number of interfaces,
406  ! kc-1 = number of interior interfaces (excluding surface and bottom)
407  ! Also, note that "K" refers to an interface, while "k" refers to the layer below.
408  ! Need at least 3 layers (2 internal interfaces) to generate a matrix, also
409  ! need number of layers to be greater than the mode number
410  if (kc >= modenum + 1) then
411  ! Set depth at surface
412  z_int(1) = 0.0
413  ! Calculate Igu, Igl, depth, and N2 at each interior interface
414  ! [excludes surface (K=1) and bottom (K=kc+1)]
415  do k=2,kc
416  igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
417  z_int(k) = z_int(k-1) + hc(k-1)
418  n2(k) = gprime(k)/(0.5*(hc(k)+hc(k-1)))
419  enddo
420  ! Set stratification for surface and bottom (setting equal to nearest interface for now)
421  n2(1) = n2(2) ; n2(kc+1) = n2(kc)
422  ! Calcualte depth at bottom
423  z_int(kc+1) = z_int(kc)+hc(kc)
424  ! check that thicknesses sum to total depth
425  if (abs(z_int(kc+1)-htot(i,j)) > 1.e-10) then
426  call mom_error(warning, "wave_structure: mismatch in total depths")
427  print *, "kc=", kc
428  print *, "z_int(kc+1)=", z_int(kc+1)
429  print *, "htot(i,j)=", htot(i,j)
430  endif
431 
432  ! Populate interior rows of tridiagonal matrix; must multiply through by
433  ! gprime to get tridiagonal matrix to the symmetrical form:
434  ! [-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,
435  ! where lam_z = lam*gprime is now a function of depth.
436  ! Frist, populate interior rows
437  do k=3,kc-1
438  row = k-1 ! indexing for TD matrix rows
439  lam_z(row) = lam*gprime(k)
440  a_diag(row) = gprime(k)*(-igu(k))
441  b_diag(row) = gprime(k)*(igu(k)+igl(k)) - lam_z(row)
442  c_diag(row) = gprime(k)*(-igl(k))
443  if(isnan(lam_z(row)))then ; print *, "Wave_structure: lam_z(row) is NAN" ; endif
444  if(isnan(a_diag(row)))then ; print *, "Wave_structure: a(k) is NAN" ; endif
445  if(isnan(b_diag(row)))then ; print *, "Wave_structure: b(k) is NAN" ; endif
446  if(isnan(c_diag(row)))then ; print *, "Wave_structure: c(k) is NAN" ; endif
447  enddo
448  ! Populate top row of tridiagonal matrix
449  k=2 ; row = k-1
450  lam_z(row) = lam*gprime(k)
451  a_diag(row) = 0.0
452  b_diag(row) = gprime(k)*(igu(k)+igl(k)) - lam_z(row)
453  c_diag(row) = gprime(k)*(-igl(k))
454  ! Populate bottom row of tridiagonal matrix
455  k=kc ; row = k-1
456  lam_z(row) = lam*gprime(k)
457  a_diag(row) = gprime(k)*(-igu(k))
458  b_diag(row) = gprime(k)*(igu(k)+igl(k)) - lam_z(row)
459  c_diag(row) = 0.0
460 
461  ! Guess a vector shape to start with (excludes surface and bottom)
462  e_guess(1:kc-1) = sin(z_int(2:kc)/htot(i,j)*pi)
463  e_guess(1:kc-1) = e_guess(1:kc-1)/sqrt(sum(e_guess(1:kc-1)**2))
464 
465  ! Perform inverse iteration with tri-diag solver
466  do itt=1,max_itt
467  call tridiag_solver(a_diag(1:kc-1),b_diag(1:kc-1),c_diag(1:kc-1), &
468  -lam_z(1:kc-1),e_guess(1:kc-1),"TDMA_H",e_itt)
469  e_guess(1:kc-1) = e_itt(1:kc-1)/sqrt(sum(e_itt(1:kc-1)**2))
470  enddo ! itt-loop
471  w_strct(2:kc) = e_guess(1:kc-1)
472  w_strct(1) = 0.0 ! rigid lid at surface
473  w_strct(kc+1) = 0.0 ! zero-flux at bottom
474 
475  ! Check to see if solver worked
476  ig_stop = 0 ; jg_stop = 0
477  if(isnan(sum(w_strct(1:kc+1))))then
478  print *, "Wave_structure: w_strct has a NAN at ig=", ig, ", jg=", jg
479  if(i<g%isc .or. i>g%iec .or. j<g%jsc .or. j>g%jec)then
480  print *, "This is occuring at a halo point."
481  endif
482  ig_stop = ig ; jg_stop = jg
483  endif
484 
485  ! Normalize vertical structure function of w such that
486  ! \int(w_strct)^2dz = a_int (a_int could be any value, e.g., 0.5)
487  nzm = kc+1 ! number of layer interfaces after merging
488  !(including surface and bottom)
489  w2avg = 0.0
490  do k=1,nzm-1
491  dz(k) = hc(k)
492  w2avg = w2avg + 0.5*(w_strct(k)**2+w_strct(k+1)**2)*dz(k)
493  enddo
494  w2avg = w2avg/htot(i,j)
495  w_strct = w_strct/sqrt(htot(i,j)*w2avg*i_a_int)
496 
497  ! Calculate vertical structure function of u (i.e. dw/dz)
498  do k=2,nzm-1
499  u_strct(k) = 0.5*((w_strct(k-1) - w_strct(k) )/dz(k-1) + &
500  (w_strct(k) - w_strct(k+1))/dz(k))
501  enddo
502  u_strct(1) = (w_strct(1) - w_strct(2) )/dz(1)
503  u_strct(nzm) = (w_strct(nzm-1)- w_strct(nzm))/dz(nzm-1)
504 
505  ! Calculate wavenumber magnitude
506  f2 = g%CoriolisBu(i,j)**2
507  !f2 = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
508  ! (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
509  kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cg_subro**2)
510 
511  ! Calculate terms in vertically integrated energy equation
512  int_dwdz2 = 0.0 ; int_w2 = 0.0 ; int_n2w2 = 0.0
513  u_strct2 = u_strct(1:nzm)**2
514  w_strct2 = w_strct(1:nzm)**2
515  ! vertical integration with Trapezoidal rule
516  do k=1,nzm-1
517  int_dwdz2 = int_dwdz2 + 0.5*(u_strct2(k)+u_strct2(k+1))*dz(k)
518  int_w2 = int_w2 + 0.5*(w_strct2(k)+w_strct2(k+1))*dz(k)
519  int_n2w2 = int_n2w2 + 0.5*(w_strct2(k)*n2(k)+w_strct2(k+1)*n2(k+1))*dz(k)
520  enddo
521 
522  ! Back-calculate amplitude from energy equation
523  if (kmag2 > 0.0) then
524  ke_term = 0.25*gv%Rho0*( (1+f2/freq**2)/kmag2*int_dwdz2 + int_w2 )
525  pe_term = 0.25*gv%Rho0*( int_n2w2/freq**2 )
526  if (en(i,j) >= 0.0) then
527  w0 = sqrt( en(i,j)/(ke_term + pe_term) )
528  else
529  call mom_error(warning, "wave_structure: En < 0.0; setting to W0 to 0.0")
530  print *, "En(i,j)=", en(i,j), " at ig=", ig, ", jg=", jg
531  w0 = 0.0
532  endif
533  ! Calculate actual vertical velocity profile and derivative
534  w_profile = w0*w_strct
535  dwdz_profile = w0*u_strct
536  ! Calculate average magnitude of actual horizontal velocity over a period
537  uavg_profile = abs(dwdz_profile) * sqrt((1+f2/freq**2)/(2.0*kmag2))
538  else
539  w_profile = 0.0
540  dwdz_profile = 0.0
541  uavg_profile = 0.0
542  endif
543 
544  ! Store values in control structure
545  cs%w_strct(i,j,1:nzm) = w_strct
546  cs%u_strct(i,j,1:nzm) = u_strct
547  cs%W_profile(i,j,1:nzm) = w_profile
548  cs%Uavg_profile(i,j,1:nzm)= uavg_profile
549  cs%z_depths(i,j,1:nzm) = z_int
550  cs%N2(i,j,1:nzm) = n2
551  cs%num_intfaces(i,j) = nzm
552 
553  !----for debugging; delete later----
554  !if(ig .eq. ig_stop .and. jg .eq. jg_stop) then
555  !print *, 'cn(ig,jg)=', cn(i,j)
556  !print *, "e_guess=", e_guess(1:kc-1)
557  !print *, "|e_guess|=", sqrt(sum(e_guess(1:kc-1)**2))
558  !print *, 'f0=', sqrt(f2)
559  !print *, 'freq=', freq
560  !print *, 'Kh=', sqrt(Kmag2)
561  !print *, 'Wave_structure: z_int(ig,jg)=', z_int(1:nzm)
562  !print *, 'Wave_structure: N2(ig,jg)=', N2(1:nzm)
563  !print *, 'gprime=', gprime(1:nzm)
564  !print *, '1/Hc=', 1/Hc
565  !print *, 'Wave_structure: a_diag(ig,jg)=', a_diag(1:kc-1)
566  !print *, 'Wave_structure: b_diag(ig,jg)=', b_diag(1:kc-1)
567  !print *, 'Wave_structure: c_diag(ig,jg)=', c_diag(1:kc-1)
568  !print *, 'Wave_structure: lam_z(ig,jg)=', lam_z(1:kc-1)
569  !print *, 'Wave_structure: w_strct(ig,jg)=', w_strct(1:nzm)
570  !print *, 'En(i,j)=', En(i,j)
571  !print *, 'Wave_structure: W_profile(ig,jg)=', W_profile(1:nzm)
572  !print *,'int_dwdz2 =',int_dwdz2
573  !print *,'int_w2 =',int_w2
574  !print *,'int_N2w2 =',int_N2w2
575  !print *,'KEterm=',KE_term
576  !print *,'PEterm=',PE_term
577  !print *, 'W0=',W0
578  !print *,'Uavg_profile=',Uavg_profile(1:nzm)
579  !open(unit=1,file='out_N2',form='formatted') ; write(1,*) N2 ; close(1)
580  !open(unit=2,file='out_z',form='formatted') ; write(2,*) z_int ; close(2)
581  !endif
582  !-----------------------------------
583 
584  else
585  ! If not enough layers, default to zero
586  nzm = kc+1
587  cs%w_strct(i,j,1:nzm) = 0.0
588  cs%u_strct(i,j,1:nzm) = 0.0
589  cs%W_profile(i,j,1:nzm) = 0.0
590  cs%Uavg_profile(i,j,1:nzm)= 0.0
591  cs%z_depths(i,j,1:nzm) = 0.0 ! could use actual values
592  cs%N2(i,j,1:nzm) = 0.0 ! could use with actual values
593  cs%num_intfaces(i,j) = nzm
594  endif ! kc >= 3 and kc > ModeNum + 1?
595  endif ! drxh_sum >= 0?
596  !else ! if at test point - delete later
597  ! return ! if at test point - delete later
598  !endif ! if at test point - delete later
599  endif ! mask2dT > 0.5?
600  else
601  ! if cn=0.0, default to zero
602  nzm = nz+1! could use actual values
603  cs%w_strct(i,j,1:nzm) = 0.0
604  cs%u_strct(i,j,1:nzm) = 0.0
605  cs%W_profile(i,j,1:nzm) = 0.0
606  cs%Uavg_profile(i,j,1:nzm)= 0.0
607  cs%z_depths(i,j,1:nzm) = 0.0 ! could use actual values
608  cs%N2(i,j,1:nzm) = 0.0 ! could use with actual values
609  cs%num_intfaces(i,j) = nzm
610  endif ; enddo ! if cn>0.0? ; i-loop
611  enddo ! j-loop
612 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ wave_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 
)
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 760 of file MOM_wave_structure.F90.

References mom_error_handler::mom_error().

760  type(time_type), intent(in) :: time !< The current model time.
761  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
762  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
763  !! parameters.
764  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
765  !! diagnostic output.
766  type(wave_structure_cs), pointer :: cs !< A pointer that is set to point to the
767  !! control structure for this module.
768 ! Arguments: Time - The current model time.
769 ! (in) G - The ocean's grid structure.
770 ! (in) param_file - A structure indicating the open file to parse for
771 ! model parameter values.
772 ! (in) diag - A structure that is used to regulate diagnostic output.
773 ! (in/out) CS - A pointer that is set to point to the control structure
774 ! for this module
775 ! This include declares and sets the variable "version".
776 #include "version_variable.h"
777  character(len=40) :: mdl = "MOM_wave_structure" ! This module's name.
778  integer :: isd, ied, jsd, jed, nz
779 
780  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
781 
782  if (associated(cs)) then
783  call mom_error(warning, "wave_structure_init called with an "// &
784  "associated control structure.")
785  return
786  else ; allocate(cs) ; endif
787 
788  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
789  "X Location of generation site for internal tide", default=1.)
790  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
791  "Y Location of generation site for internal tide", default=1.)
792 
793  cs%diag => diag
794 
795  ! Allocate memory for variable in control structure; note,
796  ! not all rows will be filled if layers get merged!
797  allocate(cs%w_strct(isd:ied,jsd:jed,nz+1))
798  allocate(cs%u_strct(isd:ied,jsd:jed,nz+1))
799  allocate(cs%W_profile(isd:ied,jsd:jed,nz+1))
800  allocate(cs%Uavg_profile(isd:ied,jsd:jed,nz+1))
801  allocate(cs%z_depths(isd:ied,jsd:jed,nz+1))
802  allocate(cs%N2(isd:ied,jsd:jed,nz+1))
803  allocate(cs%num_intfaces(isd:ied,jsd:jed))
804 
805  ! Write all relevant parameters to the model log.
806  call log_version(param_file, mdl, version, "")
807 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function: