\hypertarget{namespacemom__wave__structure}{}\section{mom\+\_\+wave\+\_\+structure Module Reference}
\label{namespacemom__wave__structure}\index{mom\+\_\+wave\+\_\+structure@{mom\+\_\+wave\+\_\+structure}}


\subsection{Detailed Description}
Vertical structure functions for first baroclinic mode wave speed. \subsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \hyperlink{structmom__wave__structure_1_1wave__structure__cs}{wave\+\_\+structure\+\_\+cs}
\begin{DoxyCompactList}\small\item\em The control structure for the M\+O\+M\+\_\+wave\+\_\+structure module. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \hyperlink{namespacemom__wave__structure_a4e0b6a0e08df15fde4d87030567b6e11}{wave\+\_\+structure} (h, tv, G, GV, US, cn, Mode\+Num, freq, CS, En, full\+\_\+halos)
\begin{DoxyCompactList}\small\item\em This subroutine determines the internal wave velocity structure for any mode. \end{DoxyCompactList}\item 
subroutine \hyperlink{namespacemom__wave__structure_ad8e6e47af44d24efcd3f3b80f4344fbd}{tridiag\+\_\+solver} (a, b, c, h, y, method, x)
\begin{DoxyCompactList}\small\item\em Solves a tri-\/diagonal system Ax=y using either the standard Thomas algorithm (T\+D\+M\+A\+\_\+T) or its more stable variant that invokes the \char`\"{}\+Hallberg substitution\char`\"{} (T\+D\+M\+A\+\_\+H). \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__wave__structure_a4dc27a0fbdbb402b9f4def03f70cfba2}{wave\+\_\+structure\+\_\+init} (Time, G, param\+\_\+file, diag, CS)
\begin{DoxyCompactList}\small\item\em Allocate memory associated with the wave structure module and read parameters. \end{DoxyCompactList}\end{DoxyCompactItemize}


\subsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespacemom__wave__structure_ad8e6e47af44d24efcd3f3b80f4344fbd}\label{namespacemom__wave__structure_ad8e6e47af44d24efcd3f3b80f4344fbd}} 
\index{mom\+\_\+wave\+\_\+structure@{mom\+\_\+wave\+\_\+structure}!tridiag\+\_\+solver@{tridiag\+\_\+solver}}
\index{tridiag\+\_\+solver@{tridiag\+\_\+solver}!mom\+\_\+wave\+\_\+structure@{mom\+\_\+wave\+\_\+structure}}
\subsubsection{\texorpdfstring{tridiag\+\_\+solver()}{tridiag\_solver()}}
{\footnotesize\ttfamily subroutine mom\+\_\+wave\+\_\+structure\+::tridiag\+\_\+solver (\begin{DoxyParamCaption}\item[{real, dimension(\+:), intent(in)}]{a,  }\item[{real, dimension(\+:), intent(in)}]{b,  }\item[{real, dimension(\+:), intent(in)}]{c,  }\item[{real, dimension(\+:), intent(in)}]{h,  }\item[{real, dimension(\+:), intent(in)}]{y,  }\item[{character(len=$\ast$), intent(in)}]{method,  }\item[{real, dimension(\+:), intent(out)}]{x }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Solves a tri-\/diagonal system Ax=y using either the standard Thomas algorithm (T\+D\+M\+A\+\_\+T) or its more stable variant that invokes the \char`\"{}\+Hallberg substitution\char`\"{} (T\+D\+M\+A\+\_\+H). 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em a} & lower diagonal with first entry equal to zero.\\
\hline
\mbox{\tt in}  & {\em b} & middle diagonal.\\
\hline
\mbox{\tt in}  & {\em c} & upper diagonal with last entry equal to zero.\\
\hline
\mbox{\tt in}  & {\em h} & vector of values that have already been added to b; used for systems of the form (e.\+g. average layer thickness in vertical diffusion case)\+: \mbox{[} -\/alpha(k-\/1/2) \mbox{]} $\ast$ e(k-\/1) + \mbox{[} alpha(k-\/1/2) + alpha(k+1/2) + h(k) \mbox{]} $\ast$ e(k) + \mbox{[} -\/alpha(k+1/2) \mbox{]} $\ast$ e(k+1) = y(k) where a(k)=\mbox{[}-\/alpha(k-\/1/2)\mbox{]}, b(k)=\mbox{[}alpha(k-\/1/2)+alpha(k+1/2) + h(k)\mbox{]}, and c(k)=\mbox{[}-\/alpha(k+1/2)\mbox{]}. Only used with T\+D\+M\+A\+\_\+H method.\\
\hline
\mbox{\tt in}  & {\em y} & vector of known values on right hand side.\\
\hline
\mbox{\tt in}  & {\em method} & A string describing the algorithm to use\\
\hline
\mbox{\tt out}  & {\em x} & vector of unknown values to solve for. \\
\hline
\end{DoxyParams}


Definition at line 563 of file M\+O\+M\+\_\+wave\+\_\+structure.\+F90.


\begin{DoxyCode}
563   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: a\textcolor{comment}{ !< lower diagonal with first entry equal to zero.}
564   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: b\textcolor{comment}{ !< middle diagonal.}
565   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: c\textcolor{comment}{ !< upper diagonal with last entry equal to zero.}
566   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{ !< vector of values that have already been added to b; used}
567 \textcolor{comment}{           !! for systems of the form (e.g. average layer thickness in vertical diffusion case):}
568 \textcolor{comment}{           !! [ -alpha(k-1/2) ]                       * e(k-1) +}
569 \textcolor{comment}{           !! [  alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k)   +}
570 \textcolor{comment}{           !! [ -alpha(k+1/2) ]                       * e(k+1) = y(k)}
571 \textcolor{comment}{           !! where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)],}
572 \textcolor{comment}{           !! and c(k)=[-alpha(k+1/2)]. Only used with TDMA\_H method.}
573   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: y\textcolor{comment}{ !< vector of known values on right hand side.}
574   \textcolor{keywordtype}{character(len=*)},   \textcolor{keywordtype}{intent(in)}  :: method\textcolor{comment}{ !< A string describing the algorithm to use}
575   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(out)} :: x\textcolor{comment}{ !< vector of unknown values to solve for.}
576   \textcolor{comment}{! Local variables}
577   \textcolor{keywordtype}{integer} :: nrow                         \textcolor{comment}{! number of rows in A matrix}
578 \textcolor{comment}{!  real, allocatable, dimension(:,:) :: A\_check ! for solution checking}
579 \textcolor{comment}{!  real, allocatable, dimension(:)   :: y\_check ! for solution checking}
580   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:)} :: c\_prime, y\_prime, q, alpha
581                                           \textcolor{comment}{! intermediate values for solvers}
582   \textcolor{keywordtype}{real}    :: q\_prime, beta                \textcolor{comment}{! intermediate values for solver}
583   \textcolor{keywordtype}{integer} :: k                            \textcolor{comment}{! row (e.g. interface) index}
584   \textcolor{keywordtype}{integer} :: i,j
585 
586   nrow = \textcolor{keyword}{size}(y)
587   \textcolor{keyword}{allocate}(c\_prime(nrow))
588   \textcolor{keyword}{allocate}(y\_prime(nrow))
589   \textcolor{keyword}{allocate}(q(nrow))
590   \textcolor{keyword}{allocate}(alpha(nrow))
591 \textcolor{comment}{!  allocate(A\_check(nrow,nrow))}
592 \textcolor{comment}{!  allocate(y\_check(nrow))}
593 
594   \textcolor{keywordflow}{if} (method == \textcolor{stringliteral}{'TDMA\_T'}) \textcolor{keywordflow}{then}
595     \textcolor{comment}{! Standard Thomas algoritim (4th variant).}
596     \textcolor{comment}{! Note: Requires A to be non-singular for accuracy/stability}
597     c\_prime(:) = 0.0       ; y\_prime(:) = 0.0
598     c\_prime(1) = c(1)/b(1) ; y\_prime(1) = y(1)/b(1)
599 
600     \textcolor{comment}{! Forward sweep}
601     \textcolor{keywordflow}{do} k=2,nrow-1
602       c\_prime(k) = c(k)/(b(k)-a(k)*c\_prime(k-1))
603 \textcolor{keywordflow}{    enddo}
604     \textcolor{comment}{!print *, 'c\_prime=', c\_prime(1:nrow)}
605     \textcolor{keywordflow}{do} k=2,nrow
606       y\_prime(k) = (y(k)-a(k)*y\_prime(k-1))/(b(k)-a(k)*c\_prime(k-1))
607 \textcolor{keywordflow}{    enddo}
608     \textcolor{comment}{!print *, 'y\_prime=', y\_prime(1:nrow)}
609     x(nrow) = y\_prime(nrow)
610 
611     \textcolor{comment}{! Backward sweep}
612     \textcolor{keywordflow}{do} k=nrow-1,1,-1
613       x(k) = y\_prime(k)-c\_prime(k)*x(k+1)
614 \textcolor{keywordflow}{    enddo}
615     \textcolor{comment}{!print *, 'x=',x(1:nrow)}
616 
617     \textcolor{comment}{! Check results - delete later}
618     \textcolor{comment}{!do j=1,nrow ; do i=1,nrow}
619     \textcolor{comment}{!  if (i==j)then ;       A\_check(i,j) = b(i)}
620     \textcolor{comment}{!  elseif (i==j+1)then ; A\_check(i,j) = a(i)}
621     \textcolor{comment}{!  elseif (i==j-1)then ; A\_check(i,j) = c(i)}
622     \textcolor{comment}{!  endif}
623     \textcolor{comment}{!enddo ; enddo}
624     \textcolor{comment}{!print *, 'A(2,1),A(2,2),A(1,2)=', A\_check(2,1), A\_check(2,2), A\_check(1,2)}
625     \textcolor{comment}{!y\_check = matmul(A\_check,x)}
626     \textcolor{comment}{!if (all(y\_check /= y))then}
627     \textcolor{comment}{!  print *, "tridiag\_solver: Uh oh, something's not right!"}
628     \textcolor{comment}{!  print *, "y=", y}
629     \textcolor{comment}{!  print *, "y\_check=", y\_check}
630     \textcolor{comment}{!endif}
631 
632   \textcolor{keywordflow}{elseif} (method == \textcolor{stringliteral}{'TDMA\_H'}) \textcolor{keywordflow}{then}
633     \textcolor{comment}{! Thomas algoritim (4th variant) w/ Hallberg substitution.}
634     \textcolor{comment}{! For a layered system where k is at interfaces, alpha\{k+1/2\} refers to}
635     \textcolor{comment}{! some property (e.g. inverse thickness for mode-structure problem) of the}
636     \textcolor{comment}{! layer below and alpha\{k-1/2\} refers to the layer above.}
637     \textcolor{comment}{! Here, alpha(k)=alpha\{k+1/2\} and alpha(k-1)=alpha\{k-1/2\}.}
638     \textcolor{comment}{! Strictly speaking, this formulation requires A to be a non-singular,}
639     \textcolor{comment}{! symmetric, diagonally dominant matrix, with h>0.}
640     \textcolor{comment}{! Need to add a check for these conditions.}
641     \textcolor{keywordflow}{do} k=1,nrow-1
642       \textcolor{keywordflow}{if} (abs(a(k+1)-c(k)) > 1.e-10*(abs(a(k+1))+abs(c(k)))) \textcolor{keywordflow}{then}
643         \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"tridiag\_solver: matrix not symmetric; need symmetry when invoking TDMA\_H"})
644 \textcolor{keywordflow}{      endif}
645 \textcolor{keywordflow}{    enddo}
646     alpha = -c
647     \textcolor{comment}{! Alpha of the bottom-most layer is not necessarily zero. Therefore,}
648     \textcolor{comment}{! back out the value from the provided b(nrow and h(nrow) values}
649     alpha(nrow)   = b(nrow)-h(nrow)-alpha(nrow-1)
650     \textcolor{comment}{! Prime other variables}
651     beta       = 1/b(1)
652     y\_prime(:) = 0.0       ; q(:) = 0.0
653     y\_prime(1) = beta*y(1) ; q(1) = beta*alpha(1)
654     q\_prime    = 1-q(1)
655 
656     \textcolor{comment}{! Forward sweep}
657     \textcolor{keywordflow}{do} k=2,nrow-1
658       beta = 1/(h(k)+alpha(k-1)*q\_prime+alpha(k))
659       \textcolor{keywordflow}{if} (isnan(beta))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Tridiag\_solver: beta is NAN"} ;\textcolor{keywordflow}{ endif}
660       q(k) = beta*alpha(k)
661       y\_prime(k) = beta*(y(k)+alpha(k-1)*y\_prime(k-1))
662       q\_prime = beta*(h(k)+alpha(k-1)*q\_prime)
663 \textcolor{keywordflow}{    enddo}
664     \textcolor{keywordflow}{if} ((h(nrow)+alpha(nrow-1)*q\_prime+alpha(nrow)) == 0.0)\textcolor{keywordflow}{then}
665       \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"Tridiag\_solver: this system is not stable."}) \textcolor{comment}{! ; overriding beta(nrow)}
666       \textcolor{comment}{! This has hard-coded dimensions: beta = 1/(1e-15) ! place holder for unstable systems - delete later}
667     \textcolor{keywordflow}{else}
668       beta = 1/(h(nrow)+alpha(nrow-1)*q\_prime+alpha(nrow))
669 \textcolor{keywordflow}{    endif}
670     y\_prime(nrow) = beta*(y(nrow)+alpha(nrow-1)*y\_prime(nrow-1))
671     x(nrow) = y\_prime(nrow)
672     \textcolor{comment}{! Backward sweep}
673     \textcolor{keywordflow}{do} k=nrow-1,1,-1
674       x(k) = y\_prime(k)+q(k)*x(k+1)
675 \textcolor{keywordflow}{    enddo}
676     \textcolor{comment}{!print *, 'yprime=',y\_prime(1:nrow)}
677     \textcolor{comment}{!print *, 'x=',x(1:nrow)}
678 \textcolor{keywordflow}{  endif}
679 
680   \textcolor{keyword}{deallocate}(c\_prime,y\_prime,q,alpha)
681 \textcolor{comment}{! deallocate(A\_check,y\_check)}
682 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__wave__structure_a4e0b6a0e08df15fde4d87030567b6e11}\label{namespacemom__wave__structure_a4e0b6a0e08df15fde4d87030567b6e11}} 
\index{mom\+\_\+wave\+\_\+structure@{mom\+\_\+wave\+\_\+structure}!wave\+\_\+structure@{wave\+\_\+structure}}
\index{wave\+\_\+structure@{wave\+\_\+structure}!mom\+\_\+wave\+\_\+structure@{mom\+\_\+wave\+\_\+structure}}
\subsubsection{\texorpdfstring{wave\+\_\+structure()}{wave\_structure()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+structure\+::wave\+\_\+structure (\begin{DoxyParamCaption}\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed), intent(in)}]{cn,  }\item[{integer, intent(in)}]{Mode\+Num,  }\item[{real, intent(in)}]{freq,  }\item[{type(\hyperlink{structmom__wave__structure_1_1wave__structure__cs}{wave\+\_\+structure\+\_\+cs}), pointer}]{CS,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed), intent(in), optional}]{En,  }\item[{logical, intent(in), optional}]{full\+\_\+halos }\end{DoxyParamCaption})}



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

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

-\/\+Igu(k)$\ast$e(k-\/1) + (Igu(k)+\+Igl(k)-\/lam)$\ast$e(k) -\/ Igl(k)$\ast$e(k+1) = 0.\+0

with rigid lid boundary conditions e(1) = e(nz+1) = 0.\+0 giving

(Igu(2)+\+Igl(2)-\/lam)$\ast$e(2) -\/ Igl(2)$\ast$e(3) = 0.\+0 -\/\+Igu(nz)$\ast$e(nz-\/1) + (Igu(nz)+\+Igl(nz)-\/lam)$\ast$e(nz) = 0.\+0

where, upon noting N2 = reduced gravity/layer thickness, we get Igl(k) = 1.\+0/(gprime(k)$\ast$H(k)) ; Igu(k) = 1.\+0/(gprime(k)$\ast$H(k-\/1))

The eigen value for this system is approximated using \char`\"{}wave\+\_\+speed.\char`\"{} This subroutine uses these eigen values (mode speeds) to estimate the corresponding eigen vectors (velocity structure) using the \char`\"{}inverse iteration with shift\char`\"{} 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$\ast$I)e = e\+\_\+guess for e Set e\+\_\+guess=e/$\vert$e$\vert$ and repeat, with each iteration refining the estimate of e


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure.\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure.\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & A structure pointing to various thermodynamic variables.\\
\hline
\mbox{\tt in}  & {\em cn} & The (non-\/rotational) mode internal gravity wave speed \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em modenum} & Mode number\\
\hline
\mbox{\tt in}  & {\em freq} & Intrinsic wave frequency \mbox{[}T-\/1 $\sim$$>$ s-\/1\mbox{]}.\\
\hline
 & {\em cs} & The control structure returned by a previous call to wave\+\_\+structure\+\_\+init.\\
\hline
\mbox{\tt in}  & {\em en} & Internal wave energy density \mbox{[}R Z3 T-\/2 $\sim$$>$ J m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em full\+\_\+halos} & If true, do the calculation over the entire computational domain. \\
\hline
\end{DoxyParams}


Definition at line 92 of file M\+O\+M\+\_\+wave\+\_\+structure.\+F90.


\begin{DoxyCode}
92   \textcolor{keywordtype}{type}(ocean\_grid\_type),                    \textcolor{keywordtype}{intent(in)}  :: g\textcolor{comment}{  !< The ocean's grid structure.}
93   \textcolor{keywordtype}{type}(verticalgrid\_type),                  \textcolor{keywordtype}{intent(in)}  :: gv\textcolor{comment}{ !< The ocean's vertical grid structure.}
94   \textcolor{keywordtype}{type}(unit\_scale\_type),                    \textcolor{keywordtype}{intent(in)}  :: us\textcolor{comment}{ !< A dimensional unit scaling type}
95   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{  !< Layer thicknesses [H ~> m or kg m-2]}
96   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),                    \textcolor{keywordtype}{intent(in)}  :: tv\textcolor{comment}{ !< A structure pointing to various}
97 \textcolor{comment}{                                                              !! thermodynamic variables.}
98   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))},         \textcolor{keywordtype}{intent(in)}  :: cn\textcolor{comment}{ !< The (non-rotational) mode internal}
99 \textcolor{comment}{                                                              !! gravity wave speed [L T-1 ~> m s-1].}
100   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}  :: modenum\textcolor{comment}{ !< Mode number}
101   \textcolor{keywordtype}{real},                                     \textcolor{keywordtype}{intent(in)}  :: freq\textcolor{comment}{ !< Intrinsic wave frequency [T-1 ~> s-1].}
102   \textcolor{keywordtype}{type}(wave\_structure\_cs),                  \textcolor{keywordtype}{pointer}     :: cs\textcolor{comment}{ !< The control structure returned by a}
103 \textcolor{comment}{                                                              !! previous call to wave\_structure\_init.}
104   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, &
105                                   \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: en\textcolor{comment}{ !< Internal wave energy density [R Z3 T-2 ~>
       J m-2]}
106   \textcolor{keywordtype}{logical},                        \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: full\_halos\textcolor{comment}{ !< If true, do the calculation}
107 \textcolor{comment}{                                                              !! over the entire computational domain.}
108   \textcolor{comment}{! Local variables}
109   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: &
110     drho\_dt, &    \textcolor{comment}{! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]}
111     drho\_ds, &    \textcolor{comment}{! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]}
112     pres, &       \textcolor{comment}{! Interface pressure [R L2 T-2 ~> Pa]}
113     t\_int, &      \textcolor{comment}{! Temperature interpolated to interfaces [degC]}
114     s\_int, &      \textcolor{comment}{! Salinity interpolated to interfaces [ppt]}
115     gprime        \textcolor{comment}{! The reduced gravity across each interface [m2 Z-1 s-2 ~> m s-2].}
116   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: &
117     igl, igu      \textcolor{comment}{! The inverse of the reduced gravity across an interface times}
118                   \textcolor{comment}{! the thickness of the layer below (Igl) or above (Igu) it [s2 m-2].}
119   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G),SZI\_(G))} :: &
120     hf, &         \textcolor{comment}{! Layer thicknesses after very thin layers are combined [Z ~> m]}
121     tf, &         \textcolor{comment}{! Layer temperatures after very thin layers are combined [degC]}
122     sf, &         \textcolor{comment}{! Layer salinities after very thin layers are combined [ppt]}
123     rf            \textcolor{comment}{! Layer densities after very thin layers are combined [R ~> kg m-3]}
124   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: &
125     hc, &         \textcolor{comment}{! A column of layer thicknesses after convective istabilities are removed [Z ~> m]}
126     tc, &         \textcolor{comment}{! A column of layer temperatures after convective istabilities are removed [degC]}
127     sc, &         \textcolor{comment}{! A column of layer salinites after convective istabilities are removed [ppt]}
128     rc, &         \textcolor{comment}{! A column of layer densities after convective istabilities are removed [R ~> kg m-3]}
129     det, ddet
130   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))} :: &
131     htot          \textcolor{comment}{! The vertical sum of the thicknesses [Z ~> m]}
132   \textcolor{keywordtype}{real} :: lam
133   \textcolor{keywordtype}{real} :: min\_h\_frac
134   \textcolor{keywordtype}{real} :: z\_to\_pres \textcolor{comment}{! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1]}
135   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
136     hmin, &        \textcolor{comment}{! Thicknesses [Z ~> m]}
137     h\_here, &      \textcolor{comment}{! A thickness [Z ~> m]}
138     hxt\_here, &    \textcolor{comment}{! A layer integrated temperature [degC Z ~> degC m]}
139     hxs\_here, &    \textcolor{comment}{! A layer integrated salinity [ppt Z ~> ppt m]}
140     hxr\_here       \textcolor{comment}{! A layer integrated density [R Z ~> kg m-2]}
141   \textcolor{keywordtype}{real} :: speed2\_tot
142   \textcolor{keywordtype}{real} :: i\_hnew   \textcolor{comment}{! The inverse of a new layer thickness [Z-1 ~> m-1]}
143   \textcolor{keywordtype}{real} :: drxh\_sum \textcolor{comment}{! The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-2]}
144   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: tol1  = 0.0001, tol2 = 0.001
145   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{pointer}, \textcolor{keywordtype}{dimension(:,:,:)} :: t => null(), s => null()
146   \textcolor{keywordtype}{real} :: g\_rho0  \textcolor{comment}{! G\_Earth/Rho0 in [m2 s-2 Z-1 R-1 ~> m4 s-2 kg-1].}
147   \textcolor{comment}{! real :: rescale, I\_rescale}
148   \textcolor{keywordtype}{integer} :: kf(szi\_(g))
149   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: max\_itt = 1 \textcolor{comment}{! number of times to iterate in solving for eigenvector}
150   \textcolor{keywordtype}{real} :: cg\_subro        \textcolor{comment}{! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1]}
151   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}    :: a\_int = 0.5 \textcolor{comment}{! value of normalized integral: \(\backslash\)int(w\_strct^2)dz = a\_int}
152   \textcolor{keywordtype}{real}               :: i\_a\_int     \textcolor{comment}{! inverse of a\_int}
153   \textcolor{keywordtype}{real}               :: f2          \textcolor{comment}{! squared Coriolis frequency [T-2 ~> s-2]}
154   \textcolor{keywordtype}{real}               :: kmag2       \textcolor{comment}{! magnitude of horizontal wave number squared}
155   \textcolor{keywordtype}{logical}            :: use\_eos     \textcolor{comment}{! If true, density is calculated from T & S using an}
156                                     \textcolor{comment}{! equation of state.}
157   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: w\_strct, u\_strct, w\_profile, uavg\_profile, z\_int, n2
158                                         \textcolor{comment}{! local representations of variables in CS; note,}
159                                         \textcolor{comment}{! not all rows will be filled if layers get merged!}
160   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: w\_strct2, u\_strct2
161                                         \textcolor{comment}{! squared values}
162   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))}   :: dz      \textcolor{comment}{! thicknesses of merged layers (same as Hc I hope)}
163   \textcolor{comment}{! real, dimension(SZK\_(G)+1) :: dWdz\_profile ! profile of dW/dz}
164   \textcolor{keywordtype}{real}                       :: w2avg   \textcolor{comment}{! average of squared vertical velocity structure funtion}
165   \textcolor{keywordtype}{real}                       :: int\_dwdz2
166   \textcolor{keywordtype}{real}                       :: int\_w2
167   \textcolor{keywordtype}{real}                       :: int\_n2w2
168   \textcolor{keywordtype}{real}                       :: ke\_term \textcolor{comment}{! terms in vertically averaged energy equation}
169   \textcolor{keywordtype}{real}                       :: pe\_term \textcolor{comment}{! terms in vertically averaged energy equation}
170   \textcolor{keywordtype}{real}                       :: w0      \textcolor{comment}{! A vertical velocity magnitude [Z T-1 ~> m s-1]}
171   \textcolor{keywordtype}{real}                       :: gp\_unscaled \textcolor{comment}{! A version of gprime rescaled to [m s-2].}
172   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)-1)} :: lam\_z   \textcolor{comment}{! product of eigen value and gprime(k); one value for each}
173                                         \textcolor{comment}{! interface (excluding surface and bottom)}
174   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)-1)} :: a\_diag, b\_diag, c\_diag
175                                         \textcolor{comment}{! diagonals of tridiagonal matrix; one value for each}
176                                         \textcolor{comment}{! interface (excluding surface and bottom)}
177   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)-1)} :: e\_guess \textcolor{comment}{! guess at eigen vector with unit amplitde (for TDMA)}
178   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)-1)} :: e\_itt   \textcolor{comment}{! improved guess at eigen vector (from TDMA)}
179   \textcolor{keywordtype}{real}    :: pi
180   \textcolor{keywordtype}{integer} :: kc
181   \textcolor{keywordtype}{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   \textcolor{comment}{!if (present(CS)) then}
187     \textcolor{keywordflow}{if} (.not. \textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_wave\_structure: "}// &
188            \textcolor{stringliteral}{"Module must be initialized before it is used."})
189   \textcolor{comment}{!endif}
190 
191   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(full\_halos)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (full\_halos) \textcolor{keywordflow}{then}
192     is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
193 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ 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  \textcolor{comment}{! The hard-coded value here might need to increase.}
200   use\_eos = \textcolor{keyword}{associated}(tv%eqn\_of\_state)
201 
202   \textcolor{comment}{! Simplifying the following could change answers at roundoff.}
203   z\_to\_pres = gv%Z\_to\_H * (gv%H\_to\_RZ * gv%g\_Earth)
204   \textcolor{comment}{! rescale = 1024.0**4 ; I\_rescale = 1.0/rescale}
205 
206   min\_h\_frac = tol1 / \textcolor{keywordtype}{real}(nz)
207 
208   \textcolor{keywordflow}{do} j=js,je
209     \textcolor{comment}{!   First merge very thin layers with the one above (or below if they are}
210     \textcolor{comment}{! at the top).  This also transposes the row order so that columns can}
211     \textcolor{comment}{! be worked upon one at a time.}
212     \textcolor{keywordflow}{do} i=is,ie ; htot(i,j) = 0.0 ;\textcolor{keywordflow}{ enddo}
213     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; htot(i,j) = htot(i,j) + h(i,j,k)*gv%H\_to\_Z ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
214 
215     \textcolor{keywordflow}{do} i=is,ie
216       hmin(i) = htot(i,j)*min\_h\_frac ; kf(i) = 1 ; h\_here(i) = 0.0
217       hxt\_here(i) = 0.0 ; hxs\_here(i) = 0.0 ; hxr\_here(i) = 0.0
218 \textcolor{keywordflow}{    enddo}
219     \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
220       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
221         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}
222           hf(kf(i),i) = h\_here(i)
223           tf(kf(i),i) = hxt\_here(i) / h\_here(i)
224           sf(kf(i),i) = hxs\_here(i) / h\_here(i)
225           kf(i) = kf(i) + 1
226 
227           \textcolor{comment}{! Start a new layer}
228           h\_here(i) = h(i,j,k)*gv%H\_to\_Z
229           hxt\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*t(i,j,k)
230           hxs\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*s(i,j,k)
231         \textcolor{keywordflow}{else}
232           h\_here(i) = h\_here(i) + h(i,j,k)*gv%H\_to\_Z
233           hxt\_here(i) = hxt\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*t(i,j,k)
234           hxs\_here(i) = hxs\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*s(i,j,k)
235 \textcolor{keywordflow}{        endif}
236 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
237       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}
238         hf(kf(i),i) = h\_here(i)
239         tf(kf(i),i) = hxt\_here(i) / h\_here(i)
240         sf(kf(i),i) = hxs\_here(i) / h\_here(i)
241 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
242     \textcolor{keywordflow}{else}
243       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
244         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}
245           hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)
246           kf(i) = kf(i) + 1
247 
248           \textcolor{comment}{! Start a new layer}
249           h\_here(i) = h(i,j,k)*gv%H\_to\_Z
250           hxr\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*gv%Rlay(k)
251         \textcolor{keywordflow}{else}
252           h\_here(i) = h\_here(i) + h(i,j,k)*gv%H\_to\_Z
253           hxr\_here(i) = hxr\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*gv%Rlay(k)
254 \textcolor{keywordflow}{        endif}
255 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
256       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}
257         hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)
258 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
259 \textcolor{keywordflow}{    endif} \textcolor{comment}{! use\_EOS?}
260 
261     \textcolor{comment}{! From this point, we can work on individual columns without causing memory}
262     \textcolor{comment}{! to have page faults.}
263     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (cn(i,j)>0.0)\textcolor{keywordflow}{then}
264       \textcolor{comment}{!----for debugging, remove later----}
265       ig = i + g%idg\_offset ; jg = j + g%jdg\_offset
266       \textcolor{comment}{!if (ig == CS%int\_tide\_source\_x .and. jg == CS%int\_tide\_source\_y) then}
267       \textcolor{comment}{!-----------------------------------}
268       \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then}
269 
270         lam = 1/(us%L\_T\_to\_m\_s**2 * cn(i,j)**2)
271 
272         \textcolor{comment}{! Calculate drxh\_sum}
273         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
274           pres(1) = 0.0
275           \textcolor{keywordflow}{do} k=2,kf(i)
276             pres(k) = pres(k-1) + z\_to\_pres*hf(k-1,i)
277             t\_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
278             s\_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
279 \textcolor{keywordflow}{          enddo}
280           \textcolor{keyword}{call }calculate\_density\_derivs(t\_int, s\_int, pres, drho\_dt, drho\_ds, &
281                                         tv%eqn\_of\_state, (/2,kf(i)/) )
282 
283           \textcolor{comment}{! Sum the reduced gravities to find out how small a density difference}
284           \textcolor{comment}{! is negligibly small.}
285           drxh\_sum = 0.0
286           \textcolor{keywordflow}{do} k=2,kf(i)
287             drxh\_sum = drxh\_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
288                 max(0.0,drho\_dt(k)*(tf(k,i)-tf(k-1,i)) + &
289                         drho\_ds(k)*(sf(k,i)-sf(k-1,i)))
290 \textcolor{keywordflow}{          enddo}
291         \textcolor{keywordflow}{else}
292           drxh\_sum = 0.0
293           \textcolor{keywordflow}{do} k=2,kf(i)
294             drxh\_sum = drxh\_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
295                               max(0.0,rf(k,i)-rf(k-1,i))
296 \textcolor{keywordflow}{          enddo}
297 \textcolor{keywordflow}{        endif} \textcolor{comment}{! use\_EOS?}
298 
299         \textcolor{comment}{!   Find gprime across each internal interface, taking care of convective}
300         \textcolor{comment}{! instabilities by merging layers.}
301         \textcolor{keywordflow}{if} (drxh\_sum >= 0.0) \textcolor{keywordflow}{then}
302           \textcolor{comment}{! Merge layers to eliminate convective instabilities or exceedingly}
303           \textcolor{comment}{! small reduced gravities.}
304           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
305             kc = 1
306             hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
307             \textcolor{keywordflow}{do} k=2,kf(i)
308               \textcolor{keywordflow}{if} ((drho\_dt(k)*(tf(k,i)-tc(kc)) + drho\_ds(k)*(sf(k,i)-sc(kc))) * &
309                   (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh\_sum) \textcolor{keywordflow}{then}
310                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}
311                 i\_hnew = 1.0 / (hc(kc) + hf(k,i))
312                 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i\_hnew
313                 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i\_hnew
314                 hc(kc) = (hc(kc) + hf(k,i))
315                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}
316                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}
317                 \textcolor{comment}{! far back we go.}
318                 \textcolor{keywordflow}{do} k2=kc,2,-1
319                   \textcolor{keywordflow}{if} ((drho\_dt(k2)*(tc(k2)-tc(k2-1)) + drho\_ds(k2)*(sc(k2)-sc(k2-1))) * &
320                       (hc(k2) + hc(k2-1)) < tol2*drxh\_sum) \textcolor{keywordflow}{then}
321                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}
322                     i\_hnew = 1.0 / (hc(kc) + hc(kc-1))
323                     tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i\_hnew
324                     sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i\_hnew
325                     hc(kc-1) = (hc(kc) + hc(kc-1))
326                     kc = kc - 1
327                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}
328 \textcolor{keywordflow}{                enddo}
329               \textcolor{keywordflow}{else}
330                 \textcolor{comment}{! Add a new layer to the column.}
331                 kc = kc + 1
332                 drho\_ds(kc) = drho\_ds(k) ; drho\_dt(kc) = drho\_dt(k)
333                 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
334 \textcolor{keywordflow}{              endif}
335 \textcolor{keywordflow}{            enddo}
336             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}
337             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-Boussinesq.}
338               gprime(k) = g\_rho0 * (drho\_dt(k)*(tc(k)-tc(k-1)) + &
339                                     drho\_ds(k)*(sc(k)-sc(k-1)))
340 \textcolor{keywordflow}{            enddo}
341           \textcolor{keywordflow}{else}  \textcolor{comment}{! .not.use\_EOS}
342             \textcolor{comment}{! Do the same with density directly...}
343             kc = 1
344             hc(1) = hf(1,i) ; rc(1) = rf(1,i)
345             \textcolor{keywordflow}{do} k=2,kf(i)
346               \textcolor{keywordflow}{if} ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh\_sum) \textcolor{keywordflow}{then}
347                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}
348                 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
349                 hc(kc) = (hc(kc) + hf(k,i))
350                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}
351                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}
352                 \textcolor{comment}{! far back we go.}
353                 \textcolor{keywordflow}{do} k2=kc,2,-1
354                   \textcolor{keywordflow}{if} ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh\_sum) \textcolor{keywordflow}{then}
355                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}
356                     rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
357                     hc(kc-1) = (hc(kc) + hc(kc-1))
358                     kc = kc - 1
359                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}
360 \textcolor{keywordflow}{                enddo}
361               \textcolor{keywordflow}{else}
362                 \textcolor{comment}{! Add a new layer to the column.}
363                 kc = kc + 1
364                 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
365 \textcolor{keywordflow}{              endif}
366 \textcolor{keywordflow}{            enddo}
367             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}
368             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-Boussinesq.}
369               gprime(k) = g\_rho0 * (rc(k)-rc(k-1))
370 \textcolor{keywordflow}{            enddo}
371 \textcolor{keywordflow}{          endif}  \textcolor{comment}{! use\_EOS?}
372 
373           \textcolor{comment}{!-----------------NOW FIND WAVE STRUCTURE-------------------------------------}
374           \textcolor{comment}{! Construct and solve tridiagonal system for the interior interfaces}
375           \textcolor{comment}{! Note that kc   = number of layers,}
376           \textcolor{comment}{!           kc+1 = nzm = number of interfaces,}
377           \textcolor{comment}{!           kc-1 = number of interior interfaces (excluding surface and bottom)}
378           \textcolor{comment}{! Also, note that "K" refers to an interface, while "k" refers to the layer below.}
379           \textcolor{comment}{! Need at least 3 layers (2 internal interfaces) to generate a matrix, also}
380           \textcolor{comment}{! need number of layers to be greater than the mode number}
381           \textcolor{keywordflow}{if} (kc >= max(3, modenum + 1)) \textcolor{keywordflow}{then}
382             \textcolor{comment}{! Set depth at surface}
383             z\_int(1) = 0.0
384             \textcolor{comment}{! Calculate Igu, Igl, depth, and N2 at each interior interface}
385             \textcolor{comment}{! [excludes surface (K=1) and bottom (K=kc+1)]}
386             \textcolor{keywordflow}{do} k=2,kc
387               igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
388               z\_int(k) = z\_int(k-1) + hc(k-1)
389               n2(k) = us%m\_to\_Z**2*gprime(k)/(0.5*(hc(k)+hc(k-1)))
390 \textcolor{keywordflow}{            enddo}
391             \textcolor{comment}{! Set stratification for surface and bottom (setting equal to nearest interface for now)}
392             n2(1) = n2(2) ; n2(kc+1) = n2(kc)
393             \textcolor{comment}{! Calcualte depth at bottom}
394             z\_int(kc+1) = z\_int(kc)+hc(kc)
395             \textcolor{comment}{! check that thicknesses sum to total depth}
396             \textcolor{keywordflow}{if} (abs(z\_int(kc+1)-htot(i,j)) > 1.e-14*htot(i,j)) \textcolor{keywordflow}{then}
397               \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"wave\_structure: mismatch in total depths"})
398 \textcolor{keywordflow}{            endif}
399 
400             \textcolor{comment}{! Note that many of the calcluation from here on revert to using vertical}
401             \textcolor{comment}{! distances in m, not Z.}
402 
403             \textcolor{comment}{! Populate interior rows of tridiagonal matrix; must multiply through by}
404             \textcolor{comment}{! gprime to get tridiagonal matrix to the symmetrical form:}
405             \textcolor{comment}{! [-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,}
406             \textcolor{comment}{! where lam\_z = lam*gprime is now a function of depth.}
407             \textcolor{comment}{! Frist, populate interior rows}
408             \textcolor{keywordflow}{do} k=3,kc-1
409               row = k-1 \textcolor{comment}{! indexing for TD matrix rows}
410               gp\_unscaled = us%m\_to\_Z*gprime(k)
411               lam\_z(row) = lam*gp\_unscaled
412               a\_diag(row) = gp\_unscaled*(-igu(k))
413               b\_diag(row) = gp\_unscaled*(igu(k)+igl(k)) - lam\_z(row)
414               c\_diag(row) = gp\_unscaled*(-igl(k))
415               \textcolor{keywordflow}{if} (isnan(lam\_z(row)))\textcolor{keywordflow}{then}  ; print *, \textcolor{stringliteral}{"Wave\_structure: lam\_z(row) is NAN"} ;\textcolor{keywordflow}{ endif}
416               \textcolor{keywordflow}{if} (isnan(a\_diag(row)))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Wave\_structure: a(k) is NAN"} ;\textcolor{keywordflow}{ endif}
417               \textcolor{keywordflow}{if} (isnan(b\_diag(row)))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Wave\_structure: b(k) is NAN"} ;\textcolor{keywordflow}{ endif}
418               \textcolor{keywordflow}{if} (isnan(c\_diag(row)))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Wave\_structure: c(k) is NAN"} ;\textcolor{keywordflow}{ endif}
419 \textcolor{keywordflow}{            enddo}
420             \textcolor{comment}{! Populate top row of tridiagonal matrix}
421             k=2 ; row = k-1 ;
422             gp\_unscaled = us%m\_to\_Z*gprime(k)
423             lam\_z(row) = lam*gp\_unscaled
424             a\_diag(row) = 0.0
425             b\_diag(row) = gp\_unscaled*(igu(k)+igl(k)) - lam\_z(row)
426             c\_diag(row) = gp\_unscaled*(-igl(k))
427             \textcolor{comment}{! Populate bottom row of tridiagonal matrix}
428             k=kc ; row = k-1
429             gp\_unscaled = us%m\_to\_Z*gprime(k)
430             lam\_z(row) = lam*gp\_unscaled
431             a\_diag(row) = gp\_unscaled*(-igu(k))
432             b\_diag(row) = gp\_unscaled*(igu(k)+igl(k)) - lam\_z(row)
433             c\_diag(row) = 0.0
434 
435             \textcolor{comment}{! Guess a vector shape to start with (excludes surface and bottom)}
436             e\_guess(1:kc-1) = sin((z\_int(2:kc)/htot(i,j)) *pi)
437             e\_guess(1:kc-1) = e\_guess(1:kc-1)/sqrt(sum(e\_guess(1:kc-1)**2))
438 
439             \textcolor{comment}{! Perform inverse iteration with tri-diag solver}
440             \textcolor{keywordflow}{do} itt=1,max\_itt
441               \textcolor{keyword}{call }tridiag\_solver(a\_diag(1:kc-1),b\_diag(1:kc-1),c\_diag(1:kc-1), &
442                                   -lam\_z(1:kc-1),e\_guess(1:kc-1),\textcolor{stringliteral}{"TDMA\_H"},e\_itt)
443               e\_guess(1:kc-1) = e\_itt(1:kc-1) / sqrt(sum(e\_itt(1:kc-1)**2))
444 \textcolor{keywordflow}{            enddo} \textcolor{comment}{! itt-loop}
445             w\_strct(2:kc) = e\_guess(1:kc-1)
446             w\_strct(1)    = 0.0 \textcolor{comment}{! rigid lid at surface}
447             w\_strct(kc+1) = 0.0 \textcolor{comment}{! zero-flux at bottom}
448 
449             \textcolor{comment}{! Check to see if solver worked}
450             ig\_stop = 0 ; jg\_stop = 0
451             \textcolor{keywordflow}{if} (isnan(sum(w\_strct(1:kc+1))))\textcolor{keywordflow}{then}
452               print *, \textcolor{stringliteral}{"Wave\_structure: w\_strct has a NAN at ig="}, ig, \textcolor{stringliteral}{", jg="}, jg
453               \textcolor{keywordflow}{if} (i<g%isc .or. i>g%iec .or. j<g%jsc .or. j>g%jec)\textcolor{keywordflow}{then}
454                 print *, \textcolor{stringliteral}{"This is occuring at a halo point."}
455 \textcolor{keywordflow}{              endif}
456               ig\_stop = ig ; jg\_stop = jg
457 \textcolor{keywordflow}{            endif}
458 
459             \textcolor{comment}{! Normalize vertical structure function of w such that}
460             \textcolor{comment}{! \(\backslash\)int(w\_strct)^2dz = a\_int (a\_int could be any value, e.g., 0.5)}
461             nzm = kc+1 \textcolor{comment}{! number of layer interfaces after merging}
462                        \textcolor{comment}{!(including surface and bottom)}
463             w2avg = 0.0
464             \textcolor{keywordflow}{do} k=1,nzm-1
465               dz(k) = us%Z\_to\_m*hc(k)
466               w2avg = w2avg + 0.5*(w\_strct(k)**2+w\_strct(k+1)**2)*dz(k)
467 \textcolor{keywordflow}{            enddo}
468             \textcolor{comment}{!### Some mathematical cancellations could occur in the next two lines.}
469             w2avg = w2avg / htot(i,j)
470             w\_strct(:) = w\_strct(:) / sqrt(htot(i,j)*w2avg*i\_a\_int)
471 
472             \textcolor{comment}{! Calculate vertical structure function of u (i.e. dw/dz)}
473             \textcolor{keywordflow}{do} k=2,nzm-1
474               u\_strct(k) = 0.5*((w\_strct(k-1) - w\_strct(k)  )/dz(k-1) + &
475                                 (w\_strct(k)   - w\_strct(k+1))/dz(k))
476 \textcolor{keywordflow}{            enddo}
477             u\_strct(1)   = (w\_strct(1)   -  w\_strct(2) )/dz(1)
478             u\_strct(nzm) = (w\_strct(nzm-1)-  w\_strct(nzm))/dz(nzm-1)
479 
480             \textcolor{comment}{! Calculate wavenumber magnitude}
481             f2 = g%CoriolisBu(i,j)**2
482             \textcolor{comment}{!f2 = 0.25*US%s\_to\_T**2 *((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &}
483             \textcolor{comment}{!    (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))}
484             kmag2 = us%m\_to\_L**2 * (freq**2 - f2) / (cn(i,j)**2 + cg\_subro**2)
485 
486             \textcolor{comment}{! Calculate terms in vertically integrated energy equation}
487             int\_dwdz2 = 0.0 ; int\_w2 = 0.0 ; int\_n2w2 = 0.0
488             u\_strct2(1:nzm) = u\_strct(1:nzm)**2
489             w\_strct2(1:nzm) = w\_strct(1:nzm)**2
490             \textcolor{comment}{! vertical integration with Trapezoidal rule}
491             \textcolor{keywordflow}{do} k=1,nzm-1
492               int\_dwdz2 = int\_dwdz2 + 0.5*(u\_strct2(k)+u\_strct2(k+1)) * us%m\_to\_Z*dz(k)
493               int\_w2    = int\_w2    + 0.5*(w\_strct2(k)+w\_strct2(k+1)) * us%m\_to\_Z*dz(k)
494               int\_n2w2  = int\_n2w2  + 0.5*(w\_strct2(k)*n2(k)+w\_strct2(k+1)*n2(k+1)) * us%m\_to\_Z*dz(k)
495 \textcolor{keywordflow}{            enddo}
496 
497             \textcolor{comment}{! Back-calculate amplitude from energy equation}
498             \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(en) .and. (freq**2*kmag2 > 0.0)) \textcolor{keywordflow}{then}
499               \textcolor{comment}{! Units here are [R}
500               ke\_term = 0.25*gv%Rho0*( ((freq**2 + f2) / (freq**2*kmag2))*int\_dwdz2 + int\_w2 )
501               pe\_term = 0.25*gv%Rho0*( int\_n2w2 / (us%s\_to\_T*freq)**2 )
502               \textcolor{keywordflow}{if} (en(i,j) >= 0.0) \textcolor{keywordflow}{then}
503                 w0 = sqrt( en(i,j) / (ke\_term + pe\_term) )
504               \textcolor{keywordflow}{else}
505                 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_structure: En < 0.0; setting to W0 to 0.0"})
506                 print *, \textcolor{stringliteral}{"En(i,j)="}, en(i,j), \textcolor{stringliteral}{" at ig="}, ig, \textcolor{stringliteral}{", jg="}, jg
507                 w0 = 0.0
508 \textcolor{keywordflow}{              endif}
509               \textcolor{comment}{! Calculate actual vertical velocity profile and derivative}
510               w\_profile(:)    = w0*w\_strct(:)
511               \textcolor{comment}{! dWdz\_profile(:) = W0*u\_strct(:)}
512               \textcolor{comment}{! Calculate average magnitude of actual horizontal velocity over a period}
513               uavg\_profile(:) = us%Z\_to\_L*abs(w0*u\_strct(:)) * sqrt((freq**2 + f2) / (2.0*freq**2*kmag2))
514             \textcolor{keywordflow}{else}
515               w\_profile(:)    = 0.0
516               \textcolor{comment}{! dWdz\_profile(:) = 0.0}
517               uavg\_profile(:) = 0.0
518 \textcolor{keywordflow}{            endif}
519 
520             \textcolor{comment}{! Store values in control structure}
521             cs%w\_strct(i,j,1:nzm)     = w\_strct(1:nzm)
522             cs%u\_strct(i,j,1:nzm)     = u\_strct(1:nzm)
523             cs%W\_profile(i,j,1:nzm)   = w\_profile(1:nzm)
524             cs%Uavg\_profile(i,j,1:nzm)= uavg\_profile(1:nzm)
525             cs%z\_depths(i,j,1:nzm)    = us%Z\_to\_m*z\_int(1:nzm)
526             cs%N2(i,j,1:nzm)          = n2(1:nzm)
527             cs%num\_intfaces(i,j)      = nzm
528           \textcolor{keywordflow}{else}
529             \textcolor{comment}{! If not enough layers, default to zero}
530             nzm = kc+1
531             cs%w\_strct(i,j,1:nzm)     = 0.0
532             cs%u\_strct(i,j,1:nzm)     = 0.0
533             cs%W\_profile(i,j,1:nzm)   = 0.0
534             cs%Uavg\_profile(i,j,1:nzm)= 0.0
535             cs%z\_depths(i,j,1:nzm)    = 0.0 \textcolor{comment}{! could use actual values}
536             cs%N2(i,j,1:nzm)          = 0.0 \textcolor{comment}{! could use with actual values}
537             cs%num\_intfaces(i,j)       = nzm
538 \textcolor{keywordflow}{          endif}  \textcolor{comment}{! kc >= 3 and kc > ModeNum + 1?}
539 \textcolor{keywordflow}{        endif} \textcolor{comment}{! drxh\_sum >= 0?}
540       \textcolor{comment}{!else     ! if at test point - delete later}
541       \textcolor{comment}{!  return ! if at test point - delete later}
542       \textcolor{comment}{!endif    ! if at test point - delete later}
543 \textcolor{keywordflow}{      endif} \textcolor{comment}{! mask2dT > 0.5?}
544     \textcolor{keywordflow}{else}
545       \textcolor{comment}{! if cn=0.0, default to zero}
546       nzm                       = nz+1\textcolor{comment}{! could use actual values}
547       cs%w\_strct(i,j,1:nzm)     = 0.0
548       cs%u\_strct(i,j,1:nzm)     = 0.0
549       cs%W\_profile(i,j,1:nzm)   = 0.0
550       cs%Uavg\_profile(i,j,1:nzm)= 0.0
551       cs%z\_depths(i,j,1:nzm)    = 0.0 \textcolor{comment}{! could use actual values}
552       cs%N2(i,j,1:nzm)          = 0.0 \textcolor{comment}{! could use with actual values}
553       cs%num\_intfaces(i,j)       = nzm
554 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! if cn>0.0? ; i-loop}
555 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-loop}
556 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__wave__structure_a4dc27a0fbdbb402b9f4def03f70cfba2}\label{namespacemom__wave__structure_a4dc27a0fbdbb402b9f4def03f70cfba2}} 
\index{mom\+\_\+wave\+\_\+structure@{mom\+\_\+wave\+\_\+structure}!wave\+\_\+structure\+\_\+init@{wave\+\_\+structure\+\_\+init}}
\index{wave\+\_\+structure\+\_\+init@{wave\+\_\+structure\+\_\+init}!mom\+\_\+wave\+\_\+structure@{mom\+\_\+wave\+\_\+structure}}
\subsubsection{\texorpdfstring{wave\+\_\+structure\+\_\+init()}{wave\_structure\_init()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+structure\+::wave\+\_\+structure\+\_\+init (\begin{DoxyParamCaption}\item[{type(time\+\_\+type), intent(in)}]{Time,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(param\+\_\+file\+\_\+type), intent(in)}]{param\+\_\+file,  }\item[{type(diag\+\_\+ctrl), intent(in), target}]{diag,  }\item[{type(\hyperlink{structmom__wave__structure_1_1wave__structure__cs}{wave\+\_\+structure\+\_\+cs}), pointer}]{CS }\end{DoxyParamCaption})}



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


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em time} & The current model time.\\
\hline
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure.\\
\hline
\mbox{\tt in}  & {\em param\+\_\+file} & A structure to parse for run-\/time parameters.\\
\hline
\mbox{\tt in}  & {\em diag} & A structure that is used to regulate diagnostic output.\\
\hline
 & {\em cs} & A pointer that is set to point to the control structure for this module. \\
\hline
\end{DoxyParams}


Definition at line 687 of file M\+O\+M\+\_\+wave\+\_\+structure.\+F90.


\begin{DoxyCode}
687   \textcolor{keywordtype}{type}(time\_type),         \textcolor{keywordtype}{intent(in)} :: time\textcolor{comment}{ !< The current model time.}
688   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{    !< The ocean's grid structure.}
689   \textcolor{keywordtype}{type}(param\_file\_type),   \textcolor{keywordtype}{intent(in)} :: param\_file\textcolor{comment}{ !< A structure to parse for run-time}
690 \textcolor{comment}{                                              !! parameters.}
691   \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)} :: diag\textcolor{comment}{ !< A structure that is used to regulate}
692 \textcolor{comment}{                                              !! diagnostic output.}
693   \textcolor{keywordtype}{type}(wave\_structure\_cs), \textcolor{keywordtype}{pointer}    :: cs\textcolor{comment}{   !< A pointer that is set to point to the}
694 \textcolor{comment}{                                              !! control structure for this module.}
695 \textcolor{comment}{! This include declares and sets the variable "version".}
696 \textcolor{preprocessor}{#include "version\_variable.h"}
697 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_wave\_structure"}  \textcolor{comment}{! This module's name.}
698   \textcolor{keywordtype}{integer} :: isd, ied, jsd, jed, nz
699 
700   isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
701 
702   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}
703     \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_structure\_init called with an "}// &
704                             \textcolor{stringliteral}{"associated control structure."})
705     \textcolor{keywordflow}{return}
706   \textcolor{keywordflow}{else} ; \textcolor{keyword}{allocate}(cs) ;\textcolor{keywordflow}{ endif}
707 
708   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"INTERNAL\_TIDE\_SOURCE\_X"}, cs%int\_tide\_source\_x, &
709                  \textcolor{stringliteral}{"X Location of generation site for internal tide"}, default=1.)
710   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"INTERNAL\_TIDE\_SOURCE\_Y"}, cs%int\_tide\_source\_y, &
711                  \textcolor{stringliteral}{"Y Location of generation site for internal tide"}, default=1.)
712 
713   cs%diag => diag
714 
715   \textcolor{comment}{! Allocate memory for variable in control structure; note,}
716   \textcolor{comment}{! not all rows will be filled if layers get merged!}
717   \textcolor{keyword}{allocate}(cs%w\_strct(isd:ied,jsd:jed,nz+1))
718   \textcolor{keyword}{allocate}(cs%u\_strct(isd:ied,jsd:jed,nz+1))
719   \textcolor{keyword}{allocate}(cs%W\_profile(isd:ied,jsd:jed,nz+1))
720   \textcolor{keyword}{allocate}(cs%Uavg\_profile(isd:ied,jsd:jed,nz+1))
721   \textcolor{keyword}{allocate}(cs%z\_depths(isd:ied,jsd:jed,nz+1))
722   \textcolor{keyword}{allocate}(cs%N2(isd:ied,jsd:jed,nz+1))
723   \textcolor{keyword}{allocate}(cs%num\_intfaces(isd:ied,jsd:jed))
724 
725   \textcolor{comment}{! Write all relevant parameters to the model log.}
726   \textcolor{keyword}{call }log\_version(param\_file, mdl, version, \textcolor{stringliteral}{""})
727 
\end{DoxyCode}
