\hypertarget{namespacemom__wave__speed}{}\doxysection{mom\+\_\+wave\+\_\+speed Module Reference}
\label{namespacemom__wave__speed}\index{mom\_wave\_speed@{mom\_wave\_speed}}


\doxysubsection{Detailed Description}
Routines for calculating baroclinic wave speeds. \doxysubsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \mbox{\hyperlink{structmom__wave__speed_1_1wave__speed__cs}{wave\+\_\+speed\+\_\+cs}}
\begin{DoxyCompactList}\small\item\em Control structure for M\+O\+M\+\_\+wave\+\_\+speed. \end{DoxyCompactList}\end{DoxyCompactItemize}
\doxysubsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \mbox{\hyperlink{namespacemom__wave__speed_a936732268d9f4097149adb82b393cf39}{wave\+\_\+speed}} (h, tv, G, GV, US, cg1, CS, full\+\_\+halos, use\+\_\+ebt\+\_\+mode, mono\+\_\+\+N2\+\_\+column\+\_\+fraction, mono\+\_\+\+N2\+\_\+depth, modal\+\_\+structure, better\+\_\+speed\+\_\+est, min\+\_\+speed, wave\+\_\+speed\+\_\+tol)
\begin{DoxyCompactList}\small\item\em Calculates the wave speed of the first baroclinic mode. \end{DoxyCompactList}\item 
subroutine \mbox{\hyperlink{namespacemom__wave__speed_a3acab282f654fdcf4a7f48f115f4047b}{tdma6}} (n, a, c, lam, y)
\begin{DoxyCompactList}\small\item\em Solve a non-\/symmetric tridiagonal problem with the sum of the upper and lower diagnonals minus a scalar contribution as the leading diagonal. This uses the Thomas algorithm rather than the Hallberg algorithm since the matrix is not symmetric. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__wave__speed_a4e969f10d0cdd765f6804973a74ff479}{wave\+\_\+speeds}} (h, tv, G, GV, US, nmodes, cn, CS, full\+\_\+halos, better\+\_\+speed\+\_\+est, min\+\_\+speed, wave\+\_\+speed\+\_\+tol)
\begin{DoxyCompactList}\small\item\em Calculates the wave speeds for the first few barolinic modes. \end{DoxyCompactList}\item 
subroutine \mbox{\hyperlink{namespacemom__wave__speed_aa07774097409ee65a75722ad690c56e4}{tridiag\+\_\+det}} (a, c, ks, ke, lam, det, ddet, row\+\_\+scale)
\begin{DoxyCompactList}\small\item\em Calculate the determinant of a tridiagonal matrix with diagonals a,b-\/lam,c and its derivative with lam, where lam is constant across rows. Only the ratio of det to its derivative and their signs are typically used, so internal rescaling by consistent factors are used to avoid over-\/ or underflow. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__wave__speed_afa7284f32921f1e5d31530633cf95022}{wave\+\_\+speed\+\_\+init}} (CS, use\+\_\+ebt\+\_\+mode, mono\+\_\+\+N2\+\_\+column\+\_\+fraction, mono\+\_\+\+N2\+\_\+depth, remap\+\_\+answers\+\_\+2018, better\+\_\+speed\+\_\+est, min\+\_\+speed, wave\+\_\+speed\+\_\+tol)
\begin{DoxyCompactList}\small\item\em Initialize control structure for M\+O\+M\+\_\+wave\+\_\+speed. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__wave__speed_a13f3b425b43466d5af43f15d26c25a59}{wave\+\_\+speed\+\_\+set\+\_\+param}} (CS, use\+\_\+ebt\+\_\+mode, mono\+\_\+\+N2\+\_\+column\+\_\+fraction, mono\+\_\+\+N2\+\_\+depth, remap\+\_\+answers\+\_\+2018, better\+\_\+speed\+\_\+est, min\+\_\+speed, wave\+\_\+speed\+\_\+tol)
\begin{DoxyCompactList}\small\item\em Sets internal parameters for M\+O\+M\+\_\+wave\+\_\+speed. \end{DoxyCompactList}\end{DoxyCompactItemize}


\doxysubsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespacemom__wave__speed_a3acab282f654fdcf4a7f48f115f4047b}\label{namespacemom__wave__speed_a3acab282f654fdcf4a7f48f115f4047b}} 
\index{mom\_wave\_speed@{mom\_wave\_speed}!tdma6@{tdma6}}
\index{tdma6@{tdma6}!mom\_wave\_speed@{mom\_wave\_speed}}
\doxysubsubsection{\texorpdfstring{tdma6()}{tdma6()}}
{\footnotesize\ttfamily subroutine mom\+\_\+wave\+\_\+speed\+::tdma6 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{n,  }\item[{real, dimension(\+:), intent(in)}]{a,  }\item[{real, dimension(\+:), intent(in)}]{c,  }\item[{real, intent(in)}]{lam,  }\item[{real, dimension(\+:), intent(inout)}]{y }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Solve a non-\/symmetric tridiagonal problem with the sum of the upper and lower diagnonals minus a scalar contribution as the leading diagonal. This uses the Thomas algorithm rather than the Hallberg algorithm since the matrix is not symmetric. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\texttt{ in}}  & {\em n} & Number of rows of matrix \\
\hline
\mbox{\texttt{ in}}  & {\em a} & Lower diagonal \mbox{[}T2 L-\/2 $\sim$$>$ s2 m-\/2\mbox{]} \\
\hline
\mbox{\texttt{ in}}  & {\em c} & Upper diagonal \mbox{[}T2 L-\/2 $\sim$$>$ s2 m-\/2\mbox{]} \\
\hline
\mbox{\texttt{ in}}  & {\em lam} & Scalar subtracted from leading diagonal \mbox{[}T2 L-\/2 $\sim$$>$ s2 m-\/2\mbox{]} \\
\hline
\mbox{\texttt{ in,out}}  & {\em y} & R\+HS on entry, result on exit \\
\hline
\end{DoxyParams}


Definition at line 593 of file M\+O\+M\+\_\+wave\+\_\+speed.\+F90.


\begin{DoxyCode}{0}
\DoxyCodeLine{593   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of rows of matrix}}
\DoxyCodeLine{594 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}    :: a\textcolor{comment}{ !< Lower diagonal   [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{595 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}    :: c\textcolor{comment}{ !< Upper diagonal   [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{596 \textcolor{keywordtype}{  real},               \textcolor{keywordtype}{intent(in)}    :: lam\textcolor{comment}{ !< Scalar subtracted from leading diagonal [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{597 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(inout)} :: y\textcolor{comment}{ !< RHS on entry, result on exit}}
\DoxyCodeLine{598 }
\DoxyCodeLine{599   \textcolor{comment}{! Local variables}}
\DoxyCodeLine{600 \textcolor{keywordtype}{  real} :: lambda     \textcolor{comment}{! A temporary variable in [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{601 \textcolor{keywordtype}{  real} :: beta(n)    \textcolor{comment}{! A temporary variable in [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{602 \textcolor{keywordtype}{  real} :: I\_beta(n)  \textcolor{comment}{! A temporary variable in [L2 T-\/2 ~> m2 s-\/2]}}
\DoxyCodeLine{603 \textcolor{keywordtype}{  real} :: yy(n)      \textcolor{comment}{! A temporary variable with the same units as y on entry.}}
\DoxyCodeLine{604   \textcolor{keywordtype}{integer} :: k, m}
\DoxyCodeLine{605 }
\DoxyCodeLine{606   lambda = lam}
\DoxyCodeLine{607   beta(1) = (a(1)+c(1)) -\/ lambda}
\DoxyCodeLine{608   \textcolor{keywordflow}{if} (beta(1)==0.) \textcolor{keywordflow}{then} \textcolor{comment}{! lam was chosen too perfectly}}
\DoxyCodeLine{609     \textcolor{comment}{! Change lambda and redo this first row}}
\DoxyCodeLine{610     lambda = (1. + 1.e-\/5) * lambda}
\DoxyCodeLine{611     beta(1) = (a(1)+c(1)) -\/ lambda}
\DoxyCodeLine{612 \textcolor{keywordflow}{  endif}}
\DoxyCodeLine{613   i\_beta(1) = 1. / beta(1)}
\DoxyCodeLine{614   yy(1) = y(1)}
\DoxyCodeLine{615   \textcolor{keywordflow}{do} k = 2, n}
\DoxyCodeLine{616     beta(k) = ( (a(k)+c(k)) -\/ lambda ) -\/ a(k) * c(k-\/1) * i\_beta(k-\/1)}
\DoxyCodeLine{617     \textcolor{comment}{! Perhaps the following 0 needs to become a tolerance to handle underflow?}}
\DoxyCodeLine{618     \textcolor{keywordflow}{if} (beta(k)==0.) \textcolor{keywordflow}{then} \textcolor{comment}{! lam was chosen too perfectly}}
\DoxyCodeLine{619       \textcolor{comment}{! Change lambda and redo everything up to row k}}
\DoxyCodeLine{620       lambda = (1. + 1.e-\/5) * lambda}
\DoxyCodeLine{621       i\_beta(1) = 1. / ( (a(1)+c(1)) -\/ lambda )}
\DoxyCodeLine{622       \textcolor{keywordflow}{do} m = 2, k}
\DoxyCodeLine{623         i\_beta(m) = 1. / ( ( (a(m)+c(m)) -\/ lambda ) -\/ a(m) * c(m-\/1) * i\_beta(m-\/1) )}
\DoxyCodeLine{624         yy(m) = y(m) + a(m) * yy(m-\/1) * i\_beta(m-\/1)}
\DoxyCodeLine{625 \textcolor{keywordflow}{      enddo}}
\DoxyCodeLine{626     \textcolor{keywordflow}{else}}
\DoxyCodeLine{627       i\_beta(k) = 1. / beta(k)}
\DoxyCodeLine{628 \textcolor{keywordflow}{    endif}}
\DoxyCodeLine{629     yy(k) = y(k) + a(k) * yy(k-\/1) * i\_beta(k-\/1)}
\DoxyCodeLine{630 \textcolor{keywordflow}{  enddo}}
\DoxyCodeLine{631   \textcolor{comment}{! The units of y change by a factor of [L2 T-\/2] in the following lines.}}
\DoxyCodeLine{632   y(n) = yy(n) * i\_beta(n)}
\DoxyCodeLine{633   \textcolor{keywordflow}{do} k = n-\/1, 1, -\/1}
\DoxyCodeLine{634     y(k) = ( yy(k) + c(k) * y(k+1) ) * i\_beta(k)}
\DoxyCodeLine{635 \textcolor{keywordflow}{  enddo}}
\DoxyCodeLine{636 }

\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__wave__speed_aa07774097409ee65a75722ad690c56e4}\label{namespacemom__wave__speed_aa07774097409ee65a75722ad690c56e4}} 
\index{mom\_wave\_speed@{mom\_wave\_speed}!tridiag\_det@{tridiag\_det}}
\index{tridiag\_det@{tridiag\_det}!mom\_wave\_speed@{mom\_wave\_speed}}
\doxysubsubsection{\texorpdfstring{tridiag\_det()}{tridiag\_det()}}
{\footnotesize\ttfamily subroutine mom\+\_\+wave\+\_\+speed\+::tridiag\+\_\+det (\begin{DoxyParamCaption}\item[{real, dimension(\+:), intent(in)}]{a,  }\item[{real, dimension(\+:), intent(in)}]{c,  }\item[{integer, intent(in)}]{ks,  }\item[{integer, intent(in)}]{ke,  }\item[{real, intent(in)}]{lam,  }\item[{real, intent(out)}]{det,  }\item[{real, intent(out)}]{ddet,  }\item[{real, intent(in), optional}]{row\+\_\+scale }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Calculate the determinant of a tridiagonal matrix with diagonals a,b-\/lam,c and its derivative with lam, where lam is constant across rows. Only the ratio of det to its derivative and their signs are typically used, so internal rescaling by consistent factors are used to avoid over-\/ or underflow. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\texttt{ in}}  & {\em a} & Lower diagonal of matrix (first entry unused) \\
\hline
\mbox{\texttt{ in}}  & {\em c} & Upper diagonal of matrix (last entry unused) \\
\hline
\mbox{\texttt{ in}}  & {\em ks} & Starting index to use in determinant \\
\hline
\mbox{\texttt{ in}}  & {\em ke} & Ending index to use in determinant \\
\hline
\mbox{\texttt{ in}}  & {\em lam} & Value subtracted from b \\
\hline
\mbox{\texttt{ out}}  & {\em det} & Determinant \\
\hline
\mbox{\texttt{ out}}  & {\em ddet} & Derivative of determinant with lam \\
\hline
\mbox{\texttt{ in}}  & {\em row\+\_\+scale} & A scaling factor of the rows of the matrix to limit the growth of the determinant \\
\hline
\end{DoxyParams}


Definition at line 1146 of file M\+O\+M\+\_\+wave\+\_\+speed.\+F90.


\begin{DoxyCode}{0}
\DoxyCodeLine{1146 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: a\textcolor{comment}{     !< Lower diagonal of matrix (first entry unused)}}
\DoxyCodeLine{1147 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: c\textcolor{comment}{     !< Upper diagonal of matrix (last entry unused)}}
\DoxyCodeLine{1148   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)} :: ks\textcolor{comment}{    !< Starting index to use in determinant}}
\DoxyCodeLine{1149   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)} :: ke\textcolor{comment}{    !< Ending index to use in determinant}}
\DoxyCodeLine{1150 \textcolor{keywordtype}{  real},               \textcolor{keywordtype}{intent(in)} :: lam\textcolor{comment}{   !< Value subtracted from b}}
\DoxyCodeLine{1151 \textcolor{keywordtype}{  real},               \textcolor{keywordtype}{intent(out)}:: det\textcolor{comment}{   !< Determinant}}
\DoxyCodeLine{1152 \textcolor{keywordtype}{  real},               \textcolor{keywordtype}{intent(out)}:: ddet\textcolor{comment}{  !< Derivative of determinant with lam}}
\DoxyCodeLine{1153 \textcolor{keywordtype}{  real},     \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: row\_scale\textcolor{comment}{ !< A scaling factor of the rows of the}}
\DoxyCodeLine{1154 \textcolor{comment}{                                      !! matrix to limit the growth of the determinant}}
\DoxyCodeLine{1155   \textcolor{comment}{! Local variables}}
\DoxyCodeLine{1156 \textcolor{keywordtype}{  real} :: detKm1, detKm2   \textcolor{comment}{! Cumulative value of the determinant for the previous two layers.}}
\DoxyCodeLine{1157 \textcolor{keywordtype}{  real} :: ddetKm1, ddetKm2 \textcolor{comment}{! Derivative of the cumulative determinant with lam for the previous two layers.}}
\DoxyCodeLine{1158 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{parameter} :: rescale = 1024.0**4 \textcolor{comment}{! max value of determinant allowed before rescaling}}
\DoxyCodeLine{1159 \textcolor{keywordtype}{  real} :: rscl      \textcolor{comment}{! A rescaling factor that is applied succesively to each row.}}
\DoxyCodeLine{1160 \textcolor{keywordtype}{  real} :: I\_rescale \textcolor{comment}{! inverse of rescale}}
\DoxyCodeLine{1161   \textcolor{keywordtype}{integer} :: k      \textcolor{comment}{! row (layer interface) index}}
\DoxyCodeLine{1162 }
\DoxyCodeLine{1163   i\_rescale = 1.0 / rescale}
\DoxyCodeLine{1164   rscl = 1.0 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(row\_scale)) rscl = row\_scale}
\DoxyCodeLine{1165 }
\DoxyCodeLine{1166   detkm1 = 1.0 ; ddetkm1 = 0.0}
\DoxyCodeLine{1167   det = (a(ks)+c(ks)) -\/ lam ; ddet = -\/1.0}
\DoxyCodeLine{1168   \textcolor{keywordflow}{do} k=ks+1,ke}
\DoxyCodeLine{1169     \textcolor{comment}{! Shift variables and rescale rows to avoid over-\/ or underflow.}}
\DoxyCodeLine{1170     detkm2 = row\_scale*detkm1 ; ddetkm2 = row\_scale*ddetkm1}
\DoxyCodeLine{1171     detkm1 = row\_scale*det    ; ddetkm1 = row\_scale*ddet}
\DoxyCodeLine{1172 }
\DoxyCodeLine{1173     det =  ((a(k)+c(k))-\/lam)*detkm1  -\/ (a(k)*c(k-\/1))*detkm2}
\DoxyCodeLine{1174     ddet = ((a(k)+c(k))-\/lam)*ddetkm1 -\/ (a(k)*c(k-\/1))*ddetkm2 -\/ detkm1}
\DoxyCodeLine{1175 }
\DoxyCodeLine{1176     \textcolor{comment}{! Rescale det \& ddet if det is getting too large or too small.}}
\DoxyCodeLine{1177     \textcolor{keywordflow}{if} (abs(det) > rescale) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1178       det = i\_rescale*det ; detkm1 = i\_rescale*detkm1}
\DoxyCodeLine{1179       ddet = i\_rescale*ddet ; ddetkm1 = i\_rescale*ddetkm1}
\DoxyCodeLine{1180     \textcolor{keywordflow}{elseif} (abs(det) < i\_rescale) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1181       det = rescale*det ; detkm1 = rescale*detkm1}
\DoxyCodeLine{1182       ddet = rescale*ddet ; ddetkm1 = rescale*ddetkm1}
\DoxyCodeLine{1183 \textcolor{keywordflow}{    endif}}
\DoxyCodeLine{1184 \textcolor{keywordflow}{  enddo}}
\DoxyCodeLine{1185 }

\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__wave__speed_a936732268d9f4097149adb82b393cf39}\label{namespacemom__wave__speed_a936732268d9f4097149adb82b393cf39}} 
\index{mom\_wave\_speed@{mom\_wave\_speed}!wave\_speed@{wave\_speed}}
\index{wave\_speed@{wave\_speed}!mom\_wave\_speed@{mom\_wave\_speed}}
\doxysubsubsection{\texorpdfstring{wave\_speed()}{wave\_speed()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+speed\+::wave\+\_\+speed (\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(out)}]{cg1,  }\item[{type(\mbox{\hyperlink{structmom__wave__speed_1_1wave__speed__cs}{wave\+\_\+speed\+\_\+cs}}), pointer}]{CS,  }\item[{logical, intent(in), optional}]{full\+\_\+halos,  }\item[{logical, intent(in), optional}]{use\+\_\+ebt\+\_\+mode,  }\item[{real, intent(in), optional}]{mono\+\_\+\+N2\+\_\+column\+\_\+fraction,  }\item[{real, intent(in), optional}]{mono\+\_\+\+N2\+\_\+depth,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(out), optional}]{modal\+\_\+structure,  }\item[{logical, intent(in), optional}]{better\+\_\+speed\+\_\+est,  }\item[{real, intent(in), optional}]{min\+\_\+speed,  }\item[{real, intent(in), optional}]{wave\+\_\+speed\+\_\+tol }\end{DoxyParamCaption})}



Calculates the wave speed of the first baroclinic mode. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\texttt{ in}}  & {\em g} & Ocean grid structure \\
\hline
\mbox{\texttt{ in}}  & {\em gv} & Vertical grid structure \\
\hline
\mbox{\texttt{ in}}  & {\em us} & A dimensional unit scaling type \\
\hline
\mbox{\texttt{ in}}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]} \\
\hline
\mbox{\texttt{ in}}  & {\em tv} & Thermodynamic variables \\
\hline
\mbox{\texttt{ out}}  & {\em cg1} & First mode internal wave speed \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]} \\
\hline
 & {\em cs} & Control structure for M\+O\+M\+\_\+wave\+\_\+speed \\
\hline
\mbox{\texttt{ in}}  & {\em full\+\_\+halos} & If true, do the calculation over the entire computational domain. \\
\hline
\mbox{\texttt{ in}}  & {\em use\+\_\+ebt\+\_\+mode} & If true, use the equivalent barotropic mode instead of the first baroclinic mode. \\
\hline
\mbox{\texttt{ in}}  & {\em mono\+\_\+n2\+\_\+column\+\_\+fraction} & The lower fraction of water column over which N2 is limited as monotonic for the purposes of calculating vertical modal structure. \\
\hline
\mbox{\texttt{ in}}  & {\em mono\+\_\+n2\+\_\+depth} & A depth below which N2 is limited as monotonic for the purposes of calculating vertical modal structure \mbox{[}Z $\sim$$>$ m\mbox{]}. \\
\hline
\mbox{\texttt{ out}}  & {\em modal\+\_\+structure} & Normalized model structure \mbox{[}nondim\mbox{]} \\
\hline
\mbox{\texttt{ in}}  & {\em better\+\_\+speed\+\_\+est} & If true, use a more robust estimate of the first mode speed as the starting point for iterations. \\
\hline
\mbox{\texttt{ in}}  & {\em min\+\_\+speed} & If present, set a floor in the first mode speed below which 0 is returned \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}. \\
\hline
\mbox{\texttt{ in}}  & {\em wave\+\_\+speed\+\_\+tol} & The fractional tolerance for finding the wave speeds \mbox{[}nondim\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 59 of file M\+O\+M\+\_\+wave\+\_\+speed.\+F90.


\begin{DoxyCode}{0}
\DoxyCodeLine{59   \textcolor{keywordtype}{type}(ocean\_grid\_type),            \textcolor{keywordtype}{intent(in)}  :: G\textcolor{comment}{  !< Ocean grid structure}}
\DoxyCodeLine{60   \textcolor{keywordtype}{type}(verticalGrid\_type),          \textcolor{keywordtype}{intent(in)}  :: GV\textcolor{comment}{ !< Vertical grid structure}}
\DoxyCodeLine{61   \textcolor{keywordtype}{type}(unit\_scale\_type),            \textcolor{keywordtype}{intent(in)}  :: US\textcolor{comment}{ !< A dimensional unit scaling type}}
\DoxyCodeLine{62 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \&}
\DoxyCodeLine{63                                     \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{  !< Layer thickness [H ~> m or kg m-\/2]}}
\DoxyCodeLine{64   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),            \textcolor{keywordtype}{intent(in)}  :: tv\textcolor{comment}{ !< Thermodynamic variables}}
\DoxyCodeLine{65 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))}, \textcolor{keywordtype}{intent(out)} :: cg1\textcolor{comment}{ !< First mode internal wave speed [L T-\/1 ~> m s-\/1]}}
\DoxyCodeLine{66   \textcolor{keywordtype}{type}(wave\_speed\_CS),              \textcolor{keywordtype}{pointer}     :: CS\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}}
\DoxyCodeLine{67   \textcolor{keywordtype}{logical},                \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: full\_halos\textcolor{comment}{ !< If true, do the calculation}}
\DoxyCodeLine{68 \textcolor{comment}{                                          !! over the entire computational domain.}}
\DoxyCodeLine{69   \textcolor{keywordtype}{logical},                \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: use\_ebt\_mode\textcolor{comment}{ !< If true, use the equivalent}}
\DoxyCodeLine{70 \textcolor{comment}{                                          !! barotropic mode instead of the first baroclinic mode.}}
\DoxyCodeLine{71 \textcolor{keywordtype}{  real},                   \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: mono\_N2\_column\_fraction\textcolor{comment}{ !< The lower fraction}}
\DoxyCodeLine{72 \textcolor{comment}{                                          !! of water column over which N2 is limited as monotonic}}
\DoxyCodeLine{73 \textcolor{comment}{                                          !! for the purposes of calculating vertical modal structure.}}
\DoxyCodeLine{74 \textcolor{keywordtype}{  real},                   \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: mono\_N2\_depth\textcolor{comment}{ !< A depth below which N2 is limited as}}
\DoxyCodeLine{75 \textcolor{comment}{                                          !! monotonic for the purposes of calculating vertical}}
\DoxyCodeLine{76 \textcolor{comment}{                                          !! modal structure [Z ~> m].}}
\DoxyCodeLine{77 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \&}
\DoxyCodeLine{78                           \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: modal\_structure\textcolor{comment}{ !< Normalized model structure [nondim]}}
\DoxyCodeLine{79   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: better\_speed\_est\textcolor{comment}{ !< If true, use a more robust estimate of the first}}
\DoxyCodeLine{80 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}}
\DoxyCodeLine{81 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: min\_speed\textcolor{comment}{ !< If present, set a floor in the first mode speed}}
\DoxyCodeLine{82 \textcolor{comment}{                                     !! below which 0 is returned [L T-\/1 ~> m s-\/1].}}
\DoxyCodeLine{83 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}}
\DoxyCodeLine{84 \textcolor{comment}{                                     !! wave speeds [nondim]}}
\DoxyCodeLine{85 }
\DoxyCodeLine{86   \textcolor{comment}{! Local variables}}
\DoxyCodeLine{87 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: \&}
\DoxyCodeLine{88     dRho\_dT, \&    \textcolor{comment}{! Partial derivative of density with temperature [R degC-\/1 ~> kg m-\/3 degC-\/1]}}
\DoxyCodeLine{89     dRho\_dS, \&    \textcolor{comment}{! Partial derivative of density with salinity [R ppt-\/1 ~> kg m-\/3 ppt-\/1]}}
\DoxyCodeLine{90     pres, \&       \textcolor{comment}{! Interface pressure [R L2 T-\/2 ~> Pa]}}
\DoxyCodeLine{91     T\_int, \&      \textcolor{comment}{! Temperature interpolated to interfaces [degC]}}
\DoxyCodeLine{92     S\_int, \&      \textcolor{comment}{! Salinity interpolated to interfaces [ppt]}}
\DoxyCodeLine{93     H\_top, \&      \textcolor{comment}{! The distance of each filtered interface from the ocean surface [Z ~> m]}}
\DoxyCodeLine{94     H\_bot, \&      \textcolor{comment}{! The distance of each filtered interface from the bottom [Z ~> m]}}
\DoxyCodeLine{95     gprime        \textcolor{comment}{! The reduced gravity across each interface [L2 Z-\/1 T-\/2 ~> m s-\/2].}}
\DoxyCodeLine{96 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: \&}
\DoxyCodeLine{97     Igl, Igu      \textcolor{comment}{! The inverse of the reduced gravity across an interface times}}
\DoxyCodeLine{98                   \textcolor{comment}{! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-\/2 ~> s2 m-\/2].}}
\DoxyCodeLine{99 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G),SZI\_(G))} :: \&}
\DoxyCodeLine{100     Hf, \&         \textcolor{comment}{! Layer thicknesses after very thin layers are combined [Z ~> m]}}
\DoxyCodeLine{101     Tf, \&         \textcolor{comment}{! Layer temperatures after very thin layers are combined [degC]}}
\DoxyCodeLine{102     Sf, \&         \textcolor{comment}{! Layer salinities after very thin layers are combined [ppt]}}
\DoxyCodeLine{103     Rf            \textcolor{comment}{! Layer densities after very thin layers are combined [R ~> kg m-\/3]}}
\DoxyCodeLine{104 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: \&}
\DoxyCodeLine{105     Hc, \&         \textcolor{comment}{! A column of layer thicknesses after convective istabilities are removed [Z ~> m]}}
\DoxyCodeLine{106     Tc, \&         \textcolor{comment}{! A column of layer temperatures after convective istabilities are removed [degC]}}
\DoxyCodeLine{107     Sc, \&         \textcolor{comment}{! A column of layer salinites after convective istabilities are removed [ppt]}}
\DoxyCodeLine{108     Rc, \&         \textcolor{comment}{! A column of layer densities after convective istabilities are removed [R ~> kg m-\/3]}}
\DoxyCodeLine{109     Hc\_H          \textcolor{comment}{! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-\/2]}}
\DoxyCodeLine{110 \textcolor{keywordtype}{  real} :: I\_Htot  \textcolor{comment}{! The inverse of the total filtered thicknesses [Z ~> m]}}
\DoxyCodeLine{111 \textcolor{keywordtype}{  real} :: det, ddet, detKm1, detKm2, ddetKm1, ddetKm2}
\DoxyCodeLine{112 \textcolor{keywordtype}{  real} :: lam     \textcolor{comment}{! The eigenvalue [T2 L-\/2 ~> s m-\/1]}}
\DoxyCodeLine{113 \textcolor{keywordtype}{  real} :: dlam    \textcolor{comment}{! The change in estimates of the eigenvalue [T2 L-\/2 ~> s m-\/1]}}
\DoxyCodeLine{114 \textcolor{keywordtype}{  real} :: lam0    \textcolor{comment}{! The first guess of the eigenvalue [T2 L-\/2 ~> s m-\/1]}}
\DoxyCodeLine{115 \textcolor{keywordtype}{  real} :: min\_h\_frac \textcolor{comment}{! [nondim]}}
\DoxyCodeLine{116 \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{117 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: \&}
\DoxyCodeLine{118     htot, hmin, \&  \textcolor{comment}{! Thicknesses [Z ~> m]}}
\DoxyCodeLine{119     H\_here, \&      \textcolor{comment}{! A thickness [Z ~> m]}}
\DoxyCodeLine{120     HxT\_here, \&    \textcolor{comment}{! A layer integrated temperature [degC Z ~> degC m]}}
\DoxyCodeLine{121     HxS\_here, \&    \textcolor{comment}{! A layer integrated salinity [ppt Z ~> ppt m]}}
\DoxyCodeLine{122     HxR\_here       \textcolor{comment}{! A layer integrated density [R Z ~> kg m-\/2]}}
\DoxyCodeLine{123 \textcolor{keywordtype}{  real} :: speed2\_tot \textcolor{comment}{! overestimate of the mode-\/1 speed squared [L2 T-\/2 ~> m2 s-\/2]}}
\DoxyCodeLine{124 \textcolor{keywordtype}{  real} :: cg1\_min2 \textcolor{comment}{! A floor in the squared first mode speed below which 0 is returned [L2 T-\/2 ~> m2 s-\/2]}}
\DoxyCodeLine{125 \textcolor{keywordtype}{  real} :: I\_Hnew   \textcolor{comment}{! The inverse of a new layer thickness [Z-\/1 ~> m-\/1]}}
\DoxyCodeLine{126 \textcolor{keywordtype}{  real} :: drxh\_sum \textcolor{comment}{! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-\/2]}}
\DoxyCodeLine{127 \textcolor{keywordtype}{  real} :: L2\_to\_Z2 \textcolor{comment}{! A scaling factor squared from units of lateral distances to depths [Z2 L-\/2 ~> 1].}}
\DoxyCodeLine{128 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{pointer}, \textcolor{keywordtype}{dimension(:,:,:)} :: T => null(), s => null()}
\DoxyCodeLine{129 \textcolor{keywordtype}{  real} :: g\_Rho0   \textcolor{comment}{! G\_Earth/Rho0 [L2 T-\/2 Z-\/1 R-\/1 ~> m4 s-\/2 kg-\/1].}}
\DoxyCodeLine{130 \textcolor{keywordtype}{  real} :: c2\_scale \textcolor{comment}{! A scaling factor for wave speeds to help control the growth of the determinant}}
\DoxyCodeLine{131                    \textcolor{comment}{! and its derivative with lam between rows of the Thomas algorithm solver.  The}}
\DoxyCodeLine{132                    \textcolor{comment}{! exact value should not matter for the final result if it is an even power of 2.}}
\DoxyCodeLine{133 \textcolor{keywordtype}{  real} :: tol\_Hfrac \textcolor{comment}{! Layers that together are smaller than this fraction of}}
\DoxyCodeLine{134                     \textcolor{comment}{! the total water column can be merged for efficiency.}}
\DoxyCodeLine{135 \textcolor{keywordtype}{  real} :: tol\_solve \textcolor{comment}{! The fractional tolerance with which to solve for the wave speeds [nondim]}}
\DoxyCodeLine{136 \textcolor{keywordtype}{  real} :: tol\_merge \textcolor{comment}{! The fractional change in estimated wave speed that is allowed}}
\DoxyCodeLine{137                     \textcolor{comment}{! when deciding to merge layers in the calculation [nondim]}}
\DoxyCodeLine{138 \textcolor{keywordtype}{  real} :: rescale, I\_rescale}
\DoxyCodeLine{139   \textcolor{keywordtype}{integer} :: kf(SZI\_(G)) \textcolor{comment}{! The number of active layers after filtering.}}
\DoxyCodeLine{140   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: max\_itt = 10}
\DoxyCodeLine{141 \textcolor{keywordtype}{  real} :: lam\_it(max\_itt), det\_it(max\_itt), ddet\_it(max\_itt)}
\DoxyCodeLine{142   \textcolor{keywordtype}{logical} :: use\_EOS    \textcolor{comment}{! If true, density is calculated from T \& S using an equation of state.}}
\DoxyCodeLine{143   \textcolor{keywordtype}{logical} :: better\_est \textcolor{comment}{! If true, use an improved estimate of the first mode internal wave speed.}}
\DoxyCodeLine{144   \textcolor{keywordtype}{logical} :: merge      \textcolor{comment}{! If true, merge the current layer with the one above.}}
\DoxyCodeLine{145   \textcolor{keywordtype}{integer} :: kc         \textcolor{comment}{! The number of layers in the column after merging}}
\DoxyCodeLine{146   \textcolor{keywordtype}{integer} :: i, j, k, k2, itt, is, ie, js, je, nz}
\DoxyCodeLine{147 \textcolor{keywordtype}{  real} :: hw, sum\_hc}
\DoxyCodeLine{148 \textcolor{keywordtype}{  real} :: gp      \textcolor{comment}{! A limited local copy of gprime [L2 Z-\/1 T-\/2 ~> m s-\/2]}}
\DoxyCodeLine{149 \textcolor{keywordtype}{  real} :: N2min   \textcolor{comment}{! A minimum buoyancy frequency [T-\/2 ~> s-\/2]}}
\DoxyCodeLine{150   \textcolor{keywordtype}{logical} :: l\_use\_ebt\_mode, calc\_modal\_structure}
\DoxyCodeLine{151 \textcolor{keywordtype}{  real} :: l\_mono\_N2\_column\_fraction, l\_mono\_N2\_depth}
\DoxyCodeLine{152 \textcolor{keywordtype}{  real} :: mode\_struct(SZK\_(G)), ms\_min, ms\_max, ms\_sq}
\DoxyCodeLine{153 }
\DoxyCodeLine{154   is = g\%isc ; ie = g\%iec ; js = g\%jsc ; je = g\%jec ; nz = g\%ke}
\DoxyCodeLine{155 }
\DoxyCodeLine{156   \textcolor{keywordflow}{if} (.not. \textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_wave\_speed: "}// \&}
\DoxyCodeLine{157            \textcolor{stringliteral}{"Module must be initialized before it is used."})}
\DoxyCodeLine{158   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(full\_halos)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (full\_halos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{159     is = g\%isd ; ie = g\%ied ; js = g\%jsd ; je = g\%jed}
\DoxyCodeLine{160 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{161 }
\DoxyCodeLine{162   l2\_to\_z2 = us\%L\_to\_Z**2}
\DoxyCodeLine{163 }
\DoxyCodeLine{164   l\_use\_ebt\_mode = cs\%use\_ebt\_mode}
\DoxyCodeLine{165   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(use\_ebt\_mode)) l\_use\_ebt\_mode = use\_ebt\_mode}
\DoxyCodeLine{166   l\_mono\_n2\_column\_fraction = cs\%mono\_N2\_column\_fraction}
\DoxyCodeLine{167   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_column\_fraction)) l\_mono\_n2\_column\_fraction = mono\_n2\_column\_fraction}
\DoxyCodeLine{168   l\_mono\_n2\_depth = cs\%mono\_N2\_depth}
\DoxyCodeLine{169   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_depth)) l\_mono\_n2\_depth = mono\_n2\_depth}
\DoxyCodeLine{170   calc\_modal\_structure = l\_use\_ebt\_mode}
\DoxyCodeLine{171   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) calc\_modal\_structure = .true.}
\DoxyCodeLine{172   \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}}
\DoxyCodeLine{173     \textcolor{keywordflow}{do} k=1,nz; \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{174       modal\_structure(i,j,k) = 0.0}
\DoxyCodeLine{175 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{176 \textcolor{keywordflow}{  endif}}
\DoxyCodeLine{177 }
\DoxyCodeLine{178   s => tv\%S ; t => tv\%T}
\DoxyCodeLine{179   g\_rho0 = gv\%g\_Earth / gv\%Rho0}
\DoxyCodeLine{180   \textcolor{comment}{! Simplifying the following could change answers at roundoff.}}
\DoxyCodeLine{181   z\_to\_pres = gv\%Z\_to\_H * (gv\%H\_to\_RZ * gv\%g\_Earth)}
\DoxyCodeLine{182   use\_eos = \textcolor{keyword}{associated}(tv\%eqn\_of\_state)}
\DoxyCodeLine{183 }
\DoxyCodeLine{184   better\_est = cs\%better\_cg1\_est ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(better\_speed\_est)) better\_est = better\_speed\_est}
\DoxyCodeLine{185 }
\DoxyCodeLine{186   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{187     tol\_solve = cs\%wave\_speed\_tol ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(wave\_speed\_tol)) tol\_solve = wave\_speed\_tol}
\DoxyCodeLine{188     tol\_hfrac  = 0.1*tol\_solve ; tol\_merge = tol\_solve / real(nz)}
\DoxyCodeLine{189   \textcolor{keywordflow}{else}}
\DoxyCodeLine{190     tol\_solve = 0.001 ; tol\_hfrac  = 0.0001 ; tol\_merge = 0.001}
\DoxyCodeLine{191 \textcolor{keywordflow}{  endif}}
\DoxyCodeLine{192 }
\DoxyCodeLine{193   \textcolor{comment}{! The rescaling below can control the growth of the determinant provided that}}
\DoxyCodeLine{194   \textcolor{comment}{! (tol\_merge*cg1\_min2/c2\_scale > I\_rescale).  For default values, this suggests a stable lower}}
\DoxyCodeLine{195   \textcolor{comment}{! bound on min\_speed of sqrt(nz/(tol\_solve*rescale)) or 3e2/1024**2 = 2.9e-\/4 m/s for 90 layers.}}
\DoxyCodeLine{196   \textcolor{comment}{! The upper bound on the rate of increase in the determinant is g'H/c2\_scale < rescale or in the}}
\DoxyCodeLine{197   \textcolor{comment}{! worst possible oceanic case of g'H < 0.5*10m/s2*1e4m = 5.e4 m2/s2 < 1024**2*c2\_scale, suggesting}}
\DoxyCodeLine{198   \textcolor{comment}{! that c2\_scale can safely be set to 1/(16*1024**2), which would decrease the stable floor on}}
\DoxyCodeLine{199   \textcolor{comment}{! min\_speed to ~6.9e-\/8 m/s for 90 layers or 2.33e-\/7 m/s for 1000 layers.}}
\DoxyCodeLine{200   cg1\_min2 = cs\%min\_speed2 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(min\_speed)) cg1\_min2 = min\_speed**2}
\DoxyCodeLine{201   rescale = 1024.0**4 ; i\_rescale = 1.0/rescale}
\DoxyCodeLine{202   c2\_scale = us\%m\_s\_to\_L\_T**2 / 4096.0**2 \textcolor{comment}{! Other powers of 2 give identical results.}}
\DoxyCodeLine{203 }
\DoxyCodeLine{204   min\_h\_frac = tol\_hfrac / real(nz)}
\DoxyCodeLine{205 \textcolor{comment}{!\$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,US,min\_h\_frac,use\_EOS,T,S,tv,\&}}
\DoxyCodeLine{206 \textcolor{comment}{!\$OMP                                  calc\_modal\_structure,l\_use\_ebt\_mode,modal\_structure, \&}}
\DoxyCodeLine{207 \textcolor{comment}{!\$OMP                                  l\_mono\_N2\_column\_fraction,l\_mono\_N2\_depth,CS,   \&}}
\DoxyCodeLine{208 \textcolor{comment}{!\$OMP                                  Z\_to\_pres,cg1,g\_Rho0,rescale,I\_rescale,L2\_to\_Z2, \&}}
\DoxyCodeLine{209 \textcolor{comment}{!\$OMP                                  better\_est,cg1\_min2,tol\_merge,tol\_solve,c2\_scale) \&}}
\DoxyCodeLine{210 \textcolor{comment}{!\$OMP                          private(htot,hmin,kf,H\_here,HxT\_here,HxS\_here,HxR\_here, \&}}
\DoxyCodeLine{211 \textcolor{comment}{!\$OMP                                  Hf,Tf,Sf,Rf,pres,T\_int,S\_int,drho\_dT,drho\_dS,   \&}}
\DoxyCodeLine{212 \textcolor{comment}{!\$OMP                                  drxh\_sum,kc,Hc,Hc\_H,tC,sc,I\_Hnew,gprime,\&}}
\DoxyCodeLine{213 \textcolor{comment}{!\$OMP                                  Rc,speed2\_tot,Igl,Igu,lam0,lam,lam\_it,dlam, \&}}
\DoxyCodeLine{214 \textcolor{comment}{!\$OMP                                  mode\_struct,sum\_hc,N2min,gp,hw,                 \&}}
\DoxyCodeLine{215 \textcolor{comment}{!\$OMP                                  ms\_min,ms\_max,ms\_sq,H\_top,H\_bot,I\_Htot,merge,   \&}}
\DoxyCodeLine{216 \textcolor{comment}{!\$OMP                                  det,ddet,detKm1,ddetKm1,detKm2,ddetKm2,det\_it,ddet\_it)}}
\DoxyCodeLine{217   \textcolor{keywordflow}{do} j=js,je}
\DoxyCodeLine{218     \textcolor{comment}{!   First merge very thin layers with the one above (or below if they are}}
\DoxyCodeLine{219     \textcolor{comment}{! at the top).  This also transposes the row order so that columns can}}
\DoxyCodeLine{220     \textcolor{comment}{! be worked upon one at a time.}}
\DoxyCodeLine{221     \textcolor{keywordflow}{do} i=is,ie ; htot(i) = 0.0 ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{222     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv\%H\_to\_Z ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{223 }
\DoxyCodeLine{224     \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{225       hmin(i) = htot(i)*min\_h\_frac ; kf(i) = 1 ; h\_here(i) = 0.0}
\DoxyCodeLine{226       hxt\_here(i) = 0.0 ; hxs\_here(i) = 0.0 ; hxr\_here(i) = 0.0}
\DoxyCodeLine{227 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{228     \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{229       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{230         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv\%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}}
\DoxyCodeLine{231           hf(kf(i),i) = h\_here(i)}
\DoxyCodeLine{232           tf(kf(i),i) = hxt\_here(i) / h\_here(i)}
\DoxyCodeLine{233           sf(kf(i),i) = hxs\_here(i) / h\_here(i)}
\DoxyCodeLine{234           kf(i) = kf(i) + 1}
\DoxyCodeLine{235 }
\DoxyCodeLine{236           \textcolor{comment}{! Start a new layer}}
\DoxyCodeLine{237           h\_here(i) = h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{238           hxt\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*t(i,j,k)}
\DoxyCodeLine{239           hxs\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*s(i,j,k)}
\DoxyCodeLine{240         \textcolor{keywordflow}{else}}
\DoxyCodeLine{241           h\_here(i) = h\_here(i) + h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{242           hxt\_here(i) = hxt\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*t(i,j,k)}
\DoxyCodeLine{243           hxs\_here(i) = hxs\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*s(i,j,k)}
\DoxyCodeLine{244 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{245 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{246       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{247         hf(kf(i),i) = h\_here(i)}
\DoxyCodeLine{248         tf(kf(i),i) = hxt\_here(i) / h\_here(i)}
\DoxyCodeLine{249         sf(kf(i),i) = hxs\_here(i) / h\_here(i)}
\DoxyCodeLine{250 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{251     \textcolor{keywordflow}{else}}
\DoxyCodeLine{252       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{253         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv\%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}}
\DoxyCodeLine{254           hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)}
\DoxyCodeLine{255           kf(i) = kf(i) + 1}
\DoxyCodeLine{256 }
\DoxyCodeLine{257           \textcolor{comment}{! Start a new layer}}
\DoxyCodeLine{258           h\_here(i) = h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{259           hxr\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*gv\%Rlay(k)}
\DoxyCodeLine{260         \textcolor{keywordflow}{else}}
\DoxyCodeLine{261           h\_here(i) = h\_here(i) + h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{262           hxr\_here(i) = hxr\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*gv\%Rlay(k)}
\DoxyCodeLine{263 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{264 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{265       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{266         hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)}
\DoxyCodeLine{267 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{268 \textcolor{keywordflow}{    endif}}
\DoxyCodeLine{269 }
\DoxyCodeLine{270     \textcolor{comment}{! From this point, we can work on individual columns without causing memory to have page faults.}}
\DoxyCodeLine{271     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g\%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then}}
\DoxyCodeLine{272       \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{273         pres(1) = 0.0 ; h\_top(1) = 0.0}
\DoxyCodeLine{274         \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{275           pres(k) = pres(k-\/1) + z\_to\_pres*hf(k-\/1,i)}
\DoxyCodeLine{276           t\_int(k) = 0.5*(tf(k,i)+tf(k-\/1,i))}
\DoxyCodeLine{277           s\_int(k) = 0.5*(sf(k,i)+sf(k-\/1,i))}
\DoxyCodeLine{278           h\_top(k) = h\_top(k-\/1) + hf(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 is negligibly small.}}
\DoxyCodeLine{284         drxh\_sum = 0.0}
\DoxyCodeLine{285         \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{286           \textcolor{comment}{! This is an estimate that is correct for the non-\/EBT mode for 2 or 3 layers, or for}}
\DoxyCodeLine{287           \textcolor{comment}{! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.}}
\DoxyCodeLine{288           \textcolor{comment}{! For a uniform stratification and a huge number of layers uniformly distributed in}}
\DoxyCodeLine{289           \textcolor{comment}{! density, this estimate is too large (as is desired) by a factor of pi\string^2/6 ~= 1.64.}}
\DoxyCodeLine{290           \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{291             i\_htot = 1.0 / (h\_top(kf(i)) + hf(kf(i),i))  \textcolor{comment}{! = 1.0 / (H\_top(K) + H\_bot(K)) for all K.}}
\DoxyCodeLine{292             h\_bot(kf(i)+1) = 0.0}
\DoxyCodeLine{293             \textcolor{keywordflow}{do} k=kf(i),2,-\/1}
\DoxyCodeLine{294               h\_bot(k) = h\_bot(k+1) + hf(k,i)}
\DoxyCodeLine{295               drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * \&}
\DoxyCodeLine{296                   max(0.0, drho\_dt(k)*(tf(k,i)-\/tf(k-\/1,i)) + drho\_ds(k)*(sf(k,i)-\/sf(k-\/1,i)))}
\DoxyCodeLine{297 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{298 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{299         \textcolor{keywordflow}{else}}
\DoxyCodeLine{300           \textcolor{comment}{! This estimate is problematic in that it goes like 1/nz for a large number of layers,}}
\DoxyCodeLine{301           \textcolor{comment}{! but it is an overestimate (as desired) for a small number of layers, by at a factor}}
\DoxyCodeLine{302           \textcolor{comment}{! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.}}
\DoxyCodeLine{303           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{304             drxh\_sum = drxh\_sum + 0.5*(hf(k-\/1,i)+hf(k,i)) * \&}
\DoxyCodeLine{305                 max(0.0, drho\_dt(k)*(tf(k,i)-\/tf(k-\/1,i)) + drho\_ds(k)*(sf(k,i)-\/sf(k-\/1,i)))}
\DoxyCodeLine{306 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{307 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{308       \textcolor{keywordflow}{else}}
\DoxyCodeLine{309         drxh\_sum = 0.0}
\DoxyCodeLine{310         \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{311           h\_top(1) = 0.0}
\DoxyCodeLine{312           \textcolor{keywordflow}{do} k=2,kf(i) ; h\_top(k) = h\_top(k-\/1) + hf(k-\/1,i) ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{313           \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{314             i\_htot = 1.0 / (h\_top(kf(i)) + hf(kf(i),i))  \textcolor{comment}{! = 1.0 / (H\_top(K) + H\_bot(K)) for all K.}}
\DoxyCodeLine{315             h\_bot(kf(i)+1) = 0.0}
\DoxyCodeLine{316             \textcolor{keywordflow}{do} k=kf(i),2,-\/1}
\DoxyCodeLine{317               h\_bot(k) = h\_bot(k+1) + hf(k,i)}
\DoxyCodeLine{318               drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * max(0.0,rf(k,i)-\/rf(k-\/1,i))}
\DoxyCodeLine{319 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{320 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{321         \textcolor{keywordflow}{else}}
\DoxyCodeLine{322           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{323             drxh\_sum = drxh\_sum + 0.5*(hf(k-\/1,i)+hf(k,i)) * max(0.0,rf(k,i)-\/rf(k-\/1,i))}
\DoxyCodeLine{324 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{325 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{326 \textcolor{keywordflow}{      endif}}
\DoxyCodeLine{327 }
\DoxyCodeLine{328       \textcolor{comment}{!   Find gprime across each internal interface, taking care of convective instabilities by}}
\DoxyCodeLine{329       \textcolor{comment}{! merging layers.  If the estimated wave speed is too small, simply return zero.}}
\DoxyCodeLine{330       \textcolor{keywordflow}{if} (g\_rho0 * drxh\_sum <= cg1\_min2) \textcolor{keywordflow}{then}}
\DoxyCodeLine{331         cg1(i,j) = 0.0}
\DoxyCodeLine{332         \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) modal\_structure(i,j,:) = 0.}
\DoxyCodeLine{333       \textcolor{keywordflow}{else}}
\DoxyCodeLine{334         \textcolor{comment}{! Merge layers to eliminate convective instabilities or exceedingly}}
\DoxyCodeLine{335         \textcolor{comment}{! small reduced gravities.  Merging layers reduces the estimated wave speed by}}
\DoxyCodeLine{336         \textcolor{comment}{! (rho(2)-\/rho(1))*h(1)*h(2) / H\_tot.}}
\DoxyCodeLine{337         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{338           kc = 1}
\DoxyCodeLine{339           hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)}
\DoxyCodeLine{340           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{341             \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{342               merge = ((drho\_dt(k)*(tf(k,i)-\/tc(kc)) + drho\_ds(k)*(sf(k,i)-\/sc(kc))) * \&}
\DoxyCodeLine{343                        ((hc(kc) * hf(k,i))*i\_htot) < 2.0 * tol\_merge*drxh\_sum)}
\DoxyCodeLine{344             \textcolor{keywordflow}{else}}
\DoxyCodeLine{345               merge = ((drho\_dt(k)*(tf(k,i)-\/tc(kc)) + drho\_ds(k)*(sf(k,i)-\/sc(kc))) * \&}
\DoxyCodeLine{346                        (hc(kc) + hf(k,i)) < 2.0 * tol\_merge*drxh\_sum)}
\DoxyCodeLine{347 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{348             \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{349               \textcolor{comment}{! Merge this layer with the one above and backtrack.}}
\DoxyCodeLine{350               i\_hnew = 1.0 / (hc(kc) + hf(k,i))}
\DoxyCodeLine{351               tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i\_hnew}
\DoxyCodeLine{352               sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i\_hnew}
\DoxyCodeLine{353               hc(kc) = (hc(kc) + hf(k,i))}
\DoxyCodeLine{354               \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}}
\DoxyCodeLine{355               \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}}
\DoxyCodeLine{356               \textcolor{comment}{! far back we go.}}
\DoxyCodeLine{357               \textcolor{keywordflow}{do} k2=kc,2,-\/1}
\DoxyCodeLine{358                 \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{359                   merge = ((drho\_dt(k2)*(tc(k2)-\/tc(k2-\/1)) + drho\_ds(k2)*(sc(k2)-\/sc(k2-\/1))) * \&}
\DoxyCodeLine{360                            ((hc(k2) * hc(k2-\/1))*i\_htot) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{361                 \textcolor{keywordflow}{else}}
\DoxyCodeLine{362                   merge = ((drho\_dt(k2)*(tc(k2)-\/tc(k2-\/1)) + drho\_ds(k2)*(sc(k2)-\/sc(k2-\/1))) * \&}
\DoxyCodeLine{363                            (hc(k2) + hc(k2-\/1)) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{364 \textcolor{keywordflow}{                endif}}
\DoxyCodeLine{365                 \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{366                   \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}}
\DoxyCodeLine{367                   i\_hnew = 1.0 / (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{368                   tc(kc-\/1) = (hc(kc)*tc(kc) + hc(kc-\/1)*tc(kc-\/1)) * i\_hnew}
\DoxyCodeLine{369                   sc(kc-\/1) = (hc(kc)*sc(kc) + hc(kc-\/1)*sc(kc-\/1)) * i\_hnew}
\DoxyCodeLine{370                   hc(kc-\/1) = (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{371                   kc = kc -\/ 1}
\DoxyCodeLine{372                 \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{373 \textcolor{keywordflow}{              enddo}}
\DoxyCodeLine{374             \textcolor{keywordflow}{else}}
\DoxyCodeLine{375               \textcolor{comment}{! Add a new layer to the column.}}
\DoxyCodeLine{376               kc = kc + 1}
\DoxyCodeLine{377               drho\_ds(kc) = drho\_ds(k) ; drho\_dt(kc) = drho\_dt(k)}
\DoxyCodeLine{378               tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)}
\DoxyCodeLine{379 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{380 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{381           \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}}
\DoxyCodeLine{382           \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-\/Boussinesq.}}
\DoxyCodeLine{383             gprime(k) = g\_rho0 * (drho\_dt(k)*(tc(k)-\/tc(k-\/1)) + drho\_ds(k)*(sc(k)-\/sc(k-\/1)))}
\DoxyCodeLine{384 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{385         \textcolor{keywordflow}{else}  \textcolor{comment}{! .not.use\_EOS}}
\DoxyCodeLine{386           \textcolor{comment}{! Do the same with density directly...}}
\DoxyCodeLine{387           kc = 1}
\DoxyCodeLine{388           hc(1) = hf(1,i) ; rc(1) = rf(1,i)}
\DoxyCodeLine{389           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{390             \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{391               merge = ((rf(k,i) -\/ rc(kc)) * ((hc(kc) * hf(k,i))*i\_htot) < 2.0*tol\_merge*drxh\_sum)}
\DoxyCodeLine{392             \textcolor{keywordflow}{else}}
\DoxyCodeLine{393               merge = ((rf(k,i) -\/ rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol\_merge*drxh\_sum)}
\DoxyCodeLine{394 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{395             \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{396               \textcolor{comment}{! Merge this layer with the one above and backtrack.}}
\DoxyCodeLine{397               rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))}
\DoxyCodeLine{398               hc(kc) = (hc(kc) + hf(k,i))}
\DoxyCodeLine{399               \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}}
\DoxyCodeLine{400               \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}}
\DoxyCodeLine{401               \textcolor{comment}{! far back we go.}}
\DoxyCodeLine{402               \textcolor{keywordflow}{do} k2=kc,2,-\/1}
\DoxyCodeLine{403                 \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{404                   merge = ((rc(k2)-\/rc(k2-\/1)) * ((hc(k2) * hc(k2-\/1))*i\_htot) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{405                 \textcolor{keywordflow}{else}}
\DoxyCodeLine{406                   merge = ((rc(k2)-\/rc(k2-\/1)) * (hc(k2)+hc(k2-\/1)) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{407 \textcolor{keywordflow}{                endif}}
\DoxyCodeLine{408                 \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{409                   \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}}
\DoxyCodeLine{410                   rc(kc-\/1) = (hc(kc)*rc(kc) + hc(kc-\/1)*rc(kc-\/1)) / (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{411                   hc(kc-\/1) = (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{412                   kc = kc -\/ 1}
\DoxyCodeLine{413                 \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{414 \textcolor{keywordflow}{              enddo}}
\DoxyCodeLine{415             \textcolor{keywordflow}{else}}
\DoxyCodeLine{416               \textcolor{comment}{! Add a new layer to the column.}}
\DoxyCodeLine{417               kc = kc + 1}
\DoxyCodeLine{418               rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)}
\DoxyCodeLine{419 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{420 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{421           \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}}
\DoxyCodeLine{422           \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-\/Boussinesq.}}
\DoxyCodeLine{423             gprime(k) = g\_rho0 * (rc(k)-\/rc(k-\/1))}
\DoxyCodeLine{424 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{425 \textcolor{keywordflow}{        endif}  \textcolor{comment}{! use\_EOS}}
\DoxyCodeLine{426 }
\DoxyCodeLine{427         \textcolor{comment}{! Sum the contributions from all of the interfaces to give an over-\/estimate}}
\DoxyCodeLine{428         \textcolor{comment}{! of the first-\/mode wave speed.  Also populate Igl and Igu which are the}}
\DoxyCodeLine{429         \textcolor{comment}{! non-\/leading diagonals of the tridiagonal matrix.}}
\DoxyCodeLine{430         \textcolor{keywordflow}{if} (kc >= 2) \textcolor{keywordflow}{then}}
\DoxyCodeLine{431           speed2\_tot = 0.0}
\DoxyCodeLine{432           \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{433             h\_top(1) = 0.0 ; h\_bot(kc+1) = 0.0}
\DoxyCodeLine{434             \textcolor{keywordflow}{do} k=2,kc+1 ; h\_top(k) = h\_top(k-\/1) + hc(k-\/1) ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{435             \textcolor{keywordflow}{do} k=kc,2,-\/1 ; h\_bot(k) = h\_bot(k+1) + hc(k) ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{436             i\_htot = 0.0 ; \textcolor{keywordflow}{if} (h\_top(kc+1) > 0.0) i\_htot = 1.0 / h\_top(kc+1)}
\DoxyCodeLine{437 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{438 }
\DoxyCodeLine{439           \textcolor{keywordflow}{if} (l\_use\_ebt\_mode) \textcolor{keywordflow}{then}}
\DoxyCodeLine{440             igu(1) = 0. \textcolor{comment}{! Neumann condition for pressure modes}}
\DoxyCodeLine{441             sum\_hc = hc(1)}
\DoxyCodeLine{442             n2min = l2\_to\_z2*gprime(2)/hc(1)}
\DoxyCodeLine{443             \textcolor{keywordflow}{do} k=2,kc}
\DoxyCodeLine{444               hw = 0.5*(hc(k-\/1)+hc(k))}
\DoxyCodeLine{445               gp = gprime(k)}
\DoxyCodeLine{446               \textcolor{keywordflow}{if} (l\_mono\_n2\_column\_fraction>0. .or. l\_mono\_n2\_depth>=0.) \textcolor{keywordflow}{then}}
\DoxyCodeLine{447                 \textcolor{keywordflow}{if} ( ((g\%bathyT(i,j)-\/sum\_hc < l\_mono\_n2\_column\_fraction*g\%bathyT(i,j)) .or. \&}
\DoxyCodeLine{448                       ((l\_mono\_n2\_depth >= 0.) .and. (sum\_hc > l\_mono\_n2\_depth))) .and. \&}
\DoxyCodeLine{449                      (l2\_to\_z2*gp > n2min*hw) ) \textcolor{keywordflow}{then}}
\DoxyCodeLine{450                   \textcolor{comment}{! Filters out regions where N2 increases with depth but only in a lower fraction}}
\DoxyCodeLine{451                   \textcolor{comment}{! of the water column or below a certain depth.}}
\DoxyCodeLine{452                   gp = us\%Z\_to\_L**2 * (n2min*hw)}
\DoxyCodeLine{453                 \textcolor{keywordflow}{else}}
\DoxyCodeLine{454                   n2min = l2\_to\_z2 * gp/hw}
\DoxyCodeLine{455 \textcolor{keywordflow}{                endif}}
\DoxyCodeLine{456 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{457               igu(k) = 1.0/(gp*hc(k))}
\DoxyCodeLine{458               igl(k-\/1) = 1.0/(gp*hc(k-\/1))}
\DoxyCodeLine{459               sum\_hc = sum\_hc + hc(k)}
\DoxyCodeLine{460               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{461                 \textcolor{comment}{! Estimate that the ebt\_mode is sqrt(2) times the speed of the flat bottom modes.}}
\DoxyCodeLine{462                 speed2\_tot = speed2\_tot + 2.0 * gprime(k)*((h\_top(k) * h\_bot(k)) * i\_htot)}
\DoxyCodeLine{463               \textcolor{keywordflow}{else} \textcolor{comment}{! The ebt\_mode wave should be faster than the flat-\/bottom mode, so 0.707 should be > 1?}}
\DoxyCodeLine{464                 speed2\_tot = speed2\_tot + gprime(k)*(hc(k-\/1)+hc(k))*0.707}
\DoxyCodeLine{465 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{466 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{467            \textcolor{comment}{!Igl(kc) = 0. ! Neumann condition for pressure modes}}
\DoxyCodeLine{468             igl(kc) = 2.*igu(kc) \textcolor{comment}{! Dirichlet condition for pressure modes}}
\DoxyCodeLine{469           \textcolor{keywordflow}{else} \textcolor{comment}{! .not. l\_use\_ebt\_mode}}
\DoxyCodeLine{470             \textcolor{keywordflow}{do} k=2,kc}
\DoxyCodeLine{471               igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-\/1))}
\DoxyCodeLine{472               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{473                 speed2\_tot = speed2\_tot + gprime(k)*((h\_top(k) * h\_bot(k)) * i\_htot)}
\DoxyCodeLine{474               \textcolor{keywordflow}{else}}
\DoxyCodeLine{475                 speed2\_tot = speed2\_tot + gprime(k)*(hc(k-\/1)+hc(k))}
\DoxyCodeLine{476 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{477 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{478 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{479 }
\DoxyCodeLine{480           \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}}
\DoxyCodeLine{481             mode\_struct(:) = 0.}
\DoxyCodeLine{482             mode\_struct(1:kc) = 1. \textcolor{comment}{! Uniform flow, first guess}}
\DoxyCodeLine{483 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{484 }
\DoxyCodeLine{485           \textcolor{comment}{! Under estimate the first eigenvalue (overestimate the speed) to start with.}}
\DoxyCodeLine{486           \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}}
\DoxyCodeLine{487             lam0 = 0.5 / speed2\_tot ; lam = lam0}
\DoxyCodeLine{488           \textcolor{keywordflow}{else}}
\DoxyCodeLine{489             lam0 = 1.0 / speed2\_tot ; lam = lam0}
\DoxyCodeLine{490 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{491           \textcolor{comment}{! Find the determinant and its derivative with lam.}}
\DoxyCodeLine{492           \textcolor{keywordflow}{do} itt=1,max\_itt}
\DoxyCodeLine{493             lam\_it(itt) = lam}
\DoxyCodeLine{494             \textcolor{keywordflow}{if} (l\_use\_ebt\_mode) \textcolor{keywordflow}{then}}
\DoxyCodeLine{495               \textcolor{comment}{! This initialization of det,ddet imply Neumann boundary conditions for horizontal}}
\DoxyCodeLine{496               \textcolor{comment}{! velocity or pressure modes, so that first 3 rows of the matrix are}}
\DoxyCodeLine{497               \textcolor{comment}{!    /   b(1)-\/lam  igl(1)      0        0     0  ...  \(\backslash\)}}
\DoxyCodeLine{498               \textcolor{comment}{!    |  igu(2)    b(2)-\/lam   igl(2)     0     0  ...  |}}
\DoxyCodeLine{499               \textcolor{comment}{!    |    0        igu(3)   b(3)-\/lam  igl(3)  0  ...  |}}
\DoxyCodeLine{500               \textcolor{comment}{! The last two rows of the pressure equation matrix are}}
\DoxyCodeLine{501               \textcolor{comment}{!    |    ...  0  igu(kc-\/1)  b(kc-\/1)-\/lam  igl(kc-\/1)  |}}
\DoxyCodeLine{502               \textcolor{comment}{!    \(\backslash\)    ...  0     0        igu(kc)     b(kc)-\/lam  /}}
\DoxyCodeLine{503               \textcolor{keyword}{call }tridiag\_det(igu, igl, 1, kc, lam, det, ddet, row\_scale=c2\_scale)}
\DoxyCodeLine{504             \textcolor{keywordflow}{else}}
\DoxyCodeLine{505               \textcolor{comment}{! This initialization of det,ddet imply Dirichlet boundary conditions for vertical}}
\DoxyCodeLine{506               \textcolor{comment}{! velocity modes, so that first 3 rows of the matrix are}}
\DoxyCodeLine{507               \textcolor{comment}{!    /  b(2)-\/lam  igl(2)      0       0     0  ...  |}}
\DoxyCodeLine{508               \textcolor{comment}{!    |  igu(3)  b(3)-\/lam   igl(3)     0     0  ...  |}}
\DoxyCodeLine{509               \textcolor{comment}{!    |    0       igu(4)  b(4)-\/lam  igl(4)  0  ...  |}}
\DoxyCodeLine{510               \textcolor{comment}{! The last three rows of the w equation matrix are}}
\DoxyCodeLine{511               \textcolor{comment}{!    |    ...   0  igu(kc-\/2)  b(kc-\/2)-\/lam  igl(kc-\/2)     0       |}}
\DoxyCodeLine{512               \textcolor{comment}{!    |    ...   0     0        igu(kc-\/1)  b(kc-\/1)-\/lam  igl(kc-\/1) |}}
\DoxyCodeLine{513               \textcolor{comment}{!    \(\backslash\)    ...   0     0           0        igu(kc)    b(kc)-\/lam  /}}
\DoxyCodeLine{514               \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lam, det, ddet, row\_scale=c2\_scale)}
\DoxyCodeLine{515 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{516             \textcolor{comment}{! Use Newton's method iteration to find a new estimate of lam.}}
\DoxyCodeLine{517             det\_it(itt) = det ; ddet\_it(itt) = ddet}
\DoxyCodeLine{518 }
\DoxyCodeLine{519             \textcolor{keywordflow}{if} ((ddet >= 0.0) .or. (-\/det > -\/0.5*lam*ddet)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{520               \textcolor{comment}{! lam was not an under-\/estimate, as intended, so Newton's method}}
\DoxyCodeLine{521               \textcolor{comment}{! may not be reliable; lam must be reduced, but not by more}}
\DoxyCodeLine{522               \textcolor{comment}{! than half.}}
\DoxyCodeLine{523               lam = 0.5 * lam}
\DoxyCodeLine{524               dlam = -\/lam}
\DoxyCodeLine{525             \textcolor{keywordflow}{else}  \textcolor{comment}{! Newton's method is OK.}}
\DoxyCodeLine{526               dlam = -\/ det / ddet}
\DoxyCodeLine{527               lam = lam + dlam}
\DoxyCodeLine{528 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{529 }
\DoxyCodeLine{530             \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}}
\DoxyCodeLine{531               \textcolor{keyword}{call }tdma6(kc, igu, igl, lam, mode\_struct)}
\DoxyCodeLine{532               ms\_min = mode\_struct(1)}
\DoxyCodeLine{533               ms\_max = mode\_struct(1)}
\DoxyCodeLine{534               ms\_sq = mode\_struct(1)**2}
\DoxyCodeLine{535               \textcolor{keywordflow}{do} k = 2,kc}
\DoxyCodeLine{536                 ms\_min = min(ms\_min, mode\_struct(k))}
\DoxyCodeLine{537                 ms\_max = max(ms\_max, mode\_struct(k))}
\DoxyCodeLine{538                 ms\_sq = ms\_sq + mode\_struct(k)**2}
\DoxyCodeLine{539 \textcolor{keywordflow}{              enddo}}
\DoxyCodeLine{540               \textcolor{keywordflow}{if} (ms\_min<0. .and. ms\_max>0.) \textcolor{keywordflow}{then} \textcolor{comment}{! Any zero crossings => lam is too high}}
\DoxyCodeLine{541                 lam = 0.5 * ( lam -\/ dlam )}
\DoxyCodeLine{542                 dlam = -\/lam}
\DoxyCodeLine{543                 mode\_struct(1:kc) = abs(mode\_struct(1:kc)) / sqrt( ms\_sq )}
\DoxyCodeLine{544               \textcolor{keywordflow}{else}}
\DoxyCodeLine{545                 mode\_struct(1:kc) = mode\_struct(1:kc) / sqrt( ms\_sq )}
\DoxyCodeLine{546 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{547 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{548 }
\DoxyCodeLine{549             \textcolor{keywordflow}{if} (abs(dlam) < tol\_solve*lam) \textcolor{keywordflow}{exit}}
\DoxyCodeLine{550 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{551 }
\DoxyCodeLine{552           cg1(i,j) = 0.0}
\DoxyCodeLine{553           \textcolor{keywordflow}{if} (lam > 0.0) cg1(i,j) = 1.0 / sqrt(lam)}
\DoxyCodeLine{554 }
\DoxyCodeLine{555           \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{556             \textcolor{keywordflow}{if} (mode\_struct(1)/=0.) \textcolor{keywordflow}{then} \textcolor{comment}{! Normalize}}
\DoxyCodeLine{557               mode\_struct(1:kc) = mode\_struct(1:kc) / mode\_struct(1)}
\DoxyCodeLine{558             \textcolor{keywordflow}{else}}
\DoxyCodeLine{559               mode\_struct(1:kc)=0.}
\DoxyCodeLine{560 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{561             \textcolor{comment}{! Note that remapping\_core\_h requires that the same units be used}}
\DoxyCodeLine{562             \textcolor{comment}{! for both the source and target grid thicknesses, here [H ~> m or kg m-\/2].}}
\DoxyCodeLine{563             \textcolor{keywordflow}{do} k = 1,kc}
\DoxyCodeLine{564               hc\_h(k) = gv\%Z\_to\_H * hc(k)}
\DoxyCodeLine{565 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{566             \textcolor{keywordflow}{if} (cs\%remap\_answers\_2018) \textcolor{keywordflow}{then}}
\DoxyCodeLine{567               \textcolor{keyword}{call }remapping\_core\_h(cs\%remapping\_CS, kc, hc\_h(:), mode\_struct, \&}
\DoxyCodeLine{568                                     nz, h(i,j,:), modal\_structure(i,j,:), \&}
\DoxyCodeLine{569                                     1.0e-\/30*gv\%m\_to\_H, 1.0e-\/10*gv\%m\_to\_H)}
\DoxyCodeLine{570             \textcolor{keywordflow}{else}}
\DoxyCodeLine{571               \textcolor{keyword}{call }remapping\_core\_h(cs\%remapping\_CS, kc, hc\_h(:), mode\_struct, \&}
\DoxyCodeLine{572                                     nz, h(i,j,:), modal\_structure(i,j,:), \&}
\DoxyCodeLine{573                                     gv\%H\_subroundoff, gv\%H\_subroundoff)}
\DoxyCodeLine{574 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{575 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{576         \textcolor{keywordflow}{else}}
\DoxyCodeLine{577           cg1(i,j) = 0.0}
\DoxyCodeLine{578           \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) modal\_structure(i,j,:) = 0.}
\DoxyCodeLine{579 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{580 \textcolor{keywordflow}{      endif} \textcolor{comment}{! cg1 /= 0.0}}
\DoxyCodeLine{581     \textcolor{keywordflow}{else}}
\DoxyCodeLine{582       cg1(i,j) = 0.0 \textcolor{comment}{! This is a land point.}}
\DoxyCodeLine{583       \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) modal\_structure(i,j,:) = 0.}
\DoxyCodeLine{584 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i-\/loop}}
\DoxyCodeLine{585 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-\/loop}}
\DoxyCodeLine{586 }

\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__wave__speed_afa7284f32921f1e5d31530633cf95022}\label{namespacemom__wave__speed_afa7284f32921f1e5d31530633cf95022}} 
\index{mom\_wave\_speed@{mom\_wave\_speed}!wave\_speed\_init@{wave\_speed\_init}}
\index{wave\_speed\_init@{wave\_speed\_init}!mom\_wave\_speed@{mom\_wave\_speed}}
\doxysubsubsection{\texorpdfstring{wave\_speed\_init()}{wave\_speed\_init()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+speed\+::wave\+\_\+speed\+\_\+init (\begin{DoxyParamCaption}\item[{type(\mbox{\hyperlink{structmom__wave__speed_1_1wave__speed__cs}{wave\+\_\+speed\+\_\+cs}}), pointer}]{CS,  }\item[{logical, intent(in), optional}]{use\+\_\+ebt\+\_\+mode,  }\item[{real, intent(in), optional}]{mono\+\_\+\+N2\+\_\+column\+\_\+fraction,  }\item[{real, intent(in), optional}]{mono\+\_\+\+N2\+\_\+depth,  }\item[{logical, intent(in), optional}]{remap\+\_\+answers\+\_\+2018,  }\item[{logical, intent(in), optional}]{better\+\_\+speed\+\_\+est,  }\item[{real, intent(in), optional}]{min\+\_\+speed,  }\item[{real, intent(in), optional}]{wave\+\_\+speed\+\_\+tol }\end{DoxyParamCaption})}



Initialize control structure for M\+O\+M\+\_\+wave\+\_\+speed. 


\begin{DoxyParams}[1]{Parameters}
 & {\em cs} & Control structure for M\+O\+M\+\_\+wave\+\_\+speed \\
\hline
\mbox{\texttt{ in}}  & {\em use\+\_\+ebt\+\_\+mode} & If true, use the equivalent barotropic mode instead of the first baroclinic mode. \\
\hline
\mbox{\texttt{ in}}  & {\em mono\+\_\+n2\+\_\+column\+\_\+fraction} & The lower fraction of water column over which N2 is limited as monotonic for the purposes of calculating the vertical modal structure. \\
\hline
\mbox{\texttt{ in}}  & {\em mono\+\_\+n2\+\_\+depth} & The depth below which N2 is limited as monotonic for the purposes of calculating the vertical modal structure \mbox{[}Z $\sim$$>$ m\mbox{]}. \\
\hline
\mbox{\texttt{ in}}  & {\em remap\+\_\+answers\+\_\+2018} & If true, use the order of arithmetic and expressions that recover the remapping answers from 2018. Otherwise use more robust but mathematically equivalent expressions. \\
\hline
\mbox{\texttt{ in}}  & {\em better\+\_\+speed\+\_\+est} & If true, use a more robust estimate of the first mode speed as the starting point for iterations. \\
\hline
\mbox{\texttt{ in}}  & {\em min\+\_\+speed} & If present, set a floor in the first mode speed below which 0 is returned \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}. \\
\hline
\mbox{\texttt{ in}}  & {\em wave\+\_\+speed\+\_\+tol} & The fractional tolerance for finding the wave speeds \mbox{[}nondim\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 1191 of file M\+O\+M\+\_\+wave\+\_\+speed.\+F90.


\begin{DoxyCode}{0}
\DoxyCodeLine{1191   \textcolor{keywordtype}{type}(wave\_speed\_CS), \textcolor{keywordtype}{pointer} :: CS\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}}
\DoxyCodeLine{1192   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: use\_ebt\_mode\textcolor{comment}{  !< If true, use the equivalent}}
\DoxyCodeLine{1193 \textcolor{comment}{                                     !! barotropic mode instead of the first baroclinic mode.}}
\DoxyCodeLine{1194 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_N2\_column\_fraction\textcolor{comment}{ !< The lower fraction of water column over}}
\DoxyCodeLine{1195 \textcolor{comment}{                                     !! which N2 is limited as monotonic for the purposes of}}
\DoxyCodeLine{1196 \textcolor{comment}{                                     !! calculating the vertical modal structure.}}
\DoxyCodeLine{1197 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_N2\_depth\textcolor{comment}{ !< The depth below which N2 is limited}}
\DoxyCodeLine{1198 \textcolor{comment}{                                     !! as monotonic for the purposes of calculating the}}
\DoxyCodeLine{1199 \textcolor{comment}{                                     !! vertical modal structure [Z ~> m].}}
\DoxyCodeLine{1200   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: remap\_answers\_2018\textcolor{comment}{ !< If true, use the order of arithmetic and expressions}}
\DoxyCodeLine{1201 \textcolor{comment}{                                     !! that recover the remapping answers from 2018.  Otherwise}}
\DoxyCodeLine{1202 \textcolor{comment}{                                     !! use more robust but mathematically equivalent expressions.}}
\DoxyCodeLine{1203   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: better\_speed\_est\textcolor{comment}{ !< If true, use a more robust estimate of the first}}
\DoxyCodeLine{1204 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}}
\DoxyCodeLine{1205 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: min\_speed\textcolor{comment}{ !< If present, set a floor in the first mode speed}}
\DoxyCodeLine{1206 \textcolor{comment}{                                     !! below which 0 is returned [L T-\/1 ~> m s-\/1].}}
\DoxyCodeLine{1207 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}}
\DoxyCodeLine{1208 \textcolor{comment}{                                     !! wave speeds [nondim]}}
\DoxyCodeLine{1209 }
\DoxyCodeLine{1210   \textcolor{comment}{! This include declares and sets the variable "version".}}
\DoxyCodeLine{1211 \textcolor{preprocessor}{\# include "version\_variable.h"}}
\DoxyCodeLine{1212 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_wave\_speed"}  \textcolor{comment}{! This module's name.}}
\DoxyCodeLine{1213 }
\DoxyCodeLine{1214   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1215     \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_speed\_init called with an "}// \&}
\DoxyCodeLine{1216                             \textcolor{stringliteral}{"associated control structure."})}
\DoxyCodeLine{1217     \textcolor{keywordflow}{return}}
\DoxyCodeLine{1218   \textcolor{keywordflow}{else} ; \textcolor{keyword}{allocate}(cs) ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{1219 }
\DoxyCodeLine{1220   \textcolor{comment}{! Write all relevant parameters to the model log.}}
\DoxyCodeLine{1221   \textcolor{keyword}{call }log\_version(mdl, version)}
\DoxyCodeLine{1222 }
\DoxyCodeLine{1223   \textcolor{keyword}{call }wave\_speed\_set\_param(cs, use\_ebt\_mode=use\_ebt\_mode, mono\_n2\_column\_fraction=mono\_n2\_column\_fraction, \&}
\DoxyCodeLine{1224                             better\_speed\_est=better\_speed\_est, min\_speed=min\_speed, wave\_speed\_tol=wave\_speed\_tol)}
\DoxyCodeLine{1225 }
\DoxyCodeLine{1226   \textcolor{keyword}{call }initialize\_remapping(cs\%remapping\_CS, \textcolor{stringliteral}{'PLM'}, boundary\_extrapolation=.false., \&}
\DoxyCodeLine{1227                             answers\_2018=cs\%remap\_answers\_2018)}
\DoxyCodeLine{1228 }

\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__wave__speed_a13f3b425b43466d5af43f15d26c25a59}\label{namespacemom__wave__speed_a13f3b425b43466d5af43f15d26c25a59}} 
\index{mom\_wave\_speed@{mom\_wave\_speed}!wave\_speed\_set\_param@{wave\_speed\_set\_param}}
\index{wave\_speed\_set\_param@{wave\_speed\_set\_param}!mom\_wave\_speed@{mom\_wave\_speed}}
\doxysubsubsection{\texorpdfstring{wave\_speed\_set\_param()}{wave\_speed\_set\_param()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+speed\+::wave\+\_\+speed\+\_\+set\+\_\+param (\begin{DoxyParamCaption}\item[{type(\mbox{\hyperlink{structmom__wave__speed_1_1wave__speed__cs}{wave\+\_\+speed\+\_\+cs}}), pointer}]{CS,  }\item[{logical, intent(in), optional}]{use\+\_\+ebt\+\_\+mode,  }\item[{real, intent(in), optional}]{mono\+\_\+\+N2\+\_\+column\+\_\+fraction,  }\item[{real, intent(in), optional}]{mono\+\_\+\+N2\+\_\+depth,  }\item[{logical, intent(in), optional}]{remap\+\_\+answers\+\_\+2018,  }\item[{logical, intent(in), optional}]{better\+\_\+speed\+\_\+est,  }\item[{real, intent(in), optional}]{min\+\_\+speed,  }\item[{real, intent(in), optional}]{wave\+\_\+speed\+\_\+tol }\end{DoxyParamCaption})}



Sets internal parameters for M\+O\+M\+\_\+wave\+\_\+speed. 


\begin{DoxyParams}[1]{Parameters}
 & {\em cs} & Control structure for M\+O\+M\+\_\+wave\+\_\+speed \\
\hline
\mbox{\texttt{ in}}  & {\em use\+\_\+ebt\+\_\+mode} & If true, use the equivalent barotropic mode instead of the first baroclinic mode. \\
\hline
\mbox{\texttt{ in}}  & {\em mono\+\_\+n2\+\_\+column\+\_\+fraction} & The lower fraction of water column over which N2 is limited as monotonic for the purposes of calculating the vertical modal structure. \\
\hline
\mbox{\texttt{ in}}  & {\em mono\+\_\+n2\+\_\+depth} & The depth below which N2 is limited as monotonic for the purposes of calculating the vertical modal structure \mbox{[}Z $\sim$$>$ m\mbox{]}. \\
\hline
\mbox{\texttt{ in}}  & {\em remap\+\_\+answers\+\_\+2018} & If true, use the order of arithmetic and expressions that recover the remapping answers from 2018. Otherwise use more robust but mathematically equivalent expressions. \\
\hline
\mbox{\texttt{ in}}  & {\em better\+\_\+speed\+\_\+est} & If true, use a more robust estimate of the first mode speed as the starting point for iterations. \\
\hline
\mbox{\texttt{ in}}  & {\em min\+\_\+speed} & If present, set a floor in the first mode speed below which 0 is returned \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}. \\
\hline
\mbox{\texttt{ in}}  & {\em wave\+\_\+speed\+\_\+tol} & The fractional tolerance for finding the wave speeds \mbox{[}nondim\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 1234 of file M\+O\+M\+\_\+wave\+\_\+speed.\+F90.


\begin{DoxyCode}{0}
\DoxyCodeLine{1234   \textcolor{keywordtype}{type}(wave\_speed\_CS), \textcolor{keywordtype}{pointer}  :: CS\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}}
\DoxyCodeLine{1235   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: use\_ebt\_mode\textcolor{comment}{  !< If true, use the equivalent}}
\DoxyCodeLine{1236 \textcolor{comment}{                                      !! barotropic mode instead of the first baroclinic mode.}}
\DoxyCodeLine{1237 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_N2\_column\_fraction\textcolor{comment}{ !< The lower fraction of water column over}}
\DoxyCodeLine{1238 \textcolor{comment}{                                      !! which N2 is limited as monotonic for the purposes of}}
\DoxyCodeLine{1239 \textcolor{comment}{                                      !! calculating the vertical modal structure.}}
\DoxyCodeLine{1240 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_N2\_depth\textcolor{comment}{ !< The depth below which N2 is limited}}
\DoxyCodeLine{1241 \textcolor{comment}{                                      !! as monotonic for the purposes of calculating the}}
\DoxyCodeLine{1242 \textcolor{comment}{                                      !! vertical modal structure [Z ~> m].}}
\DoxyCodeLine{1243   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: remap\_answers\_2018\textcolor{comment}{ !< If true, use the order of arithmetic and expressions}}
\DoxyCodeLine{1244 \textcolor{comment}{                                      !! that recover the remapping answers from 2018.  Otherwise}}
\DoxyCodeLine{1245 \textcolor{comment}{                                      !! use more robust but mathematically equivalent expressions.}}
\DoxyCodeLine{1246   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: better\_speed\_est\textcolor{comment}{ !< If true, use a more robust estimate of the first}}
\DoxyCodeLine{1247 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}}
\DoxyCodeLine{1248 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: min\_speed\textcolor{comment}{ !< If present, set a floor in the first mode speed}}
\DoxyCodeLine{1249 \textcolor{comment}{                                     !! below which 0 is returned [L T-\/1 ~> m s-\/1].}}
\DoxyCodeLine{1250 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}}
\DoxyCodeLine{1251 \textcolor{comment}{                                     !! wave speeds [nondim]}}
\DoxyCodeLine{1252 }
\DoxyCodeLine{1253   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \&}
\DoxyCodeLine{1254      \textcolor{stringliteral}{"wave\_speed\_set\_param called with an associated control structure."})}
\DoxyCodeLine{1255 }
\DoxyCodeLine{1256   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(use\_ebt\_mode)) cs\%use\_ebt\_mode = use\_ebt\_mode}
\DoxyCodeLine{1257   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_column\_fraction)) cs\%mono\_N2\_column\_fraction = mono\_n2\_column\_fraction}
\DoxyCodeLine{1258   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_depth)) cs\%mono\_N2\_depth = mono\_n2\_depth}
\DoxyCodeLine{1259   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(remap\_answers\_2018)) cs\%remap\_answers\_2018 = remap\_answers\_2018}
\DoxyCodeLine{1260   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(better\_speed\_est)) cs\%better\_cg1\_est = better\_speed\_est}
\DoxyCodeLine{1261   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(min\_speed)) cs\%min\_speed2 = min\_speed**2}
\DoxyCodeLine{1262   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(wave\_speed\_tol)) cs\%wave\_speed\_tol = wave\_speed\_tol}
\DoxyCodeLine{1263 }

\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__wave__speed_a4e969f10d0cdd765f6804973a74ff479}\label{namespacemom__wave__speed_a4e969f10d0cdd765f6804973a74ff479}} 
\index{mom\_wave\_speed@{mom\_wave\_speed}!wave\_speeds@{wave\_speeds}}
\index{wave\_speeds@{wave\_speeds}!mom\_wave\_speed@{mom\_wave\_speed}}
\doxysubsubsection{\texorpdfstring{wave\_speeds()}{wave\_speeds()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+speed\+::wave\+\_\+speeds (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), 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[{integer, intent(in)}]{nmodes,  }\item[{real, dimension(g\%isd\+:g\%ied,g\%jsd\+:g\%jed,nmodes), intent(out)}]{cn,  }\item[{type(\mbox{\hyperlink{structmom__wave__speed_1_1wave__speed__cs}{wave\+\_\+speed\+\_\+cs}}), optional, pointer}]{CS,  }\item[{logical, intent(in), optional}]{full\+\_\+halos,  }\item[{logical, intent(in), optional}]{better\+\_\+speed\+\_\+est,  }\item[{real, intent(in), optional}]{min\+\_\+speed,  }\item[{real, intent(in), optional}]{wave\+\_\+speed\+\_\+tol }\end{DoxyParamCaption})}



Calculates the wave speeds for the first few barolinic modes. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\texttt{ in}}  & {\em g} & Ocean grid structure \\
\hline
\mbox{\texttt{ in}}  & {\em gv} & Vertical grid structure \\
\hline
\mbox{\texttt{ in}}  & {\em us} & A dimensional unit scaling type \\
\hline
\mbox{\texttt{ in}}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]} \\
\hline
\mbox{\texttt{ in}}  & {\em tv} & Thermodynamic variables \\
\hline
\mbox{\texttt{ in}}  & {\em nmodes} & Number of modes \\
\hline
\mbox{\texttt{ out}}  & {\em cn} & Waves speeds \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]} \\
\hline
 & {\em cs} & Control structure for M\+O\+M\+\_\+wave\+\_\+speed \\
\hline
\mbox{\texttt{ in}}  & {\em full\+\_\+halos} & If true, do the calculation over the entire computational domain. \\
\hline
\mbox{\texttt{ in}}  & {\em better\+\_\+speed\+\_\+est} & If true, use a more robust estimate of the first mode speed as the starting point for iterations. \\
\hline
\mbox{\texttt{ in}}  & {\em min\+\_\+speed} & If present, set a floor in the first mode speed below which 0 is returned \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}. \\
\hline
\mbox{\texttt{ in}}  & {\em wave\+\_\+speed\+\_\+tol} & The fractional tolerance for finding the wave speeds \mbox{[}nondim\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 642 of file M\+O\+M\+\_\+wave\+\_\+speed.\+F90.


\begin{DoxyCode}{0}
\DoxyCodeLine{642   \textcolor{keywordtype}{type}(ocean\_grid\_type),                    \textcolor{keywordtype}{intent(in)}  :: G\textcolor{comment}{ !< Ocean grid structure}}
\DoxyCodeLine{643   \textcolor{keywordtype}{type}(verticalGrid\_type),                  \textcolor{keywordtype}{intent(in)}  :: GV\textcolor{comment}{ !< Vertical grid structure}}
\DoxyCodeLine{644   \textcolor{keywordtype}{type}(unit\_scale\_type),                    \textcolor{keywordtype}{intent(in)}  :: US\textcolor{comment}{ !< A dimensional unit scaling type}}
\DoxyCodeLine{645 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{ !< Layer thickness [H ~> m or kg m-\/2]}}
\DoxyCodeLine{646   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),                    \textcolor{keywordtype}{intent(in)}  :: tv\textcolor{comment}{ !< Thermodynamic variables}}
\DoxyCodeLine{647   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}  :: nmodes\textcolor{comment}{ !< Number of modes}}
\DoxyCodeLine{648 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(G\%isd:G\%ied,G\%jsd:G\%jed,nmodes)}, \textcolor{keywordtype}{intent(out)} :: cn\textcolor{comment}{ !< Waves speeds [L T-\/1 ~> m s-\/1]}}
\DoxyCodeLine{649   \textcolor{keywordtype}{type}(wave\_speed\_CS), \textcolor{keywordtype}{optional},            \textcolor{keywordtype}{pointer}     :: CS\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}}
\DoxyCodeLine{650   \textcolor{keywordtype}{logical},             \textcolor{keywordtype}{optional},            \textcolor{keywordtype}{intent(in)}  :: full\_halos\textcolor{comment}{ !< If true, do the calculation}}
\DoxyCodeLine{651 \textcolor{comment}{                                                                      !! over the entire computational domain.}}
\DoxyCodeLine{652   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: better\_speed\_est\textcolor{comment}{ !< If true, use a more robust estimate of the first}}
\DoxyCodeLine{653 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}}
\DoxyCodeLine{654 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: min\_speed\textcolor{comment}{ !< If present, set a floor in the first mode speed}}
\DoxyCodeLine{655 \textcolor{comment}{                                     !! below which 0 is returned [L T-\/1 ~> m s-\/1].}}
\DoxyCodeLine{656 \textcolor{keywordtype}{  real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}}
\DoxyCodeLine{657 \textcolor{comment}{                                     !! wave speeds [nondim]}}
\DoxyCodeLine{658 }
\DoxyCodeLine{659   \textcolor{comment}{! Local variables}}
\DoxyCodeLine{660 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: \&}
\DoxyCodeLine{661     dRho\_dT, \&    \textcolor{comment}{! Partial derivative of density with temperature [R degC-\/1 ~> kg m-\/3 degC-\/1]}}
\DoxyCodeLine{662     dRho\_dS, \&    \textcolor{comment}{! Partial derivative of density with salinity [R ppt-\/1 ~> kg m-\/3 ppt-\/1]}}
\DoxyCodeLine{663     pres, \&       \textcolor{comment}{! Interface pressure [R L2 T-\/2 ~> Pa]}}
\DoxyCodeLine{664     T\_int, \&      \textcolor{comment}{! Temperature interpolated to interfaces [degC]}}
\DoxyCodeLine{665     S\_int, \&      \textcolor{comment}{! Salinity interpolated to interfaces [ppt]}}
\DoxyCodeLine{666     H\_top, \&      \textcolor{comment}{! The distance of each filtered interface from the ocean surface [Z ~> m]}}
\DoxyCodeLine{667     H\_bot, \&      \textcolor{comment}{! The distance of each filtered interface from the bottom [Z ~> m]}}
\DoxyCodeLine{668     gprime        \textcolor{comment}{! The reduced gravity across each interface [L2 Z-\/1 T-\/2 ~> m s-\/2].}}
\DoxyCodeLine{669 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G),SZI\_(G))} :: \&}
\DoxyCodeLine{670     Hf, \&         \textcolor{comment}{! Layer thicknesses after very thin layers are combined [Z ~> m]}}
\DoxyCodeLine{671     Tf, \&         \textcolor{comment}{! Layer temperatures after very thin layers are combined [degC]}}
\DoxyCodeLine{672     Sf, \&         \textcolor{comment}{! Layer salinities after very thin layers are combined [ppt]}}
\DoxyCodeLine{673     Rf            \textcolor{comment}{! Layer densities after very thin layers are combined [R ~> kg m-\/3]}}
\DoxyCodeLine{674 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: \&}
\DoxyCodeLine{675     Igl, Igu, \&   \textcolor{comment}{! The inverse of the reduced gravity across an interface times}}
\DoxyCodeLine{676                   \textcolor{comment}{! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-\/2 ~> s2 m-\/2].}}
\DoxyCodeLine{677     hc, \&         \textcolor{comment}{! A column of layer thicknesses after convective istabilities are removed [Z ~> m]}}
\DoxyCodeLine{678     tc, \&         \textcolor{comment}{! A column of layer temperatures after convective istabilities are removed [degC]}}
\DoxyCodeLine{679     sc, \&         \textcolor{comment}{! A column of layer salinites after convective istabilities are removed [ppt]}}
\DoxyCodeLine{680     rc            \textcolor{comment}{! A column of layer densities after convective istabilities are removed [R ~> kg m-\/3]}}
\DoxyCodeLine{681 \textcolor{keywordtype}{  real} :: I\_Htot  \textcolor{comment}{! The inverse of the total filtered thicknesses [Z ~> m]}}
\DoxyCodeLine{682 \textcolor{keywordtype}{  real} :: c1\_thresh  \textcolor{comment}{! if c1 is below this value, don't bother calculating}}
\DoxyCodeLine{683                      \textcolor{comment}{! cn values for higher modes [L T-\/1 ~> m s-\/1]}}
\DoxyCodeLine{684 \textcolor{keywordtype}{  real} :: c2\_scale \textcolor{comment}{! A scaling factor for wave speeds to help control the growth of the determinant}}
\DoxyCodeLine{685                    \textcolor{comment}{! and its derivative with lam between rows of the Thomas algorithm solver.  The}}
\DoxyCodeLine{686                    \textcolor{comment}{! exact value should not matter for the final result if it is an even power of 2.}}
\DoxyCodeLine{687 \textcolor{keywordtype}{  real} :: det, ddet       \textcolor{comment}{! determinant \& its derivative of eigen system}}
\DoxyCodeLine{688 \textcolor{keywordtype}{  real} :: lam\_1           \textcolor{comment}{! approximate mode-\/1 eigenvalue [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{689 \textcolor{keywordtype}{  real} :: lam\_n           \textcolor{comment}{! approximate mode-\/n eigenvalue [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{690 \textcolor{keywordtype}{  real} :: dlam            \textcolor{comment}{! The change in estimates of the eigenvalue [T2 L-\/2 ~> s m-\/1]}}
\DoxyCodeLine{691 \textcolor{keywordtype}{  real} :: lamMin          \textcolor{comment}{! minimum lam value for root searching range [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{692 \textcolor{keywordtype}{  real} :: lamMax          \textcolor{comment}{! maximum lam value for root searching range [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{693 \textcolor{keywordtype}{  real} :: lamInc          \textcolor{comment}{! width of moving window for root searching [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{694 \textcolor{keywordtype}{  real} :: det\_l,det\_r     \textcolor{comment}{! determinant value at left and right of window}}
\DoxyCodeLine{695 \textcolor{keywordtype}{  real} :: ddet\_l,ddet\_r   \textcolor{comment}{! derivative of determinant at left and right of window}}
\DoxyCodeLine{696 \textcolor{keywordtype}{  real} :: det\_sub,ddet\_sub\textcolor{comment}{! derivative of determinant at subinterval endpoint}}
\DoxyCodeLine{697 \textcolor{keywordtype}{  real} :: xl,xr           \textcolor{comment}{! lam guesses at left and right of window [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{698 \textcolor{keywordtype}{  real} :: xl\_sub          \textcolor{comment}{! lam guess at left of subinterval window [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{699 \textcolor{keywordtype}{  real},\textcolor{keywordtype}{dimension(nmodes)} :: \&}
\DoxyCodeLine{700           xbl,xbr         \textcolor{comment}{! lam guesses bracketing a zero-\/crossing (root) [T2 L-\/2 ~> s2 m-\/2]}}
\DoxyCodeLine{701   \textcolor{keywordtype}{integer} :: numint       \textcolor{comment}{! number of widows (intervals) in root searching range}}
\DoxyCodeLine{702   \textcolor{keywordtype}{integer} :: nrootsfound  \textcolor{comment}{! number of extra roots found (not including 1st root)}}
\DoxyCodeLine{703 \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{704 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: \&}
\DoxyCodeLine{705     htot, hmin, \&  \textcolor{comment}{! Thicknesses [Z ~> m]}}
\DoxyCodeLine{706     H\_here, \&      \textcolor{comment}{! A thickness [Z ~> m]}}
\DoxyCodeLine{707     HxT\_here, \&    \textcolor{comment}{! A layer integrated temperature [degC Z ~> degC m]}}
\DoxyCodeLine{708     HxS\_here, \&    \textcolor{comment}{! A layer integrated salinity [ppt Z ~> ppt m]}}
\DoxyCodeLine{709     HxR\_here       \textcolor{comment}{! A layer integrated density [R Z ~> kg m-\/2]}}
\DoxyCodeLine{710 \textcolor{keywordtype}{  real} :: speed2\_tot \textcolor{comment}{! overestimate of the mode-\/1 speed squared [L2 T-\/2 ~> m2 s-\/2]}}
\DoxyCodeLine{711 \textcolor{keywordtype}{  real} :: speed2\_min \textcolor{comment}{! minimum mode speed (squared) to consider in root searching [L2 T-\/2 ~> m2 s-\/2]}}
\DoxyCodeLine{712 \textcolor{keywordtype}{  real} :: cg1\_min2 \textcolor{comment}{! A floor in the squared first mode speed below which 0 is returned [L2 T-\/2 ~> m2 s-\/2]}}
\DoxyCodeLine{713 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{parameter} :: reduct\_factor = 0.5}
\DoxyCodeLine{714                      \textcolor{comment}{! A factor used in setting speed2\_min [nondim]}}
\DoxyCodeLine{715 \textcolor{keywordtype}{  real} :: I\_Hnew   \textcolor{comment}{! The inverse of a new layer thickness [Z-\/1 ~> m-\/1]}}
\DoxyCodeLine{716 \textcolor{keywordtype}{  real} :: drxh\_sum \textcolor{comment}{! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-\/2]}}
\DoxyCodeLine{717 \textcolor{keywordtype}{  real}, \textcolor{keywordtype}{pointer}, \textcolor{keywordtype}{dimension(:,:,:)} :: T => null(), s => null()}
\DoxyCodeLine{718 \textcolor{keywordtype}{  real} :: g\_Rho0   \textcolor{comment}{! G\_Earth/Rho0 [L2 T-\/2 Z-\/1 R-\/1 ~> m4 s-\/2 kg-\/1].}}
\DoxyCodeLine{719 \textcolor{keywordtype}{  real} :: tol\_Hfrac  \textcolor{comment}{! Layers that together are smaller than this fraction of}}
\DoxyCodeLine{720                      \textcolor{comment}{! the total water column can be merged for efficiency.}}
\DoxyCodeLine{721 \textcolor{keywordtype}{  real} :: min\_h\_frac \textcolor{comment}{! tol\_Hfrac divided by the total number of layers [nondim].}}
\DoxyCodeLine{722 \textcolor{keywordtype}{  real} :: tol\_solve  \textcolor{comment}{! The fractional tolerance with which to solve for the wave speeds [nondim].}}
\DoxyCodeLine{723 \textcolor{keywordtype}{  real} :: tol\_merge  \textcolor{comment}{! The fractional change in estimated wave speed that is allowed}}
\DoxyCodeLine{724                      \textcolor{comment}{! when deciding to merge layers in the calculation [nondim]}}
\DoxyCodeLine{725   \textcolor{keywordtype}{integer} :: kf(SZI\_(G)) \textcolor{comment}{! The number of active layers after filtering.}}
\DoxyCodeLine{726   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: max\_itt = 10}
\DoxyCodeLine{727   \textcolor{keywordtype}{logical} :: use\_EOS    \textcolor{comment}{! If true, density is calculated from T \& S using the equation of state.}}
\DoxyCodeLine{728   \textcolor{keywordtype}{logical} :: better\_est \textcolor{comment}{! If true, use an improved estimate of the first mode internal wave speed.}}
\DoxyCodeLine{729   \textcolor{keywordtype}{logical} :: merge      \textcolor{comment}{! If true, merge the current layer with the one above.}}
\DoxyCodeLine{730   \textcolor{keywordtype}{integer} :: nsub       \textcolor{comment}{! number of subintervals used for root finding}}
\DoxyCodeLine{731   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: sub\_it\_max = 4}
\DoxyCodeLine{732                         \textcolor{comment}{! maximum number of times to subdivide interval}}
\DoxyCodeLine{733                         \textcolor{comment}{! for root finding (\# intervals = 2**sub\_it\_max)}}
\DoxyCodeLine{734   \textcolor{keywordtype}{logical} :: sub\_rootfound \textcolor{comment}{! if true, subdivision has located root}}
\DoxyCodeLine{735   \textcolor{keywordtype}{integer} :: kc         \textcolor{comment}{! The number of layers in the column after merging}}
\DoxyCodeLine{736   \textcolor{keywordtype}{integer} :: sub, sub\_it}
\DoxyCodeLine{737   \textcolor{keywordtype}{integer} :: i, j, k, k2, itt, is, ie, js, je, nz, row, iint, m, ig, jg}
\DoxyCodeLine{738 }
\DoxyCodeLine{739   is = g\%isc ; ie = g\%iec ; js = g\%jsc ; je = g\%jec ; nz = g\%ke}
\DoxyCodeLine{740 }
\DoxyCodeLine{741   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{742     \textcolor{keywordflow}{if} (.not. \textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_wave\_speed: "}// \&}
\DoxyCodeLine{743            \textcolor{stringliteral}{"Module must be initialized before it is used."})}
\DoxyCodeLine{744 \textcolor{keywordflow}{  endif}}
\DoxyCodeLine{745 }
\DoxyCodeLine{746   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(full\_halos)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (full\_halos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{747     is = g\%isd ; ie = g\%ied ; js = g\%jsd ; je = g\%jed}
\DoxyCodeLine{748 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{749 }
\DoxyCodeLine{750   s => tv\%S ; t => tv\%T}
\DoxyCodeLine{751   g\_rho0 = gv\%g\_Earth / gv\%Rho0}
\DoxyCodeLine{752   \textcolor{comment}{! Simplifying the following could change answers at roundoff.}}
\DoxyCodeLine{753   z\_to\_pres = gv\%Z\_to\_H * (gv\%H\_to\_RZ * gv\%g\_Earth)}
\DoxyCodeLine{754   use\_eos = \textcolor{keyword}{associated}(tv\%eqn\_of\_state)}
\DoxyCodeLine{755   c1\_thresh = 0.01*us\%m\_s\_to\_L\_T}
\DoxyCodeLine{756   c2\_scale = us\%m\_s\_to\_L\_T**2 / 4096.0**2 \textcolor{comment}{! Other powers of 2 give identical results.}}
\DoxyCodeLine{757 }
\DoxyCodeLine{758   better\_est = .false. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) better\_est = cs\%better\_cg1\_est}
\DoxyCodeLine{759   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(better\_speed\_est)) better\_est = better\_speed\_est}
\DoxyCodeLine{760   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{761     tol\_solve = 0.001 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) tol\_solve = cs\%wave\_speed\_tol}
\DoxyCodeLine{762     \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(wave\_speed\_tol)) tol\_solve = wave\_speed\_tol}
\DoxyCodeLine{763     tol\_hfrac  = 0.1*tol\_solve ; tol\_merge = tol\_solve / real(nz)}
\DoxyCodeLine{764   \textcolor{keywordflow}{else}}
\DoxyCodeLine{765     tol\_solve = 0.001 ; tol\_hfrac  = 0.0001 ; tol\_merge = 0.001}
\DoxyCodeLine{766 \textcolor{keywordflow}{  endif}}
\DoxyCodeLine{767   cg1\_min2 = 0.0 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) cg1\_min2 = cs\%min\_speed2}
\DoxyCodeLine{768   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(min\_speed)) cg1\_min2 = min\_speed**2}
\DoxyCodeLine{769 }
\DoxyCodeLine{770   \textcolor{comment}{! Zero out all wave speeds.  Values over land or for columns that are too weakly stratified}}
\DoxyCodeLine{771   \textcolor{comment}{! are not changed from this zero value.}}
\DoxyCodeLine{772   cn(:,:,:) = 0.0}
\DoxyCodeLine{773 }
\DoxyCodeLine{774   min\_h\_frac = tol\_hfrac / real(nz)}
\DoxyCodeLine{775   \textcolor{comment}{!\$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,min\_h\_frac,use\_EOS,T,S, \&}}
\DoxyCodeLine{776   \textcolor{comment}{!\$OMP                                     Z\_to\_pres,tv,cn,g\_Rho0,nmodes,cg1\_min2,better\_est, \&}}
\DoxyCodeLine{777   \textcolor{comment}{!\$OMP                                     c1\_thresh,tol\_solve,tol\_merge)}}
\DoxyCodeLine{778   \textcolor{keywordflow}{do} j=js,je}
\DoxyCodeLine{779     \textcolor{comment}{!   First merge very thin layers with the one above (or below if they are}}
\DoxyCodeLine{780     \textcolor{comment}{! at the top).  This also transposes the row order so that columns can}}
\DoxyCodeLine{781     \textcolor{comment}{! be worked upon one at a time.}}
\DoxyCodeLine{782     \textcolor{keywordflow}{do} i=is,ie ; htot(i) = 0.0 ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{783     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv\%H\_to\_Z ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{784 }
\DoxyCodeLine{785     \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{786       hmin(i) = htot(i)*min\_h\_frac ; kf(i) = 1 ; h\_here(i) = 0.0}
\DoxyCodeLine{787       hxt\_here(i) = 0.0 ; hxs\_here(i) = 0.0 ; hxr\_here(i) = 0.0}
\DoxyCodeLine{788 \textcolor{keywordflow}{    enddo}}
\DoxyCodeLine{789     \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{790       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{791         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv\%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}}
\DoxyCodeLine{792           hf(kf(i),i) = h\_here(i)}
\DoxyCodeLine{793           tf(kf(i),i) = hxt\_here(i) / h\_here(i)}
\DoxyCodeLine{794           sf(kf(i),i) = hxs\_here(i) / h\_here(i)}
\DoxyCodeLine{795           kf(i) = kf(i) + 1}
\DoxyCodeLine{796 }
\DoxyCodeLine{797           \textcolor{comment}{! Start a new layer}}
\DoxyCodeLine{798           h\_here(i) = h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{799           hxt\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*t(i,j,k)}
\DoxyCodeLine{800           hxs\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*s(i,j,k)}
\DoxyCodeLine{801         \textcolor{keywordflow}{else}}
\DoxyCodeLine{802           h\_here(i) = h\_here(i) + h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{803           hxt\_here(i) = hxt\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*t(i,j,k)}
\DoxyCodeLine{804           hxs\_here(i) = hxs\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*s(i,j,k)}
\DoxyCodeLine{805 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{806 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{807       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{808         hf(kf(i),i) = h\_here(i)}
\DoxyCodeLine{809         tf(kf(i),i) = hxt\_here(i) / h\_here(i)}
\DoxyCodeLine{810         sf(kf(i),i) = hxs\_here(i) / h\_here(i)}
\DoxyCodeLine{811 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{812     \textcolor{keywordflow}{else}}
\DoxyCodeLine{813       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{814         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv\%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}}
\DoxyCodeLine{815           hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)}
\DoxyCodeLine{816           kf(i) = kf(i) + 1}
\DoxyCodeLine{817 }
\DoxyCodeLine{818           \textcolor{comment}{! Start a new layer}}
\DoxyCodeLine{819           h\_here(i) = h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{820           hxr\_here(i) = (h(i,j,k)*gv\%H\_to\_Z)*gv\%Rlay(k)}
\DoxyCodeLine{821         \textcolor{keywordflow}{else}}
\DoxyCodeLine{822           h\_here(i) = h\_here(i) + h(i,j,k)*gv\%H\_to\_Z}
\DoxyCodeLine{823           hxr\_here(i) = hxr\_here(i) + (h(i,j,k)*gv\%H\_to\_Z)*gv\%Rlay(k)}
\DoxyCodeLine{824 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{825 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{826       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{827         hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)}
\DoxyCodeLine{828 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{829 \textcolor{keywordflow}{    endif}}
\DoxyCodeLine{830 }
\DoxyCodeLine{831     \textcolor{comment}{! From this point, we can work on individual columns without causing memory to have page faults.}}
\DoxyCodeLine{832     \textcolor{keywordflow}{do} i=is,ie}
\DoxyCodeLine{833       \textcolor{keywordflow}{if} (g\%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then}}
\DoxyCodeLine{834         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{835           pres(1) = 0.0 ; h\_top(1) = 0.0}
\DoxyCodeLine{836           \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{837             pres(k) = pres(k-\/1) + z\_to\_pres*hf(k-\/1,i)}
\DoxyCodeLine{838             t\_int(k) = 0.5*(tf(k,i)+tf(k-\/1,i))}
\DoxyCodeLine{839             s\_int(k) = 0.5*(sf(k,i)+sf(k-\/1,i))}
\DoxyCodeLine{840             h\_top(k) = h\_top(k-\/1) + hf(k-\/1,i)}
\DoxyCodeLine{841 \textcolor{keywordflow}{          enddo}}
\DoxyCodeLine{842           \textcolor{keyword}{call }calculate\_density\_derivs(t\_int, s\_int, pres, drho\_dt, drho\_ds, \&}
\DoxyCodeLine{843                                         tv\%eqn\_of\_state, (/2,kf(i)/) )}
\DoxyCodeLine{844 }
\DoxyCodeLine{845           \textcolor{comment}{! Sum the reduced gravities to find out how small a density difference is negligibly small.}}
\DoxyCodeLine{846           drxh\_sum = 0.0}
\DoxyCodeLine{847           \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{848             \textcolor{comment}{! This is an estimate that is correct for the non-\/EBT mode for 2 or 3 layers, or for}}
\DoxyCodeLine{849             \textcolor{comment}{! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.}}
\DoxyCodeLine{850             \textcolor{comment}{! For a uniform stratification and a huge number of layers uniformly distributed in}}
\DoxyCodeLine{851             \textcolor{comment}{! density, this estimate is too large (as is desired) by a factor of pi\string^2/6 ~= 1.64.}}
\DoxyCodeLine{852             \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{853               i\_htot = 1.0 / (h\_top(kf(i)) + hf(kf(i),i))  \textcolor{comment}{! = 1.0 / (H\_top(K) + H\_bot(K)) for all K.}}
\DoxyCodeLine{854               h\_bot(kf(i)+1) = 0.0}
\DoxyCodeLine{855               \textcolor{keywordflow}{do} k=kf(i),2,-\/1}
\DoxyCodeLine{856                 h\_bot(k) = h\_bot(k+1) + hf(k,i)}
\DoxyCodeLine{857                 drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * \&}
\DoxyCodeLine{858                     max(0.0, drho\_dt(k)*(tf(k,i)-\/tf(k-\/1,i)) + drho\_ds(k)*(sf(k,i)-\/sf(k-\/1,i)))}
\DoxyCodeLine{859 \textcolor{keywordflow}{              enddo}}
\DoxyCodeLine{860 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{861           \textcolor{keywordflow}{else}}
\DoxyCodeLine{862             \textcolor{comment}{! This estimate is problematic in that it goes like 1/nz for a large number of layers,}}
\DoxyCodeLine{863             \textcolor{comment}{! but it is an overestimate (as desired) for a small number of layers, by at a factor}}
\DoxyCodeLine{864             \textcolor{comment}{! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.}}
\DoxyCodeLine{865             \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{866               drxh\_sum = drxh\_sum + 0.5*(hf(k-\/1,i)+hf(k,i)) * \&}
\DoxyCodeLine{867                   max(0.0, drho\_dt(k)*(tf(k,i)-\/tf(k-\/1,i)) + drho\_ds(k)*(sf(k,i)-\/sf(k-\/1,i)))}
\DoxyCodeLine{868 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{869 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{870         \textcolor{keywordflow}{else}}
\DoxyCodeLine{871           drxh\_sum = 0.0}
\DoxyCodeLine{872           \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{873             h\_top(1) = 0.0}
\DoxyCodeLine{874             \textcolor{keywordflow}{do} k=2,kf(i) ; h\_top(k) = h\_top(k-\/1) + hf(k-\/1,i) ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{875               \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}}
\DoxyCodeLine{876               i\_htot = 1.0 / (h\_top(kf(i)) + hf(kf(i),i))  \textcolor{comment}{! = 1.0 / (H\_top(K) + H\_bot(K)) for all K.}}
\DoxyCodeLine{877               h\_bot(kf(i)+1) = 0.0}
\DoxyCodeLine{878               \textcolor{keywordflow}{do} k=kf(i),2,-\/1}
\DoxyCodeLine{879                 h\_bot(k) = h\_bot(k+1) + hf(k,i)}
\DoxyCodeLine{880                 drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * max(0.0,rf(k,i)-\/rf(k-\/1,i))}
\DoxyCodeLine{881 \textcolor{keywordflow}{              enddo}}
\DoxyCodeLine{882 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{883           \textcolor{keywordflow}{else}}
\DoxyCodeLine{884             \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{885               drxh\_sum = drxh\_sum + 0.5*(hf(k-\/1,i)+hf(k,i)) * max(0.0,rf(k,i)-\/rf(k-\/1,i))}
\DoxyCodeLine{886 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{887 \textcolor{keywordflow}{          endif}}
\DoxyCodeLine{888 \textcolor{keywordflow}{        endif}}
\DoxyCodeLine{889 }
\DoxyCodeLine{890         \textcolor{comment}{!   Find gprime across each internal interface, taking care of convective}}
\DoxyCodeLine{891         \textcolor{comment}{! instabilities by merging layers.}}
\DoxyCodeLine{892         \textcolor{keywordflow}{if} (g\_rho0 * drxh\_sum > cg1\_min2) \textcolor{keywordflow}{then}}
\DoxyCodeLine{893           \textcolor{comment}{! Merge layers to eliminate convective instabilities or exceedingly}}
\DoxyCodeLine{894           \textcolor{comment}{! small reduced gravities.  Merging layers reduces the estimated wave speed by}}
\DoxyCodeLine{895           \textcolor{comment}{! (rho(2)-\/rho(1))*h(1)*h(2) / H\_tot.}}
\DoxyCodeLine{896           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}}
\DoxyCodeLine{897             kc = 1}
\DoxyCodeLine{898             hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)}
\DoxyCodeLine{899             \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{900               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{901                 merge = ((drho\_dt(k)*(tf(k,i)-\/tc(kc)) + drho\_ds(k)*(sf(k,i)-\/sc(kc))) * \&}
\DoxyCodeLine{902                          ((hc(kc) * hf(k,i))*i\_htot) < 2.0 * tol\_merge*drxh\_sum)}
\DoxyCodeLine{903               \textcolor{keywordflow}{else}}
\DoxyCodeLine{904                 merge = ((drho\_dt(k)*(tf(k,i)-\/tc(kc)) + drho\_ds(k)*(sf(k,i)-\/sc(kc))) * \&}
\DoxyCodeLine{905                          (hc(kc) + hf(k,i)) < 2.0 * tol\_merge*drxh\_sum)}
\DoxyCodeLine{906 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{907               \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{908                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}}
\DoxyCodeLine{909                 i\_hnew = 1.0 / (hc(kc) + hf(k,i))}
\DoxyCodeLine{910                 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i\_hnew}
\DoxyCodeLine{911                 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i\_hnew}
\DoxyCodeLine{912                 hc(kc) = (hc(kc) + hf(k,i))}
\DoxyCodeLine{913                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}}
\DoxyCodeLine{914                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}}
\DoxyCodeLine{915                 \textcolor{comment}{! far back we go.}}
\DoxyCodeLine{916                 \textcolor{keywordflow}{do} k2=kc,2,-\/1}
\DoxyCodeLine{917                   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{918                     merge = ((drho\_dt(k2)*(tc(k2)-\/tc(k2-\/1)) + drho\_ds(k2)*(sc(k2)-\/sc(k2-\/1))) * \&}
\DoxyCodeLine{919                              ((hc(k2) * hc(k2-\/1))*i\_htot) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{920                   \textcolor{keywordflow}{else}}
\DoxyCodeLine{921                     merge = ((drho\_dt(k2)*(tc(k2)-\/tc(k2-\/1)) + drho\_ds(k2)*(sc(k2)-\/sc(k2-\/1))) * \&}
\DoxyCodeLine{922                              (hc(k2) + hc(k2-\/1)) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{923 \textcolor{keywordflow}{                  endif}}
\DoxyCodeLine{924                   \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{925                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}}
\DoxyCodeLine{926                     i\_hnew = 1.0 / (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{927                     tc(kc-\/1) = (hc(kc)*tc(kc) + hc(kc-\/1)*tc(kc-\/1)) * i\_hnew}
\DoxyCodeLine{928                     sc(kc-\/1) = (hc(kc)*sc(kc) + hc(kc-\/1)*sc(kc-\/1)) * i\_hnew}
\DoxyCodeLine{929                     hc(kc-\/1) = (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{930                     kc = kc -\/ 1}
\DoxyCodeLine{931                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{932 \textcolor{keywordflow}{                enddo}}
\DoxyCodeLine{933               \textcolor{keywordflow}{else}}
\DoxyCodeLine{934                 \textcolor{comment}{! Add a new layer to the column.}}
\DoxyCodeLine{935                 kc = kc + 1}
\DoxyCodeLine{936                 drho\_ds(kc) = drho\_ds(k) ; drho\_dt(kc) = drho\_dt(k)}
\DoxyCodeLine{937                 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)}
\DoxyCodeLine{938 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{939 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{940             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}}
\DoxyCodeLine{941             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-\/Boussinesq.}}
\DoxyCodeLine{942               gprime(k) = g\_rho0 * (drho\_dt(k)*(tc(k)-\/tc(k-\/1)) + drho\_ds(k)*(sc(k)-\/sc(k-\/1)))}
\DoxyCodeLine{943 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{944           \textcolor{keywordflow}{else}  \textcolor{comment}{! .not.use\_EOS}}
\DoxyCodeLine{945             \textcolor{comment}{! Do the same with density directly...}}
\DoxyCodeLine{946             kc = 1}
\DoxyCodeLine{947             hc(1) = hf(1,i) ; rc(1) = rf(1,i)}
\DoxyCodeLine{948             \textcolor{keywordflow}{do} k=2,kf(i)}
\DoxyCodeLine{949               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{950                 merge = ((rf(k,i) -\/ rc(kc)) * ((hc(kc) * hf(k,i))*i\_htot) < 2.0 * tol\_merge*drxh\_sum)}
\DoxyCodeLine{951               \textcolor{keywordflow}{else}}
\DoxyCodeLine{952                 merge = ((rf(k,i) -\/ rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol\_merge*drxh\_sum)}
\DoxyCodeLine{953 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{954               \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{955                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}}
\DoxyCodeLine{956                 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))}
\DoxyCodeLine{957                 hc(kc) = (hc(kc) + hf(k,i))}
\DoxyCodeLine{958                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}}
\DoxyCodeLine{959                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}}
\DoxyCodeLine{960                 \textcolor{comment}{! far back we go.}}
\DoxyCodeLine{961                 \textcolor{keywordflow}{do} k2=kc,2,-\/1}
\DoxyCodeLine{962                   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{963                     merge = ((rc(k2)-\/rc(k2-\/1)) * ((hc(k2) * hc(k2-\/1))*i\_htot) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{964                   \textcolor{keywordflow}{else}}
\DoxyCodeLine{965                     merge = ((rc(k2)-\/rc(k2-\/1)) * (hc(k2)+hc(k2-\/1)) < tol\_merge*drxh\_sum)}
\DoxyCodeLine{966 \textcolor{keywordflow}{                  endif}}
\DoxyCodeLine{967                   \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}}
\DoxyCodeLine{968                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}}
\DoxyCodeLine{969                     rc(kc-\/1) = (hc(kc)*rc(kc) + hc(kc-\/1)*rc(kc-\/1)) / (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{970                     hc(kc-\/1) = (hc(kc) + hc(kc-\/1))}
\DoxyCodeLine{971                     kc = kc -\/ 1}
\DoxyCodeLine{972                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}}
\DoxyCodeLine{973 \textcolor{keywordflow}{                enddo}}
\DoxyCodeLine{974               \textcolor{keywordflow}{else}}
\DoxyCodeLine{975                 \textcolor{comment}{! Add a new layer to the column.}}
\DoxyCodeLine{976                 kc = kc + 1}
\DoxyCodeLine{977                 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)}
\DoxyCodeLine{978 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{979 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{980             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}}
\DoxyCodeLine{981             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-\/Boussinesq.}}
\DoxyCodeLine{982               gprime(k) = g\_rho0 * (rc(k)-\/rc(k-\/1))}
\DoxyCodeLine{983 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{984 \textcolor{keywordflow}{          endif}  \textcolor{comment}{! use\_EOS}}
\DoxyCodeLine{985 }
\DoxyCodeLine{986           \textcolor{comment}{!-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/NOW FIND WAVE SPEEDS-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/}}
\DoxyCodeLine{987           \textcolor{comment}{! ig = i + G\%idg\_offset ; jg = j + G\%jdg\_offset}}
\DoxyCodeLine{988           \textcolor{comment}{!   Sum the contributions from all of the interfaces to give an over-\/estimate}}
\DoxyCodeLine{989           \textcolor{comment}{! of the first-\/mode wave speed.  Also populate Igl and Igu which are the}}
\DoxyCodeLine{990           \textcolor{comment}{! non-\/leading diagonals of the tridiagonal matrix.}}
\DoxyCodeLine{991           \textcolor{keywordflow}{if} (kc >= 2) \textcolor{keywordflow}{then}}
\DoxyCodeLine{992             \textcolor{comment}{! initialize speed2\_tot}}
\DoxyCodeLine{993             speed2\_tot = 0.0}
\DoxyCodeLine{994             \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{995               h\_top(1) = 0.0 ; h\_bot(kc+1) = 0.0}
\DoxyCodeLine{996               \textcolor{keywordflow}{do} k=2,kc+1 ; h\_top(k) = h\_top(k-\/1) + hc(k-\/1) ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{997               \textcolor{keywordflow}{do} k=kc,2,-\/1 ; h\_bot(k) = h\_bot(k+1) + hc(k) ;\textcolor{keywordflow}{ enddo}}
\DoxyCodeLine{998               i\_htot = 0.0 ; \textcolor{keywordflow}{if} (h\_top(kc+1) > 0.0) i\_htot = 1.0 / h\_top(kc+1)}
\DoxyCodeLine{999 \textcolor{keywordflow}{            endif}}
\DoxyCodeLine{1000 }
\DoxyCodeLine{1001             \textcolor{comment}{! Calculate Igu, Igl, depth, and N2 at each interior interface}}
\DoxyCodeLine{1002             \textcolor{comment}{! [excludes surface (K=1) and bottom (K=kc+1)]}}
\DoxyCodeLine{1003             \textcolor{keywordflow}{do} k=2,kc}
\DoxyCodeLine{1004               igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-\/1))}
\DoxyCodeLine{1005               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1006                 speed2\_tot = speed2\_tot + gprime(k)*((h\_top(k) * h\_bot(k)) * i\_htot)}
\DoxyCodeLine{1007               \textcolor{keywordflow}{else}}
\DoxyCodeLine{1008                 speed2\_tot = speed2\_tot + gprime(k)*(hc(k-\/1)+hc(k))}
\DoxyCodeLine{1009 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{1010 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{1011 }
\DoxyCodeLine{1012             \textcolor{comment}{! Under estimate the first eigenvalue (overestimate the speed) to start with.}}
\DoxyCodeLine{1013             lam\_1 = 1.0 / speed2\_tot}
\DoxyCodeLine{1014 }
\DoxyCodeLine{1015             \textcolor{comment}{! Find the first eigen value}}
\DoxyCodeLine{1016             \textcolor{keywordflow}{do} itt=1,max\_itt}
\DoxyCodeLine{1017               \textcolor{comment}{! calculate the determinant of (A-\/lam\_1*I)}}
\DoxyCodeLine{1018               \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lam\_1, det, ddet, row\_scale=c2\_scale)}
\DoxyCodeLine{1019 }
\DoxyCodeLine{1020               \textcolor{comment}{! If possible, use Newton's method iteration to find a new estimate of lam\_1}}
\DoxyCodeLine{1021               \textcolor{comment}{!det = det\_it(itt) ; ddet = ddet\_it(itt)}}
\DoxyCodeLine{1022               \textcolor{keywordflow}{if} ((ddet >= 0.0) .or. (-\/det > -\/0.5*lam\_1*ddet)) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1023                 \textcolor{comment}{! lam\_1 was not an under-\/estimate, as intended, so Newton's method}}
\DoxyCodeLine{1024                 \textcolor{comment}{! may not be reliable; lam\_1 must be reduced, but not by more than half.}}
\DoxyCodeLine{1025                 lam\_1 = 0.5 * lam\_1}
\DoxyCodeLine{1026                 dlam = -\/lam\_1}
\DoxyCodeLine{1027               \textcolor{keywordflow}{else}  \textcolor{comment}{! Newton's method is OK.}}
\DoxyCodeLine{1028                 dlam = -\/ det / ddet}
\DoxyCodeLine{1029                 lam\_1 = lam\_1 + dlam}
\DoxyCodeLine{1030 \textcolor{keywordflow}{              endif}}
\DoxyCodeLine{1031 }
\DoxyCodeLine{1032               \textcolor{keywordflow}{if} (abs(dlam) < tol\_solve*lam\_1) \textcolor{keywordflow}{exit}}
\DoxyCodeLine{1033 \textcolor{keywordflow}{            enddo}}
\DoxyCodeLine{1034 }
\DoxyCodeLine{1035             \textcolor{keywordflow}{if} (lam\_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam\_1)}
\DoxyCodeLine{1036 }
\DoxyCodeLine{1037             \textcolor{comment}{! Find other eigen values if c1 is of significant magnitude, > cn\_thresh}}
\DoxyCodeLine{1038             nrootsfound = 0    \textcolor{comment}{! number of extra roots found (not including 1st root)}}
\DoxyCodeLine{1039             \textcolor{keywordflow}{if} (nmodes>1 .and. kc>=nmodes+1 .and. cn(i,j,1)>c1\_thresh) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1040               \textcolor{comment}{! Set the the range to look for the other desired eigen values}}
\DoxyCodeLine{1041               \textcolor{comment}{! set min value just greater than the 1st root (found above)}}
\DoxyCodeLine{1042               lammin = lam\_1*(1.0 + tol\_solve)}
\DoxyCodeLine{1043               \textcolor{comment}{! set max value based on a low guess at wavespeed for highest mode}}
\DoxyCodeLine{1044               speed2\_min = (reduct\_factor*cn(i,j,1)/real(nmodes))**2}
\DoxyCodeLine{1045               lammax = 1.0 / speed2\_min}
\DoxyCodeLine{1046               \textcolor{comment}{! set width of interval (not sure about this -\/ BDM)}}
\DoxyCodeLine{1047               laminc = 0.5*lam\_1}
\DoxyCodeLine{1048               \textcolor{comment}{! set number of intervals within search range}}
\DoxyCodeLine{1049               numint = nint((lammax -\/ lammin)/laminc)}
\DoxyCodeLine{1050 }
\DoxyCodeLine{1051               \textcolor{comment}{!   Find intervals containing zero-\/crossings (roots) of the determinant}}
\DoxyCodeLine{1052               \textcolor{comment}{! that are beyond the first root}}
\DoxyCodeLine{1053 }
\DoxyCodeLine{1054               \textcolor{comment}{! find det\_l of first interval (det at left endpoint)}}
\DoxyCodeLine{1055               \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lammin, det\_l, ddet\_l, row\_scale=c2\_scale)}
\DoxyCodeLine{1056               \textcolor{comment}{! move interval window looking for zero-\/crossings************************}}
\DoxyCodeLine{1057               \textcolor{keywordflow}{do} iint=1,numint}
\DoxyCodeLine{1058                 xr = lammin + laminc * iint}
\DoxyCodeLine{1059                 xl = xr -\/ laminc}
\DoxyCodeLine{1060                 \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, xr, det\_r, ddet\_r, row\_scale=c2\_scale)}
\DoxyCodeLine{1061                 \textcolor{keywordflow}{if} (det\_l*det\_r < 0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! if function changes sign}}
\DoxyCodeLine{1062                   \textcolor{keywordflow}{if} (det\_l*ddet\_l < 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! if function at left is headed to zero}}
\DoxyCodeLine{1063                     nrootsfound = nrootsfound + 1}
\DoxyCodeLine{1064                     xbl(nrootsfound) = xl}
\DoxyCodeLine{1065                     xbr(nrootsfound) = xr}
\DoxyCodeLine{1066                   \textcolor{keywordflow}{else}}
\DoxyCodeLine{1067                     \textcolor{comment}{!   function changes sign but has a local max/min in interval,}}
\DoxyCodeLine{1068                     \textcolor{comment}{! try subdividing interval as many times as necessary (or sub\_it\_max).}}
\DoxyCodeLine{1069                     \textcolor{comment}{! loop that increases number of subintervals:}}
\DoxyCodeLine{1070                     \textcolor{comment}{!call MOM\_error(WARNING, "determinant changes sign"// \&}}
\DoxyCodeLine{1071                     \textcolor{comment}{!            "but has a local max/min in interval;"//\&}}
\DoxyCodeLine{1072                     \textcolor{comment}{!            " reduce increment in lam.")}}
\DoxyCodeLine{1073                     \textcolor{comment}{! begin subdivision loop -\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/}}
\DoxyCodeLine{1074                     sub\_rootfound = .false. \textcolor{comment}{! initialize}}
\DoxyCodeLine{1075                     \textcolor{keywordflow}{do} sub\_it=1,sub\_it\_max}
\DoxyCodeLine{1076                       nsub = 2**sub\_it \textcolor{comment}{! number of subintervals; nsub=2,4,8,...}}
\DoxyCodeLine{1077                       \textcolor{comment}{! loop over each subinterval:}}
\DoxyCodeLine{1078                       \textcolor{keywordflow}{do} sub=1,nsub-\/1,2 \textcolor{comment}{! only check odds; sub = 1; 1,3; 1,3,5,7;...}}
\DoxyCodeLine{1079                         xl\_sub = xl + laminc/(nsub)*sub}
\DoxyCodeLine{1080                         \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, xl\_sub, det\_sub, ddet\_sub, \&}
\DoxyCodeLine{1081                                          row\_scale=c2\_scale)}
\DoxyCodeLine{1082                         \textcolor{keywordflow}{if} (det\_sub*det\_r < 0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! if function changes sign}}
\DoxyCodeLine{1083                           \textcolor{keywordflow}{if} (det\_sub*ddet\_sub < 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! if function at left is headed to zero}}
\DoxyCodeLine{1084                             sub\_rootfound = .true.}
\DoxyCodeLine{1085                             nrootsfound = nrootsfound + 1}
\DoxyCodeLine{1086                             xbl(nrootsfound) = xl\_sub}
\DoxyCodeLine{1087                             xbr(nrootsfound) = xr}
\DoxyCodeLine{1088                             \textcolor{keywordflow}{exit} \textcolor{comment}{! exit sub loop}}
\DoxyCodeLine{1089 \textcolor{keywordflow}{                          endif} \textcolor{comment}{! headed toward zero}}
\DoxyCodeLine{1090 \textcolor{keywordflow}{                        endif} \textcolor{comment}{! sign change}}
\DoxyCodeLine{1091 \textcolor{keywordflow}{                      enddo} \textcolor{comment}{! sub-\/loop}}
\DoxyCodeLine{1092                       \textcolor{keywordflow}{if} (sub\_rootfound) \textcolor{keywordflow}{exit} \textcolor{comment}{! root has been found, exit sub\_it loop}}
\DoxyCodeLine{1093                       \textcolor{comment}{!   Otherwise, function changes sign but has a local max/min in one of the}}
\DoxyCodeLine{1094                       \textcolor{comment}{! sub intervals, try subdividing again unless sub\_it\_max has been reached.}}
\DoxyCodeLine{1095                       \textcolor{keywordflow}{if} (sub\_it == sub\_it\_max) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1096                         \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_speed: root not found "}// \&}
\DoxyCodeLine{1097                                        \textcolor{stringliteral}{" after sub\_it\_max subdivisions of original"}// \&}
\DoxyCodeLine{1098                                        \textcolor{stringliteral}{" interval."})}
\DoxyCodeLine{1099 \textcolor{keywordflow}{                      endif} \textcolor{comment}{! sub\_it == sub\_it\_max}}
\DoxyCodeLine{1100 \textcolor{keywordflow}{                    enddo} \textcolor{comment}{! sub\_it-\/loop-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/-\/}}
\DoxyCodeLine{1101 \textcolor{keywordflow}{                  endif} \textcolor{comment}{! det\_l*ddet\_l < 0.0}}
\DoxyCodeLine{1102 \textcolor{keywordflow}{                endif} \textcolor{comment}{! det\_l*det\_r < 0.0}}
\DoxyCodeLine{1103                 \textcolor{comment}{! exit iint-\/loop if all desired roots have been found}}
\DoxyCodeLine{1104                 \textcolor{keywordflow}{if} (nrootsfound >= nmodes-\/1) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1105                   \textcolor{comment}{! exit if all additional roots found}}
\DoxyCodeLine{1106                   \textcolor{keywordflow}{exit}}
\DoxyCodeLine{1107                 \textcolor{keywordflow}{elseif} (iint == numint) \textcolor{keywordflow}{then}}
\DoxyCodeLine{1108                   \textcolor{comment}{! oops, lamMax not large enough -\/ could add code to increase (BDM)}}
\DoxyCodeLine{1109                   \textcolor{comment}{! set unfound modes to zero for now (BDM)}}
\DoxyCodeLine{1110                   \textcolor{comment}{!   cn(i,j,nrootsfound+2:nmodes) = 0.0}}
\DoxyCodeLine{1111                 \textcolor{keywordflow}{else}}
\DoxyCodeLine{1112                   \textcolor{comment}{! else shift interval and keep looking until nmodes or numint is reached}}
\DoxyCodeLine{1113                   det\_l = det\_r}
\DoxyCodeLine{1114                   ddet\_l = ddet\_r}
\DoxyCodeLine{1115 \textcolor{keywordflow}{                endif}}
\DoxyCodeLine{1116 \textcolor{keywordflow}{              enddo} \textcolor{comment}{! iint-\/loop}}
\DoxyCodeLine{1117 }
\DoxyCodeLine{1118               \textcolor{comment}{! Use Newton's method to find the roots within the identified windows}}
\DoxyCodeLine{1119               \textcolor{keywordflow}{do} m=1,nrootsfound \textcolor{comment}{! loop over the root-\/containing widows (excluding 1st mode)}}
\DoxyCodeLine{1120                 lam\_n = xbl(m) \textcolor{comment}{! first guess is left edge of window}}
\DoxyCodeLine{1121                 \textcolor{keywordflow}{do} itt=1,max\_itt}
\DoxyCodeLine{1122                   \textcolor{comment}{! calculate the determinant of (A-\/lam\_n*I)}}
\DoxyCodeLine{1123                   \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lam\_n, det, ddet, row\_scale=c2\_scale)}
\DoxyCodeLine{1124                   \textcolor{comment}{! Use Newton's method to find a new estimate of lam\_n}}
\DoxyCodeLine{1125                   dlam = -\/ det / ddet}
\DoxyCodeLine{1126                   lam\_n = lam\_n + dlam}
\DoxyCodeLine{1127                   \textcolor{keywordflow}{if} (abs(dlam) < tol\_solve*lam\_1)  \textcolor{keywordflow}{exit}}
\DoxyCodeLine{1128 \textcolor{keywordflow}{                enddo} \textcolor{comment}{! itt-\/loop}}
\DoxyCodeLine{1129                 \textcolor{comment}{! calculate nth mode speed}}
\DoxyCodeLine{1130                 \textcolor{keywordflow}{if} (lam\_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam\_n)}
\DoxyCodeLine{1131 \textcolor{keywordflow}{              enddo} \textcolor{comment}{! n-\/loop}}
\DoxyCodeLine{1132 \textcolor{keywordflow}{            endif} \textcolor{comment}{! if nmodes>1 .and. kc>nmodes .and. c1>c1\_thresh}}
\DoxyCodeLine{1133 \textcolor{keywordflow}{          endif} \textcolor{comment}{! if more than 2 layers}}
\DoxyCodeLine{1134 \textcolor{keywordflow}{        endif} \textcolor{comment}{! if drxh\_sum < 0}}
\DoxyCodeLine{1135 \textcolor{keywordflow}{      endif} \textcolor{comment}{! if not land}}
\DoxyCodeLine{1136 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! i-\/loop}}
\DoxyCodeLine{1137 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-\/loop}}
\DoxyCodeLine{1138 }

\end{DoxyCode}
