\hypertarget{namespacemom__wave__structure}{}\doxysection{mom\+\_\+wave\+\_\+structure Module Reference}
\label{namespacemom__wave__structure}\index{mom\_wave\_structure@{mom\_wave\_structure}}


\doxysubsection{Detailed Description}
Vertical structure functions for first baroclinic mode wave speed. \doxysubsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \mbox{\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}
\doxysubsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \mbox{\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 \mbox{\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 \mbox{\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}


\doxysubsection{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}}
\doxysubsubsection{\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{\texttt{ in}}  & {\em a} & lower diagonal with first entry equal to zero. \\
\hline
\mbox{\texttt{ in}}  & {\em b} & middle diagonal. \\
\hline
\mbox{\texttt{ in}}  & {\em c} & upper diagonal with last entry equal to zero. \\
\hline
\mbox{\texttt{ 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{\texttt{ in}}  & {\em y} & vector of known values on right hand side. \\
\hline
\mbox{\texttt{ in}}  & {\em method} & A string describing the algorithm to use \\
\hline
\mbox{\texttt{ 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}{0}
\DoxyCodeLine{563 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: a\textcolor{comment}{ !< lower diagonal with first entry equal to zero.}}
\DoxyCodeLine{564 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: b\textcolor{comment}{ !< middle diagonal.}}
\DoxyCodeLine{565 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: c\textcolor{comment}{ !< upper diagonal with last entry equal to zero.}}
\DoxyCodeLine{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}}
\DoxyCodeLine{567 \textcolor{comment}{           !! for systems of the form (e.g. average layer thickness in vertical diffusion case):}}
\DoxyCodeLine{568 \textcolor{comment}{           !! [ -\/alpha(k-\/1/2) ]                       * e(k-\/1) +}}
\DoxyCodeLine{569 \textcolor{comment}{           !! [  alpha(k-\/1/2) + alpha(k+1/2) + h(k) ] * e(k)   +}}
\DoxyCodeLine{570 \textcolor{comment}{           !! [ -\/alpha(k+1/2) ]                       * e(k+1) = y(k)}}
\DoxyCodeLine{571 \textcolor{comment}{           !! where a(k)=[-\/alpha(k-\/1/2)], b(k)=[alpha(k-\/1/2)+alpha(k+1/2) + h(k)],}}
\DoxyCodeLine{572 \textcolor{comment}{           !! and c(k)=[-\/alpha(k+1/2)]. Only used with TDMA\_H method.}}
\DoxyCodeLine{573 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}  :: y\textcolor{comment}{ !< vector of known values on right hand side.}}
\DoxyCodeLine{574   \textcolor{keywordtype}{character(len=*)},   \textcolor{keywordtype}{intent(in)}  :: method\textcolor{comment}{ !< A string describing the algorithm to use}}
\DoxyCodeLine{575 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(out)} :: x\textcolor{comment}{ !< vector of unknown values to solve for.}}
\DoxyCodeLine{576   \textcolor{comment}{! Local variables}}
\DoxyCodeLine{577   \textcolor{keywordtype}{integer} :: nrow                         \textcolor{comment}{! number of rows in A matrix}}
\DoxyCodeLine{578 \textcolor{comment}{!  real, allocatable, dimension(:,:) :: A\_check ! for solution checking}}
\DoxyCodeLine{579 \textcolor{comment}{!  real, allocatable, dimension(:)   :: y\_check ! for solution checking}}
\DoxyCodeLine{580 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:)} :: c\_prime, y\_prime, q, alpha}
\DoxyCodeLine{581                                           \textcolor{comment}{! intermediate values for solvers}}
\DoxyCodeLine{582 \textcolor{keywordtype}{  real}    :: Q\_prime, beta                \textcolor{comment}{! intermediate values for solver}}
\DoxyCodeLine{583   \textcolor{keywordtype}{integer} :: k                            \textcolor{comment}{! row (e.g. interface) index}}
\DoxyCodeLine{584   \textcolor{keywordtype}{integer} :: i,j}
\DoxyCodeLine{585 }
\DoxyCodeLine{586   nrow = \textcolor{keyword}{size}(y)}
\DoxyCodeLine{587   \textcolor{keyword}{allocate}(c\_prime(nrow))}
\DoxyCodeLine{588   \textcolor{keyword}{allocate}(y\_prime(nrow))}
\DoxyCodeLine{589   \textcolor{keyword}{allocate}(q(nrow))}
\DoxyCodeLine{590   \textcolor{keyword}{allocate}(alpha(nrow))}
\DoxyCodeLine{591 \textcolor{comment}{!  allocate(A\_check(nrow,nrow))}}
\DoxyCodeLine{592 \textcolor{comment}{!  allocate(y\_check(nrow))}}
\DoxyCodeLine{593 }
\DoxyCodeLine{594   \textcolor{keywordflow}{if} (method == \textcolor{stringliteral}{'TDMA\_T'}) \textcolor{keywordflow}{then}}
\DoxyCodeLine{595     \textcolor{comment}{! Standard Thomas algoritim (4th variant).}}
\DoxyCodeLine{596     \textcolor{comment}{! Note: Requires A to be non-\/singular for accuracy/stability}}
\DoxyCodeLine{597     c\_prime(:) = 0.0       ; y\_prime(:) = 0.0}
\DoxyCodeLine{598     c\_prime(1) = c(1)/b(1) ; y\_prime(1) = y(1)/b(1)}
\DoxyCodeLine{599 }
\DoxyCodeLine{600     \textcolor{comment}{! Forward sweep}}
\DoxyCodeLine{601     \textcolor{keywordflow}{do} k=2,nrow-\/1}
\DoxyCodeLine{602       c\_prime(k) = c(k)/(b(k)-\/a(k)*c\_prime(k-\/1))}
\DoxyCodeLine{603 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{604     \textcolor{comment}{!print *, 'c\_prime=', c\_prime(1:nrow)}}
\DoxyCodeLine{605     \textcolor{keywordflow}{do} k=2,nrow}
\DoxyCodeLine{606       y\_prime(k) = (y(k)-\/a(k)*y\_prime(k-\/1))/(b(k)-\/a(k)*c\_prime(k-\/1))}
\DoxyCodeLine{607 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{608     \textcolor{comment}{!print *, 'y\_prime=', y\_prime(1:nrow)}}
\DoxyCodeLine{609     x(nrow) = y\_prime(nrow)}
\DoxyCodeLine{610 }
\DoxyCodeLine{611     \textcolor{comment}{! Backward sweep}}
\DoxyCodeLine{612     \textcolor{keywordflow}{do} k=nrow-\/1,1,-\/1}
\DoxyCodeLine{613       x(k) = y\_prime(k)-\/c\_prime(k)*x(k+1)}
\DoxyCodeLine{614 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{615     \textcolor{comment}{!print *, 'x=',x(1:nrow)}}
\DoxyCodeLine{616 }
\DoxyCodeLine{617     \textcolor{comment}{! Check results -\/ delete later}}
\DoxyCodeLine{618     \textcolor{comment}{!do j=1,nrow ; do i=1,nrow}}
\DoxyCodeLine{619     \textcolor{comment}{!  if (i==j)then ;       A\_check(i,j) = b(i)}}
\DoxyCodeLine{620     \textcolor{comment}{!  elseif (i==j+1)then ; A\_check(i,j) = a(i)}}
\DoxyCodeLine{621     \textcolor{comment}{!  elseif (i==j-\/1)then ; A\_check(i,j) = c(i)}}
\DoxyCodeLine{622     \textcolor{comment}{!  endif}}
\DoxyCodeLine{623     \textcolor{comment}{!enddo ; enddo}}
\DoxyCodeLine{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)}}
\DoxyCodeLine{625     \textcolor{comment}{!y\_check = matmul(A\_check,x)}}
\DoxyCodeLine{626     \textcolor{comment}{!if (all(y\_check /= y))then}}
\DoxyCodeLine{627     \textcolor{comment}{!  print *, "tridiag\_solver: Uh oh, something's not right!"}}
\DoxyCodeLine{628     \textcolor{comment}{!  print *, "y=", y}}
\DoxyCodeLine{629     \textcolor{comment}{!  print *, "y\_check=", y\_check}}
\DoxyCodeLine{630     \textcolor{comment}{!endif}}
\DoxyCodeLine{631 }
\DoxyCodeLine{632   \textcolor{keywordflow}{elseif} (method == \textcolor{stringliteral}{'TDMA\_H'}) \textcolor{keywordflow}{then}}
\DoxyCodeLine{633     \textcolor{comment}{! Thomas algoritim (4th variant) w/ Hallberg substitution.}}
\DoxyCodeLine{634     \textcolor{comment}{! For a layered system where k is at interfaces, alpha\{k+1/2\} refers to}}
\DoxyCodeLine{635     \textcolor{comment}{! some property (e.g. inverse thickness for mode-\/structure problem) of the}}
\DoxyCodeLine{636     \textcolor{comment}{! layer below and alpha\{k-\/1/2\} refers to the layer above.}}
\DoxyCodeLine{637     \textcolor{comment}{! Here, alpha(k)=alpha\{k+1/2\} and alpha(k-\/1)=alpha\{k-\/1/2\}.}}
\DoxyCodeLine{638     \textcolor{comment}{! Strictly speaking, this formulation requires A to be a non-\/singular,}}
\DoxyCodeLine{639     \textcolor{comment}{! symmetric, diagonally dominant matrix, with h>0.}}
\DoxyCodeLine{640     \textcolor{comment}{! Need to add a check for these conditions.}}
\DoxyCodeLine{641     \textcolor{keywordflow}{do} k=1,nrow-\/1}
\DoxyCodeLine{642       \textcolor{keywordflow}{if} (abs(a(k+1)-\/c(k)) > 1.e-\/10*(abs(a(k+1))+abs(c(k)))) \textcolor{keywordflow}{then}}
\DoxyCodeLine{643         \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"tridiag\_solver: matrix not symmetric; need symmetry when invoking TDMA\_H"})}
\DoxyCodeLine{644 \textcolor{keywordflow}{      endif}}
\DoxyCodeLine{645 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{646     alpha = -\/c}
\DoxyCodeLine{647     \textcolor{comment}{! Alpha of the bottom-\/most layer is not necessarily zero. Therefore,}}
\DoxyCodeLine{648     \textcolor{comment}{! back out the value from the provided b(nrow and h(nrow) values}}
\DoxyCodeLine{649     alpha(nrow)   = b(nrow)-\/h(nrow)-\/alpha(nrow-\/1)}
\DoxyCodeLine{650     \textcolor{comment}{! Prime other variables}}
\DoxyCodeLine{651     beta       = 1/b(1)}
\DoxyCodeLine{652     y\_prime(:) = 0.0       ; q(:) = 0.0}
\DoxyCodeLine{653     y\_prime(1) = beta*y(1) ; q(1) = beta*alpha(1)}
\DoxyCodeLine{654     q\_prime    = 1-\/q(1)}
\DoxyCodeLine{655 }
\DoxyCodeLine{656     \textcolor{comment}{! Forward sweep}}
\DoxyCodeLine{657     \textcolor{keywordflow}{do} k=2,nrow-\/1}
\DoxyCodeLine{658       beta = 1/(h(k)+alpha(k-\/1)*q\_prime+alpha(k))}
\DoxyCodeLine{659       \textcolor{keywordflow}{if} (isnan(beta))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Tridiag\_solver: beta is NAN"} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{660       q(k) = beta*alpha(k)}
\DoxyCodeLine{661       y\_prime(k) = beta*(y(k)+alpha(k-\/1)*y\_prime(k-\/1))}
\DoxyCodeLine{662       q\_prime = beta*(h(k)+alpha(k-\/1)*q\_prime)}
\DoxyCodeLine{663 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{664     \textcolor{keywordflow}{if} ((h(nrow)+alpha(nrow-\/1)*q\_prime+alpha(nrow)) == 0.0)\textcolor{keywordflow}{then}}
\DoxyCodeLine{665       \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"Tridiag\_solver: this system is not stable."}) \textcolor{comment}{! ; overriding beta(nrow)}}
\DoxyCodeLine{666       \textcolor{comment}{! This has hard-\/coded dimensions: beta = 1/(1e-\/15) ! place holder for unstable systems -\/ delete later}}
\DoxyCodeLine{667     \textcolor{keywordflow}{else}}
\DoxyCodeLine{668       beta = 1/(h(nrow)+alpha(nrow-\/1)*q\_prime+alpha(nrow))}
\DoxyCodeLine{669 \textcolor{keywordflow}{    endif}}
\DoxyCodeLine{670     y\_prime(nrow) = beta*(y(nrow)+alpha(nrow-\/1)*y\_prime(nrow-\/1))}
\DoxyCodeLine{671     x(nrow) = y\_prime(nrow)}
\DoxyCodeLine{672     \textcolor{comment}{! Backward sweep}}
\DoxyCodeLine{673     \textcolor{keywordflow}{do} k=nrow-\/1,1,-\/1}
\DoxyCodeLine{674       x(k) = y\_prime(k)+q(k)*x(k+1)}
\DoxyCodeLine{675 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{676     \textcolor{comment}{!print *, 'yprime=',y\_prime(1:nrow)}}
\DoxyCodeLine{677     \textcolor{comment}{!print *, 'x=',x(1:nrow)}}
\DoxyCodeLine{678 \textcolor{keywordflow}{  endif}}
\DoxyCodeLine{679 }
\DoxyCodeLine{680   \textcolor{keyword}{deallocate}(c\_prime,y\_prime,q,alpha)}
\DoxyCodeLine{681 \textcolor{comment}{! deallocate(A\_check,y\_check)}}
\DoxyCodeLine{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}}
\doxysubsubsection{\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(\mbox{\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{\texttt{ in}}  & {\em g} & The ocean\textquotesingle{}s grid structure. \\
\hline
\mbox{\texttt{ in}}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure. \\
\hline
\mbox{\texttt{ in}}  & {\em us} & A dimensional unit scaling type \\
\hline
\mbox{\texttt{ in}}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]} \\
\hline
\mbox{\texttt{ in}}  & {\em tv} & A structure pointing to various thermodynamic variables. \\
\hline
\mbox{\texttt{ in}}  & {\em cn} & The (non-\/rotational) mode internal gravity wave speed \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}. \\
\hline
\mbox{\texttt{ in}}  & {\em modenum} & Mode number \\
\hline
\mbox{\texttt{ 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{\texttt{ in}}  & {\em en} & Internal wave energy density \mbox{[}R Z3 T-\/2 $\sim$$>$ J m-\/2\mbox{]} \\
\hline
\mbox{\texttt{ 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}{0}
\DoxyCodeLine{92   \textcolor{keywordtype}{type}(ocean\_grid\_type),                    \textcolor{keywordtype}{intent(in)}  :: G\textcolor{comment}{  !< The ocean's grid structure.}}
\DoxyCodeLine{93   \textcolor{keywordtype}{type}(verticalGrid\_type),                  \textcolor{keywordtype}{intent(in)}  :: GV\textcolor{comment}{ !< The ocean's vertical grid structure.}}
\DoxyCodeLine{94   \textcolor{keywordtype}{type}(unit\_scale\_type),                    \textcolor{keywordtype}{intent(in)}  :: US\textcolor{comment}{ !< A dimensional unit scaling type}}
\DoxyCodeLine{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]}}
\DoxyCodeLine{96   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),                    \textcolor{keywordtype}{intent(in)}  :: tv\textcolor{comment}{ !< A structure pointing to various}}
\DoxyCodeLine{97 \textcolor{comment}{                                                              !! thermodynamic variables.}}
\DoxyCodeLine{98 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))},         \textcolor{keywordtype}{intent(in)}  :: cn\textcolor{comment}{ !< The (non-\/rotational) mode internal}}
\DoxyCodeLine{99 \textcolor{comment}{                                                              !! gravity wave speed [L T-\/1 ~> m s-\/1].}}
\DoxyCodeLine{100   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}  :: ModeNum\textcolor{comment}{ !< Mode number}}
\DoxyCodeLine{101 \textcolor{keywordtype}{  real},                                     \textcolor{keywordtype}{intent(in)}  :: freq\textcolor{comment}{ !< Intrinsic wave frequency [T-\/1 ~> s-\/1].}}
\DoxyCodeLine{102   \textcolor{keywordtype}{type}(wave\_structure\_CS),                  \textcolor{keywordtype}{pointer}     :: CS\textcolor{comment}{ !< The control structure returned by a}}
\DoxyCodeLine{103 \textcolor{comment}{                                                              !! previous call to wave\_structure\_init.}}
\DoxyCodeLine{104 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, \&}
\DoxyCodeLine{105                                   \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: En\textcolor{comment}{ !< Internal wave energy density [R Z3 T-\/2 ~> J m-\/2]}}
\DoxyCodeLine{106   \textcolor{keywordtype}{logical},                        \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: full\_halos\textcolor{comment}{ !< If true, do the calculation}}
\DoxyCodeLine{107 \textcolor{comment}{                                                              !! over the entire computational domain.}}
\DoxyCodeLine{108   \textcolor{comment}{! Local variables}}
\DoxyCodeLine{109 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: \&}
\DoxyCodeLine{110     dRho\_dT, \&    \textcolor{comment}{! Partial derivative of density with temperature [R degC-\/1 ~> kg m-\/3 degC-\/1]}}
\DoxyCodeLine{111     dRho\_dS, \&    \textcolor{comment}{! Partial derivative of density with salinity [R ppt-\/1 ~> kg m-\/3 ppt-\/1]}}
\DoxyCodeLine{112     pres, \&       \textcolor{comment}{! Interface pressure [R L2 T-\/2 ~> Pa]}}
\DoxyCodeLine{113     T\_int, \&      \textcolor{comment}{! Temperature interpolated to interfaces [degC]}}
\DoxyCodeLine{114     S\_int, \&      \textcolor{comment}{! Salinity interpolated to interfaces [ppt]}}
\DoxyCodeLine{115     gprime        \textcolor{comment}{! The reduced gravity across each interface [m2 Z-\/1 s-\/2 ~> m s-\/2].}}
\DoxyCodeLine{116 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: \&}
\DoxyCodeLine{117     Igl, Igu      \textcolor{comment}{! The inverse of the reduced gravity across an interface times}}
\DoxyCodeLine{118                   \textcolor{comment}{! the thickness of the layer below (Igl) or above (Igu) it [s2 m-\/2].}}
\DoxyCodeLine{119 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G),SZI\_(G))} :: \&}
\DoxyCodeLine{120     Hf, \&         \textcolor{comment}{! Layer thicknesses after very thin layers are combined [Z ~> m]}}
\DoxyCodeLine{121     Tf, \&         \textcolor{comment}{! Layer temperatures after very thin layers are combined [degC]}}
\DoxyCodeLine{122     Sf, \&         \textcolor{comment}{! Layer salinities after very thin layers are combined [ppt]}}
\DoxyCodeLine{123     Rf            \textcolor{comment}{! Layer densities after very thin layers are combined [R ~> kg m-\/3]}}
\DoxyCodeLine{124 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: \&}
\DoxyCodeLine{125     Hc, \&         \textcolor{comment}{! A column of layer thicknesses after convective istabilities are removed [Z ~> m]}}
\DoxyCodeLine{126     Tc, \&         \textcolor{comment}{! A column of layer temperatures after convective istabilities are removed [degC]}}
\DoxyCodeLine{127     Sc, \&         \textcolor{comment}{! A column of layer salinites after convective istabilities are removed [ppt]}}
\DoxyCodeLine{128     Rc, \&         \textcolor{comment}{! A column of layer densities after convective istabilities are removed [R ~> kg m-\/3]}}
\DoxyCodeLine{129     det, ddet}
\DoxyCodeLine{130 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))} :: \&}
\DoxyCodeLine{131     htot          \textcolor{comment}{! The vertical sum of the thicknesses [Z ~> m]}}
\DoxyCodeLine{132 \textcolor{keywordtype}{  real} :: lam}
\DoxyCodeLine{133 \textcolor{keywordtype}{  real} :: min\_h\_frac}
\DoxyCodeLine{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]}}
\DoxyCodeLine{135 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: \&}
\DoxyCodeLine{136     hmin, \&        \textcolor{comment}{! Thicknesses [Z ~> m]}}
\DoxyCodeLine{137     H\_here, \&      \textcolor{comment}{! A thickness [Z ~> m]}}
\DoxyCodeLine{138     HxT\_here, \&    \textcolor{comment}{! A layer integrated temperature [degC Z ~> degC m]}}
\DoxyCodeLine{139     HxS\_here, \&    \textcolor{comment}{! A layer integrated salinity [ppt Z ~> ppt m]}}
\DoxyCodeLine{140     HxR\_here       \textcolor{comment}{! A layer integrated density [R Z ~> kg m-\/2]}}
\DoxyCodeLine{141 \textcolor{keywordtype}{  real} :: speed2\_tot}
\DoxyCodeLine{142 \textcolor{keywordtype}{  real} :: I\_Hnew   \textcolor{comment}{! The inverse of a new layer thickness [Z-\/1 ~> m-\/1]}}
\DoxyCodeLine{143 \textcolor{keywordtype}{  real} :: drxh\_sum \textcolor{comment}{! The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-\/2]}}
\DoxyCodeLine{144 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{parameter} :: tol1  = 0.0001, tol2 = 0.001}
\DoxyCodeLine{145 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{pointer}, \textcolor{keywordtype}{dimension(:,:,:)} :: T => null(), s => null()}
\DoxyCodeLine{146 \textcolor{keywordtype}{  real} :: g\_Rho0  \textcolor{comment}{! G\_Earth/Rho0 in [m2 s-\/2 Z-\/1 R-\/1 ~> m4 s-\/2 kg-\/1].}}
\DoxyCodeLine{147   \textcolor{comment}{! real :: rescale, I\_rescale}}
\DoxyCodeLine{148   \textcolor{keywordtype}{integer} :: kf(SZI\_(G))}
\DoxyCodeLine{149   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: max\_itt = 1 \textcolor{comment}{! number of times to iterate in solving for eigenvector}}
\DoxyCodeLine{150 \textcolor{keywordtype}{  real} :: cg\_subRO        \textcolor{comment}{! A tiny wave speed to prevent division by zero [L T-\/1 ~> m s-\/1]}}
\DoxyCodeLine{151 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{parameter}    :: a\_int = 0.5 \textcolor{comment}{! value of normalized integral: \(\backslash\)int(w\_strct\string^2)dz = a\_int}}
\DoxyCodeLine{152 \textcolor{keywordtype}{  real}               :: I\_a\_int     \textcolor{comment}{! inverse of a\_int}}
\DoxyCodeLine{153 \textcolor{keywordtype}{  real}               :: f2          \textcolor{comment}{! squared Coriolis frequency [T-\/2 ~> s-\/2]}}
\DoxyCodeLine{154 \textcolor{keywordtype}{  real}               :: Kmag2       \textcolor{comment}{! magnitude of horizontal wave number squared}}
\DoxyCodeLine{155   \textcolor{keywordtype}{logical}            :: use\_EOS     \textcolor{comment}{! If true, density is calculated from T \& S using an}}
\DoxyCodeLine{156                                     \textcolor{comment}{! equation of state.}}
\DoxyCodeLine{157 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: w\_strct, u\_strct, W\_profile, Uavg\_profile, z\_int, N2}
\DoxyCodeLine{158                                         \textcolor{comment}{! local representations of variables in CS; note,}}
\DoxyCodeLine{159                                         \textcolor{comment}{! not all rows will be filled if layers get merged!}}
\DoxyCodeLine{160 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: w\_strct2, u\_strct2}
\DoxyCodeLine{161                                         \textcolor{comment}{! squared values}}
\DoxyCodeLine{162 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G))}   :: dz      \textcolor{comment}{! thicknesses of merged layers (same as Hc I hope)}}
\DoxyCodeLine{163   \textcolor{comment}{! real, dimension(SZK\_(G)+1) :: dWdz\_profile ! profile of dW/dz}}
\DoxyCodeLine{164 \textcolor{keywordtype}{  real}                       :: w2avg   \textcolor{comment}{! average of squared vertical velocity structure funtion}}
\DoxyCodeLine{165 \textcolor{keywordtype}{  real}                       :: int\_dwdz2}
\DoxyCodeLine{166 \textcolor{keywordtype}{  real}                       :: int\_w2}
\DoxyCodeLine{167 \textcolor{keywordtype}{  real}                       :: int\_N2w2}
\DoxyCodeLine{168 \textcolor{keywordtype}{  real}                       :: KE\_term \textcolor{comment}{! terms in vertically averaged energy equation}}
\DoxyCodeLine{169 \textcolor{keywordtype}{  real}                       :: PE\_term \textcolor{comment}{! terms in vertically averaged energy equation}}
\DoxyCodeLine{170 \textcolor{keywordtype}{  real}                       :: W0      \textcolor{comment}{! A vertical velocity magnitude [Z T-\/1 ~> m s-\/1]}}
\DoxyCodeLine{171 \textcolor{keywordtype}{  real}                       :: gp\_unscaled \textcolor{comment}{! A version of gprime rescaled to [m s-\/2].}}
\DoxyCodeLine{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}}
\DoxyCodeLine{173                                         \textcolor{comment}{! interface (excluding surface and bottom)}}
\DoxyCodeLine{174 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)-\/1)} :: a\_diag, b\_diag, c\_diag}
\DoxyCodeLine{175                                         \textcolor{comment}{! diagonals of tridiagonal matrix; one value for each}}
\DoxyCodeLine{176                                         \textcolor{comment}{! interface (excluding surface and bottom)}}
\DoxyCodeLine{177 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)-\/1)} :: e\_guess \textcolor{comment}{! guess at eigen vector with unit amplitde (for TDMA)}}
\DoxyCodeLine{178 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)-\/1)} :: e\_itt   \textcolor{comment}{! improved guess at eigen vector (from TDMA)}}
\DoxyCodeLine{179 \textcolor{keywordtype}{  real}    :: Pi}
\DoxyCodeLine{180   \textcolor{keywordtype}{integer} :: kc}
\DoxyCodeLine{181   \textcolor{keywordtype}{integer} :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig\_stop, jg\_stop}
\DoxyCodeLine{182 }
\DoxyCodeLine{183   is = g\%isc ; ie = g\%iec ; js = g\%jsc ; je = g\%jec ; nz = g\%ke}
\DoxyCodeLine{184   i\_a\_int = 1/a\_int}
\DoxyCodeLine{185 }
\DoxyCodeLine{186   \textcolor{comment}{!if (present(CS)) then}}
\DoxyCodeLine{187     \textcolor{keywordflow}{if} (.not. \textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_wave\_structure: "}// \&}
\DoxyCodeLine{188            \textcolor{stringliteral}{"Module must be initialized before it is used."})}
\DoxyCodeLine{189   \textcolor{comment}{!endif}}
\DoxyCodeLine{190 }
\DoxyCodeLine{191   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(full\_halos)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (full\_halos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{192     is = g\%isd ; ie = g\%ied ; js = g\%jsd ; je = g\%jed}
\DoxyCodeLine{193 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{194 }
\DoxyCodeLine{195   pi = (4.0*atan(1.0))}
\DoxyCodeLine{196 }
\DoxyCodeLine{197   s => tv\%S ; t => tv\%T}
\DoxyCodeLine{198   g\_rho0 = us\%L\_T\_to\_m\_s**2 * gv\%g\_Earth / gv\%Rho0}
\DoxyCodeLine{199   cg\_subro = 1e-\/100*us\%m\_s\_to\_L\_T  \textcolor{comment}{! The hard-\/coded value here might need to increase.}}
\DoxyCodeLine{200   use\_eos = \textcolor{keyword}{associated}(tv\%eqn\_of\_state)}
\DoxyCodeLine{201 }
\DoxyCodeLine{202   \textcolor{comment}{! Simplifying the following could change answers at roundoff.}}
\DoxyCodeLine{203   z\_to\_pres = gv\%Z\_to\_H * (gv\%H\_to\_RZ * gv\%g\_Earth)}
\DoxyCodeLine{204   \textcolor{comment}{! rescale = 1024.0**4 ; I\_rescale = 1.0/rescale}}
\DoxyCodeLine{205 }
\DoxyCodeLine{206   min\_h\_frac = tol1 / real(nz)}
\DoxyCodeLine{207 }
\DoxyCodeLine{208   \textcolor{keywordflow}{do} j=js,je}
\DoxyCodeLine{209     \textcolor{comment}{!   First merge very thin layers with the one above (or below if they are}}
\DoxyCodeLine{210     \textcolor{comment}{! at the top).  This also transposes the row order so that columns can}}
\DoxyCodeLine{211     \textcolor{comment}{! be worked upon one at a time.}}
\DoxyCodeLine{212     \textcolor{keywordflow}{do} i=is,ie ; htot(i,j) = 0.0 ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{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}}
\DoxyCodeLine{214 }
\DoxyCodeLine{215     \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{216       hmin(i) = htot(i,j)*min\_h\_frac ; kf(i) = 1 ; h\_here(i) = 0.0}
\DoxyCodeLine{217       hxt\_here(i) = 0.0 ; hxs\_here(i) = 0.0 ; hxr\_here(i) = 0.0}
\DoxyCodeLine{218 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{219     \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{220       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{221         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv\%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}}
\DoxyCodeLine{222           hf(kf(i),i) = h\_here(i)}
\DoxyCodeLine{223           tf(kf(i),i) = hxt\_here(i) / h\_here(i)}
\DoxyCodeLine{224           sf(kf(i),i) = hxs\_here(i) / h\_here(i)}
\DoxyCodeLine{225           kf(i) = kf(i) + 1}
\DoxyCodeLine{226 }
\DoxyCodeLine{227           \textcolor{comment}{! Start a new layer}}
\DoxyCodeLine{228           h\_here(i) = h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{229           hxt\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*t(i,j,k)}
\DoxyCodeLine{230           hxs\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*s(i,j,k)}
\DoxyCodeLine{231         \textcolor{keywordflow}{else}}
\DoxyCodeLine{232           h\_here(i) = h\_here(i) + h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{233           hxt\_here(i) = hxt\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*t(i,j,k)}
\DoxyCodeLine{234           hxs\_here(i) = hxs\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*s(i,j,k)}
\DoxyCodeLine{235 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{236 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{237       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{238         hf(kf(i),i) = h\_here(i)}
\DoxyCodeLine{239         tf(kf(i),i) = hxt\_here(i) / h\_here(i)}
\DoxyCodeLine{240         sf(kf(i),i) = hxs\_here(i) / h\_here(i)}
\DoxyCodeLine{241 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{242     \textcolor{keywordflow}{else}}
\DoxyCodeLine{243       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{244         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv\%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}}
\DoxyCodeLine{245           hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)}
\DoxyCodeLine{246           kf(i) = kf(i) + 1}
\DoxyCodeLine{247 }
\DoxyCodeLine{248           \textcolor{comment}{! Start a new layer}}
\DoxyCodeLine{249           h\_here(i) = h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{250           hxr\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*gv\%Rlay(k)}
\DoxyCodeLine{251         \textcolor{keywordflow}{else}}
\DoxyCodeLine{252           h\_here(i) = h\_here(i) + h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{253           hxr\_here(i) = hxr\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*gv\%Rlay(k)}
\DoxyCodeLine{254 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{255 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{256       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{257         hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)}
\DoxyCodeLine{258 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{259 \textcolor{keywordflow}{    endif} \textcolor{comment}{! use\_EOS?}}
\DoxyCodeLine{260 }
\DoxyCodeLine{261     \textcolor{comment}{! From this point, we can work on individual columns without causing memory}}
\DoxyCodeLine{262     \textcolor{comment}{! to have page faults.}}
\DoxyCodeLine{263     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (cn(i,j)>0.0)\textcolor{keywordflow}{then}}
\DoxyCodeLine{264       \textcolor{comment}{!-\/-\/-\/-\/for debugging, remove later-\/-\/-\/-\/}}
\DoxyCodeLine{265       ig = i + g\%idg\_offset ; jg = j + g\%jdg\_offset}
\DoxyCodeLine{266       \textcolor{comment}{!if (ig == CS\%int\_tide\_source\_x .and. jg == CS\%int\_tide\_source\_y) then}}
\DoxyCodeLine{267       \textcolor{comment}{!-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/}}
\DoxyCodeLine{268       \textcolor{keywordflow}{if} (g\%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then}}
\DoxyCodeLine{269 }
\DoxyCodeLine{270         lam = 1/(us\%L\_T\_to\_m\_s**2 * cn(i,j)**2)}
\DoxyCodeLine{271 }
\DoxyCodeLine{272         \textcolor{comment}{! Calculate drxh\_sum}}
\DoxyCodeLine{273         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{274           pres(1) = 0.0}
\DoxyCodeLine{275           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{276             pres(k) = pres(k-\/1) + z\_to\_pres*hf(k-\/1,i)}
\DoxyCodeLine{277             t\_int(k) = 0.5*(tf(k,i)+tf(k-\/1,i))}
\DoxyCodeLine{278             s\_int(k) = 0.5*(sf(k,i)+sf(k-\/1,i))}
\DoxyCodeLine{279 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{280           \textcolor{keyword}{call }calculate\_density\_derivs(t\_int, s\_int, pres, drho\_dt, drho\_ds, \&}
\DoxyCodeLine{281                                         tv\%eqn\_of\_state, (/2,kf(i)/) )}
\DoxyCodeLine{282 }
\DoxyCodeLine{283           \textcolor{comment}{! Sum the reduced gravities to find out how small a density difference}}
\DoxyCodeLine{284           \textcolor{comment}{! is negligibly small.}}
\DoxyCodeLine{285           drxh\_sum = 0.0}
\DoxyCodeLine{286           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{287             drxh\_sum = drxh\_sum + 0.5*(hf(k-\/1,i)+hf(k,i)) * \&}
\DoxyCodeLine{288                 max(0.0,drho\_dt(k)*(tf(k,i)-\/tf(k-\/1,i)) + \&}
\DoxyCodeLine{289                         drho\_ds(k)*(sf(k,i)-\/sf(k-\/1,i)))}
\DoxyCodeLine{290 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{291         \textcolor{keywordflow}{else}}
\DoxyCodeLine{292           drxh\_sum = 0.0}
\DoxyCodeLine{293           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{294             drxh\_sum = drxh\_sum + 0.5*(hf(k-\/1,i)+hf(k,i)) * \&}
\DoxyCodeLine{295                               max(0.0,rf(k,i)-\/rf(k-\/1,i))}
\DoxyCodeLine{296 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{297 \textcolor{keywordflow}{        endif} \textcolor{comment}{! use\_EOS?}}
\DoxyCodeLine{298 }
\DoxyCodeLine{299         \textcolor{comment}{!   Find gprime across each internal interface, taking care of convective}}
\DoxyCodeLine{300         \textcolor{comment}{! instabilities by merging layers.}}
\DoxyCodeLine{301         \textcolor{keywordflow}{if} (drxh\_sum >= 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{302           \textcolor{comment}{! Merge layers to eliminate convective instabilities or exceedingly}}
\DoxyCodeLine{303           \textcolor{comment}{! small reduced gravities.}}
\DoxyCodeLine{304           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{305             kc = 1}
\DoxyCodeLine{306             hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)}
\DoxyCodeLine{307             \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{308               \textcolor{keywordflow}{if} ((drho\_dt(k)*(tf(k,i)-\/tc(kc)) + drho\_ds(k)*(sf(k,i)-\/sc(kc))) * \&}
\DoxyCodeLine{309                   (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh\_sum) \textcolor{keywordflow}{then}}
\DoxyCodeLine{310                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}}
\DoxyCodeLine{311                 i\_hnew = 1.0 / (hc(kc) + hf(k,i))}
\DoxyCodeLine{312                 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i\_hnew}
\DoxyCodeLine{313                 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i\_hnew}
\DoxyCodeLine{314                 hc(kc) = (hc(kc) + hf(k,i))}
\DoxyCodeLine{315                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}}
\DoxyCodeLine{316                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}}
\DoxyCodeLine{317                 \textcolor{comment}{! far back we go.}}
\DoxyCodeLine{318                 \textcolor{keywordflow}{do} k2=kc,2,-\/1}
\DoxyCodeLine{319                   \textcolor{keywordflow}{if} ((drho\_dt(k2)*(tc(k2)-\/tc(k2-\/1)) + drho\_ds(k2)*(sc(k2)-\/sc(k2-\/1))) * \&}
\DoxyCodeLine{320                       (hc(k2) + hc(k2-\/1)) < tol2*drxh\_sum) \textcolor{keywordflow}{then}}
\DoxyCodeLine{321                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}}
\DoxyCodeLine{322                     i\_hnew = 1.0 / (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{323                     tc(kc-\/1) = (hc(kc)*tc(kc) + hc(kc-\/1)*tc(kc-\/1)) * i\_hnew}
\DoxyCodeLine{324                     sc(kc-\/1) = (hc(kc)*sc(kc) + hc(kc-\/1)*sc(kc-\/1)) * i\_hnew}
\DoxyCodeLine{325                     hc(kc-\/1) = (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{326                     kc = kc -\/ 1}
\DoxyCodeLine{327                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{328 \textcolor{keywordflow}{                enddo}}
\DoxyCodeLine{329               \textcolor{keywordflow}{else}}
\DoxyCodeLine{330                 \textcolor{comment}{! Add a new layer to the column.}}
\DoxyCodeLine{331                 kc = kc + 1}
\DoxyCodeLine{332                 drho\_ds(kc) = drho\_ds(k) ; drho\_dt(kc) = drho\_dt(k)}
\DoxyCodeLine{333                 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)}
\DoxyCodeLine{334 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{335 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{336             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}}
\DoxyCodeLine{337             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-\/Boussinesq.}}
\DoxyCodeLine{338               gprime(k) = g\_rho0 * (drho\_dt(k)*(tc(k)-\/tc(k-\/1)) + \&}
\DoxyCodeLine{339                                     drho\_ds(k)*(sc(k)-\/sc(k-\/1)))}
\DoxyCodeLine{340 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{341           \textcolor{keywordflow}{else}  \textcolor{comment}{! .not.use\_EOS}}
\DoxyCodeLine{342             \textcolor{comment}{! Do the same with density directly...}}
\DoxyCodeLine{343             kc = 1}
\DoxyCodeLine{344             hc(1) = hf(1,i) ; rc(1) = rf(1,i)}
\DoxyCodeLine{345             \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{346               \textcolor{keywordflow}{if} ((rf(k,i) -\/ rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh\_sum) \textcolor{keywordflow}{then}}
\DoxyCodeLine{347                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}}
\DoxyCodeLine{348                 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))}
\DoxyCodeLine{349                 hc(kc) = (hc(kc) + hf(k,i))}
\DoxyCodeLine{350                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}}
\DoxyCodeLine{351                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}}
\DoxyCodeLine{352                 \textcolor{comment}{! far back we go.}}
\DoxyCodeLine{353                 \textcolor{keywordflow}{do} k2=kc,2,-\/1}
\DoxyCodeLine{354                   \textcolor{keywordflow}{if} ((rc(k2)-\/rc(k2-\/1)) * (hc(k2)+hc(k2-\/1)) < tol2*drxh\_sum) \textcolor{keywordflow}{then}}
\DoxyCodeLine{355                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}}
\DoxyCodeLine{356                     rc(kc-\/1) = (hc(kc)*rc(kc) + hc(kc-\/1)*rc(kc-\/1)) / (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{357                     hc(kc-\/1) = (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{358                     kc = kc -\/ 1}
\DoxyCodeLine{359                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{360 \textcolor{keywordflow}{                enddo}}
\DoxyCodeLine{361               \textcolor{keywordflow}{else}}
\DoxyCodeLine{362                 \textcolor{comment}{! Add a new layer to the column.}}
\DoxyCodeLine{363                 kc = kc + 1}
\DoxyCodeLine{364                 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)}
\DoxyCodeLine{365 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{366 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{367             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}}
\DoxyCodeLine{368             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-\/Boussinesq.}}
\DoxyCodeLine{369               gprime(k) = g\_rho0 * (rc(k)-\/rc(k-\/1))}
\DoxyCodeLine{370 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{371 \textcolor{keywordflow}{          endif}  \textcolor{comment}{! use\_EOS?}}
\DoxyCodeLine{372 }
\DoxyCodeLine{373           \textcolor{comment}{!-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/NOW FIND WAVE STRUCTURE-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/}}
\DoxyCodeLine{374           \textcolor{comment}{! Construct and solve tridiagonal system for the interior interfaces}}
\DoxyCodeLine{375           \textcolor{comment}{! Note that kc   = number of layers,}}
\DoxyCodeLine{376           \textcolor{comment}{!           kc+1 = nzm = number of interfaces,}}
\DoxyCodeLine{377           \textcolor{comment}{!           kc-\/1 = number of interior interfaces (excluding surface and bottom)}}
\DoxyCodeLine{378           \textcolor{comment}{! Also, note that "K" refers to an interface, while "k" refers to the layer below.}}
\DoxyCodeLine{379           \textcolor{comment}{! Need at least 3 layers (2 internal interfaces) to generate a matrix, also}}
\DoxyCodeLine{380           \textcolor{comment}{! need number of layers to be greater than the mode number}}
\DoxyCodeLine{381           \textcolor{keywordflow}{if} (kc >= max(3, modenum + 1)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{382             \textcolor{comment}{! Set depth at surface}}
\DoxyCodeLine{383             z\_int(1) = 0.0}
\DoxyCodeLine{384             \textcolor{comment}{! Calculate Igu, Igl, depth, and N2 at each interior interface}}
\DoxyCodeLine{385             \textcolor{comment}{! [excludes surface (K=1) and bottom (K=kc+1)]}}
\DoxyCodeLine{386             \textcolor{keywordflow}{do} k=2,kc}
\DoxyCodeLine{387               igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-\/1))}
\DoxyCodeLine{388               z\_int(k) = z\_int(k-\/1) + hc(k-\/1)}
\DoxyCodeLine{389               n2(k) = us\%m\_to\_Z**2*gprime(k)/(0.5*(hc(k)+hc(k-\/1)))}
\DoxyCodeLine{390 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{391             \textcolor{comment}{! Set stratification for surface and bottom (setting equal to nearest interface for now)}}
\DoxyCodeLine{392             n2(1) = n2(2) ; n2(kc+1) = n2(kc)}
\DoxyCodeLine{393             \textcolor{comment}{! Calcualte depth at bottom}}
\DoxyCodeLine{394             z\_int(kc+1) = z\_int(kc)+hc(kc)}
\DoxyCodeLine{395             \textcolor{comment}{! check that thicknesses sum to total depth}}
\DoxyCodeLine{396             \textcolor{keywordflow}{if} (abs(z\_int(kc+1)-\/htot(i,j)) > 1.e-\/14*htot(i,j)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{397               \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"wave\_structure: mismatch in total depths"})}
\DoxyCodeLine{398 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{399 }
\DoxyCodeLine{400             \textcolor{comment}{! Note that many of the calcluation from here on revert to using vertical}}
\DoxyCodeLine{401             \textcolor{comment}{! distances in m, not Z.}}
\DoxyCodeLine{402 }
\DoxyCodeLine{403             \textcolor{comment}{! Populate interior rows of tridiagonal matrix; must multiply through by}}
\DoxyCodeLine{404             \textcolor{comment}{! gprime to get tridiagonal matrix to the symmetrical form:}}
\DoxyCodeLine{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,}}
\DoxyCodeLine{406             \textcolor{comment}{! where lam\_z = lam*gprime is now a function of depth.}}
\DoxyCodeLine{407             \textcolor{comment}{! Frist, populate interior rows}}
\DoxyCodeLine{408             \textcolor{keywordflow}{do} k=3,kc-\/1}
\DoxyCodeLine{409               row = k-\/1 \textcolor{comment}{! indexing for TD matrix rows}}
\DoxyCodeLine{410               gp\_unscaled = us\%m\_to\_Z*gprime(k)}
\DoxyCodeLine{411               lam\_z(row) = lam*gp\_unscaled}
\DoxyCodeLine{412               a\_diag(row) = gp\_unscaled*(-\/igu(k))}
\DoxyCodeLine{413               b\_diag(row) = gp\_unscaled*(igu(k)+igl(k)) -\/ lam\_z(row)}
\DoxyCodeLine{414               c\_diag(row) = gp\_unscaled*(-\/igl(k))}
\DoxyCodeLine{415               \textcolor{keywordflow}{if} (isnan(lam\_z(row)))\textcolor{keywordflow}{then}  ; print *, \textcolor{stringliteral}{"Wave\_structure: lam\_z(row) is NAN"} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{416               \textcolor{keywordflow}{if} (isnan(a\_diag(row)))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Wave\_structure: a(k) is NAN"} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{417               \textcolor{keywordflow}{if} (isnan(b\_diag(row)))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Wave\_structure: b(k) is NAN"} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{418               \textcolor{keywordflow}{if} (isnan(c\_diag(row)))\textcolor{keywordflow}{then} ; print *, \textcolor{stringliteral}{"Wave\_structure: c(k) is NAN"} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{419 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{420             \textcolor{comment}{! Populate top row of tridiagonal matrix}}
\DoxyCodeLine{421             k=2 ; row = k-\/1 ;}
\DoxyCodeLine{422             gp\_unscaled = us\%m\_to\_Z*gprime(k)}
\DoxyCodeLine{423             lam\_z(row) = lam*gp\_unscaled}
\DoxyCodeLine{424             a\_diag(row) = 0.0}
\DoxyCodeLine{425             b\_diag(row) = gp\_unscaled*(igu(k)+igl(k)) -\/ lam\_z(row)}
\DoxyCodeLine{426             c\_diag(row) = gp\_unscaled*(-\/igl(k))}
\DoxyCodeLine{427             \textcolor{comment}{! Populate bottom row of tridiagonal matrix}}
\DoxyCodeLine{428             k=kc ; row = k-\/1}
\DoxyCodeLine{429             gp\_unscaled = us\%m\_to\_Z*gprime(k)}
\DoxyCodeLine{430             lam\_z(row) = lam*gp\_unscaled}
\DoxyCodeLine{431             a\_diag(row) = gp\_unscaled*(-\/igu(k))}
\DoxyCodeLine{432             b\_diag(row) = gp\_unscaled*(igu(k)+igl(k)) -\/ lam\_z(row)}
\DoxyCodeLine{433             c\_diag(row) = 0.0}
\DoxyCodeLine{434 }
\DoxyCodeLine{435             \textcolor{comment}{! Guess a vector shape to start with (excludes surface and bottom)}}
\DoxyCodeLine{436             e\_guess(1:kc-\/1) = sin((z\_int(2:kc)/htot(i,j)) *pi)}
\DoxyCodeLine{437             e\_guess(1:kc-\/1) = e\_guess(1:kc-\/1)/sqrt(sum(e\_guess(1:kc-\/1)**2))}
\DoxyCodeLine{438 }
\DoxyCodeLine{439             \textcolor{comment}{! Perform inverse iteration with tri-\/diag solver}}
\DoxyCodeLine{440             \textcolor{keywordflow}{do} itt=1,max\_itt}
\DoxyCodeLine{441               \textcolor{keyword}{call }tridiag\_solver(a\_diag(1:kc-\/1),b\_diag(1:kc-\/1),c\_diag(1:kc-\/1), \&}
\DoxyCodeLine{442                                   -\/lam\_z(1:kc-\/1),e\_guess(1:kc-\/1),\textcolor{stringliteral}{"TDMA\_H"},e\_itt)}
\DoxyCodeLine{443               e\_guess(1:kc-\/1) = e\_itt(1:kc-\/1) / sqrt(sum(e\_itt(1:kc-\/1)**2))}
\DoxyCodeLine{444 \textcolor{keywordflow}{            enddo} \textcolor{comment}{! itt-\/loop}}
\DoxyCodeLine{445             w\_strct(2:kc) = e\_guess(1:kc-\/1)}
\DoxyCodeLine{446             w\_strct(1)    = 0.0 \textcolor{comment}{! rigid lid at surface}}
\DoxyCodeLine{447             w\_strct(kc+1) = 0.0 \textcolor{comment}{! zero-\/flux at bottom}}
\DoxyCodeLine{448 }
\DoxyCodeLine{449             \textcolor{comment}{! Check to see if solver worked}}
\DoxyCodeLine{450             ig\_stop = 0 ; jg\_stop = 0}
\DoxyCodeLine{451             \textcolor{keywordflow}{if} (isnan(sum(w\_strct(1:kc+1))))\textcolor{keywordflow}{then}}
\DoxyCodeLine{452               print *, \textcolor{stringliteral}{"Wave\_structure: w\_strct has a NAN at ig="}, ig, \textcolor{stringliteral}{", jg="}, jg}
\DoxyCodeLine{453               \textcolor{keywordflow}{if} (i<g\%isc .or. i>g\%iec .or. j<g\%jsc .or. j>g\%jec)\textcolor{keywordflow}{then}}
\DoxyCodeLine{454                 print *, \textcolor{stringliteral}{"This is occuring at a halo point."}}
\DoxyCodeLine{455 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{456               ig\_stop = ig ; jg\_stop = jg}
\DoxyCodeLine{457 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{458 }
\DoxyCodeLine{459             \textcolor{comment}{! Normalize vertical structure function of w such that}}
\DoxyCodeLine{460             \textcolor{comment}{! \(\backslash\)int(w\_strct)\string^2dz = a\_int (a\_int could be any value, e.g., 0.5)}}
\DoxyCodeLine{461             nzm = kc+1 \textcolor{comment}{! number of layer interfaces after merging}}
\DoxyCodeLine{462                        \textcolor{comment}{!(including surface and bottom)}}
\DoxyCodeLine{463             w2avg = 0.0}
\DoxyCodeLine{464             \textcolor{keywordflow}{do} k=1,nzm-\/1}
\DoxyCodeLine{465               dz(k) = us\%Z\_to\_m*hc(k)}
\DoxyCodeLine{466               w2avg = w2avg + 0.5*(w\_strct(k)**2+w\_strct(k+1)**2)*dz(k)}
\DoxyCodeLine{467 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{468             \textcolor{comment}{!\#\#\# Some mathematical cancellations could occur in the next two lines.}}
\DoxyCodeLine{469             w2avg = w2avg / htot(i,j)}
\DoxyCodeLine{470             w\_strct(:) = w\_strct(:) / sqrt(htot(i,j)*w2avg*i\_a\_int)}
\DoxyCodeLine{471 }
\DoxyCodeLine{472             \textcolor{comment}{! Calculate vertical structure function of u (i.e. dw/dz)}}
\DoxyCodeLine{473             \textcolor{keywordflow}{do} k=2,nzm-\/1}
\DoxyCodeLine{474               u\_strct(k) = 0.5*((w\_strct(k-\/1) -\/ w\_strct(k)  )/dz(k-\/1) + \&}
\DoxyCodeLine{475                                 (w\_strct(k)   -\/ w\_strct(k+1))/dz(k))}
\DoxyCodeLine{476 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{477             u\_strct(1)   = (w\_strct(1)   -\/  w\_strct(2) )/dz(1)}
\DoxyCodeLine{478             u\_strct(nzm) = (w\_strct(nzm-\/1)-\/  w\_strct(nzm))/dz(nzm-\/1)}
\DoxyCodeLine{479 }
\DoxyCodeLine{480             \textcolor{comment}{! Calculate wavenumber magnitude}}
\DoxyCodeLine{481             f2 = g\%CoriolisBu(i,j)**2}
\DoxyCodeLine{482             \textcolor{comment}{!f2 = 0.25*US\%s\_to\_T**2 *((G\%CoriolisBu(I,J)**2 + G\%CoriolisBu(I-\/1,J-\/1)**2) + \&}}
\DoxyCodeLine{483             \textcolor{comment}{!    (G\%CoriolisBu(I,J-\/1)**2 + G\%CoriolisBu(I-\/1,J)**2))}}
\DoxyCodeLine{484             kmag2 = us\%m\_to\_L**2 * (freq**2 -\/ f2) / (cn(i,j)**2 + cg\_subro**2)}
\DoxyCodeLine{485 }
\DoxyCodeLine{486             \textcolor{comment}{! Calculate terms in vertically integrated energy equation}}
\DoxyCodeLine{487             int\_dwdz2 = 0.0 ; int\_w2 = 0.0 ; int\_n2w2 = 0.0}
\DoxyCodeLine{488             u\_strct2(1:nzm) = u\_strct(1:nzm)**2}
\DoxyCodeLine{489             w\_strct2(1:nzm) = w\_strct(1:nzm)**2}
\DoxyCodeLine{490             \textcolor{comment}{! vertical integration with Trapezoidal rule}}
\DoxyCodeLine{491             \textcolor{keywordflow}{do} k=1,nzm-\/1}
\DoxyCodeLine{492               int\_dwdz2 = int\_dwdz2 + 0.5*(u\_strct2(k)+u\_strct2(k+1)) * us\%m\_to\_Z*dz(k)}
\DoxyCodeLine{493               int\_w2    = int\_w2    + 0.5*(w\_strct2(k)+w\_strct2(k+1)) * us\%m\_to\_Z*dz(k)}
\DoxyCodeLine{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)}
\DoxyCodeLine{495 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{496 }
\DoxyCodeLine{497             \textcolor{comment}{! Back-\/calculate amplitude from energy equation}}
\DoxyCodeLine{498             \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(en) .and. (freq**2*kmag2 > 0.0)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{499               \textcolor{comment}{! Units here are [R}}
\DoxyCodeLine{500               ke\_term = 0.25*gv\%Rho0*( ((freq**2 + f2) / (freq**2*kmag2))*int\_dwdz2 + int\_w2 )}
\DoxyCodeLine{501               pe\_term = 0.25*gv\%Rho0*( int\_n2w2 / (us\%s\_to\_T*freq)**2 )}
\DoxyCodeLine{502               \textcolor{keywordflow}{if} (en(i,j) >= 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{503                 w0 = sqrt( en(i,j) / (ke\_term + pe\_term) )}
\DoxyCodeLine{504               \textcolor{keywordflow}{else}}
\DoxyCodeLine{505                 \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_structure: En < 0.0; setting to W0 to 0.0"})}
\DoxyCodeLine{506                 print *, \textcolor{stringliteral}{"En(i,j)="}, en(i,j), \textcolor{stringliteral}{" at ig="}, ig, \textcolor{stringliteral}{", jg="}, jg}
\DoxyCodeLine{507                 w0 = 0.0}
\DoxyCodeLine{508 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{509               \textcolor{comment}{! Calculate actual vertical velocity profile and derivative}}
\DoxyCodeLine{510               w\_profile(:)    = w0*w\_strct(:)}
\DoxyCodeLine{511               \textcolor{comment}{! dWdz\_profile(:) = W0*u\_strct(:)}}
\DoxyCodeLine{512               \textcolor{comment}{! Calculate average magnitude of actual horizontal velocity over a period}}
\DoxyCodeLine{513               uavg\_profile(:) = us\%Z\_to\_L*abs(w0*u\_strct(:)) * sqrt((freq**2 + f2) / (2.0*freq**2*kmag2))}
\DoxyCodeLine{514             \textcolor{keywordflow}{else}}
\DoxyCodeLine{515               w\_profile(:)    = 0.0}
\DoxyCodeLine{516               \textcolor{comment}{! dWdz\_profile(:) = 0.0}}
\DoxyCodeLine{517               uavg\_profile(:) = 0.0}
\DoxyCodeLine{518 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{519 }
\DoxyCodeLine{520             \textcolor{comment}{! Store values in control structure}}
\DoxyCodeLine{521             cs\%w\_strct(i,j,1:nzm)     = w\_strct(1:nzm)}
\DoxyCodeLine{522             cs\%u\_strct(i,j,1:nzm)     = u\_strct(1:nzm)}
\DoxyCodeLine{523             cs\%W\_profile(i,j,1:nzm)   = w\_profile(1:nzm)}
\DoxyCodeLine{524             cs\%Uavg\_profile(i,j,1:nzm)= uavg\_profile(1:nzm)}
\DoxyCodeLine{525             cs\%z\_depths(i,j,1:nzm)    = us\%Z\_to\_m*z\_int(1:nzm)}
\DoxyCodeLine{526             cs\%N2(i,j,1:nzm)          = n2(1:nzm)}
\DoxyCodeLine{527             cs\%num\_intfaces(i,j)      = nzm}
\DoxyCodeLine{528           \textcolor{keywordflow}{else}}
\DoxyCodeLine{529             \textcolor{comment}{! If not enough layers, default to zero}}
\DoxyCodeLine{530             nzm = kc+1}
\DoxyCodeLine{531             cs\%w\_strct(i,j,1:nzm)     = 0.0}
\DoxyCodeLine{532             cs\%u\_strct(i,j,1:nzm)     = 0.0}
\DoxyCodeLine{533             cs\%W\_profile(i,j,1:nzm)   = 0.0}
\DoxyCodeLine{534             cs\%Uavg\_profile(i,j,1:nzm)= 0.0}
\DoxyCodeLine{535             cs\%z\_depths(i,j,1:nzm)    = 0.0 \textcolor{comment}{! could use actual values}}
\DoxyCodeLine{536             cs\%N2(i,j,1:nzm)          = 0.0 \textcolor{comment}{! could use with actual values}}
\DoxyCodeLine{537             cs\%num\_intfaces(i,j)       = nzm}
\DoxyCodeLine{538 \textcolor{keywordflow}{          endif}  \textcolor{comment}{! kc >= 3 and kc > ModeNum + 1?}}
\DoxyCodeLine{539 \textcolor{keywordflow}{        endif} \textcolor{comment}{! drxh\_sum >= 0?}}
\DoxyCodeLine{540       \textcolor{comment}{!else     ! if at test point -\/ delete later}}
\DoxyCodeLine{541       \textcolor{comment}{!  return ! if at test point -\/ delete later}}
\DoxyCodeLine{542       \textcolor{comment}{!endif    ! if at test point -\/ delete later}}
\DoxyCodeLine{543 \textcolor{keywordflow}{      endif} \textcolor{comment}{! mask2dT > 0.5?}}
\DoxyCodeLine{544     \textcolor{keywordflow}{else}}
\DoxyCodeLine{545       \textcolor{comment}{! if cn=0.0, default to zero}}
\DoxyCodeLine{546       nzm                       = nz+1\textcolor{comment}{! could use actual values}}
\DoxyCodeLine{547       cs\%w\_strct(i,j,1:nzm)     = 0.0}
\DoxyCodeLine{548       cs\%u\_strct(i,j,1:nzm)     = 0.0}
\DoxyCodeLine{549       cs\%W\_profile(i,j,1:nzm)   = 0.0}
\DoxyCodeLine{550       cs\%Uavg\_profile(i,j,1:nzm)= 0.0}
\DoxyCodeLine{551       cs\%z\_depths(i,j,1:nzm)    = 0.0 \textcolor{comment}{! could use actual values}}
\DoxyCodeLine{552       cs\%N2(i,j,1:nzm)          = 0.0 \textcolor{comment}{! could use with actual values}}
\DoxyCodeLine{553       cs\%num\_intfaces(i,j)       = nzm}
\DoxyCodeLine{554 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! if cn>0.0? ; i-\/loop}}
\DoxyCodeLine{555 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-\/loop}}
\DoxyCodeLine{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}}
\doxysubsubsection{\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(\mbox{\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{\texttt{ in}}  & {\em time} & The current model time. \\
\hline
\mbox{\texttt{ in}}  & {\em g} & The ocean\textquotesingle{}s grid structure. \\
\hline
\mbox{\texttt{ in}}  & {\em param\+\_\+file} & A structure to parse for run-\/time parameters. \\
\hline
\mbox{\texttt{ 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}{0}
\DoxyCodeLine{687   \textcolor{keywordtype}{type}(time\_type),         \textcolor{keywordtype}{intent(in)} :: Time\textcolor{comment}{ !< The current model time.}}
\DoxyCodeLine{688   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)} :: G\textcolor{comment}{    !< The ocean's grid structure.}}
\DoxyCodeLine{689   \textcolor{keywordtype}{type}(param\_file\_type),   \textcolor{keywordtype}{intent(in)} :: param\_file\textcolor{comment}{ !< A structure to parse for run-\/time}}
\DoxyCodeLine{690 \textcolor{comment}{                                              !! parameters.}}
\DoxyCodeLine{691   \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)} :: diag\textcolor{comment}{ !< A structure that is used to regulate}}
\DoxyCodeLine{692 \textcolor{comment}{                                              !! diagnostic output.}}
\DoxyCodeLine{693   \textcolor{keywordtype}{type}(wave\_structure\_CS), \textcolor{keywordtype}{pointer}    :: CS\textcolor{comment}{   !< A pointer that is set to point to the}}
\DoxyCodeLine{694 \textcolor{comment}{                                              !! control structure for this module.}}
\DoxyCodeLine{695 \textcolor{comment}{! This include declares and sets the variable "version".}}
\DoxyCodeLine{696 \textcolor{preprocessor}{\#include "version\_variable.h"}}
\DoxyCodeLine{697 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_wave\_structure"}  \textcolor{comment}{! This module's name.}}
\DoxyCodeLine{698   \textcolor{keywordtype}{integer} :: isd, ied, jsd, jed, nz}
\DoxyCodeLine{699 }
\DoxyCodeLine{700   isd = g\%isd ; ied = g\%ied ; jsd = g\%jsd ; jed = g\%jed ; nz = g\%ke}
\DoxyCodeLine{701 }
\DoxyCodeLine{702   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{703     \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_structure\_init called with an "}// \&}
\DoxyCodeLine{704                             \textcolor{stringliteral}{"associated control structure."})}
\DoxyCodeLine{705     \textcolor{keywordflow}{return}}
\DoxyCodeLine{706   \textcolor{keywordflow}{else} ; \textcolor{keyword}{allocate}(cs) ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{707 }
\DoxyCodeLine{708   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"INTERNAL\_TIDE\_SOURCE\_X"}, cs\%int\_tide\_source\_x, \&}
\DoxyCodeLine{709                  \textcolor{stringliteral}{"X Location of generation site for internal tide"}, default=1.)}
\DoxyCodeLine{710   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"INTERNAL\_TIDE\_SOURCE\_Y"}, cs\%int\_tide\_source\_y, \&}
\DoxyCodeLine{711                  \textcolor{stringliteral}{"Y Location of generation site for internal tide"}, default=1.)}
\DoxyCodeLine{712 }
\DoxyCodeLine{713   cs\%diag => diag}
\DoxyCodeLine{714 }
\DoxyCodeLine{715   \textcolor{comment}{! Allocate memory for variable in control structure; note,}}
\DoxyCodeLine{716   \textcolor{comment}{! not all rows will be filled if layers get merged!}}
\DoxyCodeLine{717   \textcolor{keyword}{allocate}(cs\%w\_strct(isd:ied,jsd:jed,nz+1))}
\DoxyCodeLine{718   \textcolor{keyword}{allocate}(cs\%u\_strct(isd:ied,jsd:jed,nz+1))}
\DoxyCodeLine{719   \textcolor{keyword}{allocate}(cs\%W\_profile(isd:ied,jsd:jed,nz+1))}
\DoxyCodeLine{720   \textcolor{keyword}{allocate}(cs\%Uavg\_profile(isd:ied,jsd:jed,nz+1))}
\DoxyCodeLine{721   \textcolor{keyword}{allocate}(cs\%z\_depths(isd:ied,jsd:jed,nz+1))}
\DoxyCodeLine{722   \textcolor{keyword}{allocate}(cs\%N2(isd:ied,jsd:jed,nz+1))}
\DoxyCodeLine{723   \textcolor{keyword}{allocate}(cs\%num\_intfaces(isd:ied,jsd:jed))}
\DoxyCodeLine{724 }
\DoxyCodeLine{725   \textcolor{comment}{! Write all relevant parameters to the model log.}}
\DoxyCodeLine{726   \textcolor{keyword}{call }log\_version(param\_file, mdl, version, \textcolor{stringliteral}{""})}
\DoxyCodeLine{727 }

\end{DoxyCode}
