\hypertarget{namespacemom__wave__speed}{}\section{mom\+\_\+wave\+\_\+speed Module Reference}
\label{namespacemom__wave__speed}\index{mom\+\_\+wave\+\_\+speed@{mom\+\_\+wave\+\_\+speed}}


\subsection{Detailed Description}
Routines for calculating baroclinic wave speeds. \subsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \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}
\subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \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 \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 \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 \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 \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 \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}


\subsection{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}}
\subsubsection{\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{\tt in}  & {\em n} & Number of rows of matrix\\
\hline
\mbox{\tt in}  & {\em a} & Lower diagonal \mbox{[}T2 L-\/2 $\sim$$>$ s2 m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em c} & Upper diagonal \mbox{[}T2 L-\/2 $\sim$$>$ s2 m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em lam} & Scalar subtracted from leading diagonal \mbox{[}T2 L-\/2 $\sim$$>$ s2 m-\/2\mbox{]}\\
\hline
\mbox{\tt 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}
593   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{ !< Number of rows of matrix}
594   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}    :: a\textcolor{comment}{ !< Lower diagonal   [T2 L-2 ~> s2 m-2]}
595   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)}    :: c\textcolor{comment}{ !< Upper diagonal   [T2 L-2 ~> s2 m-2]}
596   \textcolor{keywordtype}{real},               \textcolor{keywordtype}{intent(in)}    :: lam\textcolor{comment}{ !< Scalar subtracted from leading diagonal [T2 L-2 ~> s2 m-2]}
597   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(inout)} :: y\textcolor{comment}{ !< RHS on entry, result on exit}
598 
599   \textcolor{comment}{! Local variables}
600   \textcolor{keywordtype}{real} :: lambda     \textcolor{comment}{! A temporary variable in [T2 L-2 ~> s2 m-2]}
601   \textcolor{keywordtype}{real} :: beta(n)    \textcolor{comment}{! A temporary variable in [T2 L-2 ~> s2 m-2]}
602   \textcolor{keywordtype}{real} :: i\_beta(n)  \textcolor{comment}{! A temporary variable in [L2 T-2 ~> m2 s-2]}
603   \textcolor{keywordtype}{real} :: yy(n)      \textcolor{comment}{! A temporary variable with the same units as y on entry.}
604   \textcolor{keywordtype}{integer} :: k, m
605 
606   lambda = lam
607   beta(1) = (a(1)+c(1)) - lambda
608   \textcolor{keywordflow}{if} (beta(1)==0.) \textcolor{keywordflow}{then} \textcolor{comment}{! lam was chosen too perfectly}
609     \textcolor{comment}{! Change lambda and redo this first row}
610     lambda = (1. + 1.e-5) * lambda
611     beta(1) = (a(1)+c(1)) - lambda
612 \textcolor{keywordflow}{  endif}
613   i\_beta(1) = 1. / beta(1)
614   yy(1) = y(1)
615   \textcolor{keywordflow}{do} k = 2, n
616     beta(k) = ( (a(k)+c(k)) - lambda ) - a(k) * c(k-1) * i\_beta(k-1)
617     \textcolor{comment}{! Perhaps the following 0 needs to become a tolerance to handle underflow?}
618     \textcolor{keywordflow}{if} (beta(k)==0.) \textcolor{keywordflow}{then} \textcolor{comment}{! lam was chosen too perfectly}
619       \textcolor{comment}{! Change lambda and redo everything up to row k}
620       lambda = (1. + 1.e-5) * lambda
621       i\_beta(1) = 1. / ( (a(1)+c(1)) - lambda )
622       \textcolor{keywordflow}{do} m = 2, k
623         i\_beta(m) = 1. / ( ( (a(m)+c(m)) - lambda ) - a(m) * c(m-1) * i\_beta(m-1) )
624         yy(m) = y(m) + a(m) * yy(m-1) * i\_beta(m-1)
625 \textcolor{keywordflow}{      enddo}
626     \textcolor{keywordflow}{else}
627       i\_beta(k) = 1. / beta(k)
628 \textcolor{keywordflow}{    endif}
629     yy(k) = y(k) + a(k) * yy(k-1) * i\_beta(k-1)
630 \textcolor{keywordflow}{  enddo}
631   \textcolor{comment}{! The units of y change by a factor of [L2 T-2] in the following lines.}
632   y(n) = yy(n) * i\_beta(n)
633   \textcolor{keywordflow}{do} k = n-1, 1, -1
634     y(k) = ( yy(k) + c(k) * y(k+1) ) * i\_beta(k)
635 \textcolor{keywordflow}{  enddo}
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}}
\subsubsection{\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{\tt in}  & {\em a} & Lower diagonal of matrix (first entry unused)\\
\hline
\mbox{\tt in}  & {\em c} & Upper diagonal of matrix (last entry unused)\\
\hline
\mbox{\tt in}  & {\em ks} & Starting index to use in determinant\\
\hline
\mbox{\tt in}  & {\em ke} & Ending index to use in determinant\\
\hline
\mbox{\tt in}  & {\em lam} & Value subtracted from b\\
\hline
\mbox{\tt out}  & {\em det} & Determinant\\
\hline
\mbox{\tt out}  & {\em ddet} & Derivative of determinant with lam\\
\hline
\mbox{\tt 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}
1146   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: a\textcolor{comment}{     !< Lower diagonal of matrix (first entry unused)}
1147   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)}, \textcolor{keywordtype}{intent(in)} :: c\textcolor{comment}{     !< Upper diagonal of matrix (last entry unused)}
1148   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)} :: ks\textcolor{comment}{    !< Starting index to use in determinant}
1149   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)} :: ke\textcolor{comment}{    !< Ending index to use in determinant}
1150   \textcolor{keywordtype}{real},               \textcolor{keywordtype}{intent(in)} :: lam\textcolor{comment}{   !< Value subtracted from b}
1151   \textcolor{keywordtype}{real},               \textcolor{keywordtype}{intent(out)}:: det\textcolor{comment}{   !< Determinant}
1152   \textcolor{keywordtype}{real},               \textcolor{keywordtype}{intent(out)}:: ddet\textcolor{comment}{  !< Derivative of determinant with lam}
1153   \textcolor{keywordtype}{real},     \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: row\_scale\textcolor{comment}{ !< A scaling factor of the rows of the}
1154 \textcolor{comment}{                                      !! matrix to limit the growth of the determinant}
1155   \textcolor{comment}{! Local variables}
1156   \textcolor{keywordtype}{real} :: detkm1, detkm2   \textcolor{comment}{! Cumulative value of the determinant for the previous two layers.}
1157   \textcolor{keywordtype}{real} :: ddetkm1, ddetkm2 \textcolor{comment}{! Derivative of the cumulative determinant with lam for the previous two layers.}
1158   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: rescale = 1024.0**4 \textcolor{comment}{! max value of determinant allowed before rescaling}
1159   \textcolor{keywordtype}{real} :: rscl      \textcolor{comment}{! A rescaling factor that is applied succesively to each row.}
1160   \textcolor{keywordtype}{real} :: i\_rescale \textcolor{comment}{! inverse of rescale}
1161   \textcolor{keywordtype}{integer} :: k      \textcolor{comment}{! row (layer interface) index}
1162 
1163   i\_rescale = 1.0 / rescale
1164   rscl = 1.0 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(row\_scale)) rscl = row\_scale
1165 
1166   detkm1 = 1.0 ; ddetkm1 = 0.0
1167   det = (a(ks)+c(ks)) - lam ; ddet = -1.0
1168   \textcolor{keywordflow}{do} k=ks+1,ke
1169     \textcolor{comment}{! Shift variables and rescale rows to avoid over- or underflow.}
1170     detkm2 = row\_scale*detkm1 ; ddetkm2 = row\_scale*ddetkm1
1171     detkm1 = row\_scale*det    ; ddetkm1 = row\_scale*ddet
1172 
1173     det =  ((a(k)+c(k))-lam)*detkm1  - (a(k)*c(k-1))*detkm2
1174     ddet = ((a(k)+c(k))-lam)*ddetkm1 - (a(k)*c(k-1))*ddetkm2 - detkm1
1175 
1176     \textcolor{comment}{! Rescale det & ddet if det is getting too large or too small.}
1177     \textcolor{keywordflow}{if} (abs(det) > rescale) \textcolor{keywordflow}{then}
1178       det = i\_rescale*det ; detkm1 = i\_rescale*detkm1
1179       ddet = i\_rescale*ddet ; ddetkm1 = i\_rescale*ddetkm1
1180     \textcolor{keywordflow}{elseif} (abs(det) < i\_rescale) \textcolor{keywordflow}{then}
1181       det = rescale*det ; detkm1 = rescale*detkm1
1182       ddet = rescale*ddet ; ddetkm1 = rescale*ddetkm1
1183 \textcolor{keywordflow}{    endif}
1184 \textcolor{keywordflow}{  enddo}
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}}
\subsubsection{\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(\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{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & Thermodynamic variables\\
\hline
\mbox{\tt 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{\tt in}  & {\em full\+\_\+halos} & If true, do the calculation over the entire computational domain.\\
\hline
\mbox{\tt in}  & {\em use\+\_\+ebt\+\_\+mode} & If true, use the equivalent barotropic mode instead of the first baroclinic mode.\\
\hline
\mbox{\tt 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{\tt 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{\tt out}  & {\em modal\+\_\+structure} & Normalized model structure \mbox{[}nondim\mbox{]}\\
\hline
\mbox{\tt 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{\tt 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{\tt 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}
59   \textcolor{keywordtype}{type}(ocean\_grid\_type),            \textcolor{keywordtype}{intent(in)}  :: g\textcolor{comment}{  !< Ocean grid structure}
60   \textcolor{keywordtype}{type}(verticalgrid\_type),          \textcolor{keywordtype}{intent(in)}  :: gv\textcolor{comment}{ !< Vertical grid structure}
61   \textcolor{keywordtype}{type}(unit\_scale\_type),            \textcolor{keywordtype}{intent(in)}  :: us\textcolor{comment}{ !< A dimensional unit scaling type}
62   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
63                                     \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{  !< Layer thickness [H ~> m or kg m-2]}
64   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),            \textcolor{keywordtype}{intent(in)}  :: tv\textcolor{comment}{ !< Thermodynamic variables}
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]}
66   \textcolor{keywordtype}{type}(wave\_speed\_cs),              \textcolor{keywordtype}{pointer}     :: cs\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}
67   \textcolor{keywordtype}{logical},                \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: full\_halos\textcolor{comment}{ !< If true, do the calculation}
68 \textcolor{comment}{                                          !! over the entire computational domain.}
69   \textcolor{keywordtype}{logical},                \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: use\_ebt\_mode\textcolor{comment}{ !< If true, use the equivalent}
70 \textcolor{comment}{                                          !! barotropic mode instead of the first baroclinic mode.}
71   \textcolor{keywordtype}{real},                   \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: mono\_n2\_column\_fraction\textcolor{comment}{ !< The lower fraction}
72 \textcolor{comment}{                                          !! of water column over which N2 is limited as monotonic}
73 \textcolor{comment}{                                          !! for the purposes of calculating vertical modal structure.}
74   \textcolor{keywordtype}{real},                   \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: mono\_n2\_depth\textcolor{comment}{ !< A depth below which N2 is limited as}
75 \textcolor{comment}{                                          !! monotonic for the purposes of calculating vertical}
76 \textcolor{comment}{                                          !! modal structure [Z ~> m].}
77   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
78                           \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: modal\_structure\textcolor{comment}{ !< Normalized model structure [nondim]}
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}
80 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}
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}
82 \textcolor{comment}{                                     !! below which 0 is returned [L T-1 ~> m s-1].}
83   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}
84 \textcolor{comment}{                                     !! wave speeds [nondim]}
85 
86   \textcolor{comment}{! Local variables}
87   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: &
88     drho\_dt, &    \textcolor{comment}{! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]}
89     drho\_ds, &    \textcolor{comment}{! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]}
90     pres, &       \textcolor{comment}{! Interface pressure [R L2 T-2 ~> Pa]}
91     t\_int, &      \textcolor{comment}{! Temperature interpolated to interfaces [degC]}
92     s\_int, &      \textcolor{comment}{! Salinity interpolated to interfaces [ppt]}
93     h\_top, &      \textcolor{comment}{! The distance of each filtered interface from the ocean surface [Z ~> m]}
94     h\_bot, &      \textcolor{comment}{! The distance of each filtered interface from the bottom [Z ~> m]}
95     gprime        \textcolor{comment}{! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].}
96   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: &
97     igl, igu      \textcolor{comment}{! The inverse of the reduced gravity across an interface times}
98                   \textcolor{comment}{! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].}
99   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G),SZI\_(G))} :: &
100     hf, &         \textcolor{comment}{! Layer thicknesses after very thin layers are combined [Z ~> m]}
101     tf, &         \textcolor{comment}{! Layer temperatures after very thin layers are combined [degC]}
102     sf, &         \textcolor{comment}{! Layer salinities after very thin layers are combined [ppt]}
103     rf            \textcolor{comment}{! Layer densities after very thin layers are combined [R ~> kg m-3]}
104   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: &
105     hc, &         \textcolor{comment}{! A column of layer thicknesses after convective istabilities are removed [Z ~> m]}
106     tc, &         \textcolor{comment}{! A column of layer temperatures after convective istabilities are removed [degC]}
107     sc, &         \textcolor{comment}{! A column of layer salinites after convective istabilities are removed [ppt]}
108     rc, &         \textcolor{comment}{! A column of layer densities after convective istabilities are removed [R ~> kg m-3]}
109     hc\_h          \textcolor{comment}{! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2]}
110   \textcolor{keywordtype}{real} :: i\_htot  \textcolor{comment}{! The inverse of the total filtered thicknesses [Z ~> m]}
111   \textcolor{keywordtype}{real} :: det, ddet, detkm1, detkm2, ddetkm1, ddetkm2
112   \textcolor{keywordtype}{real} :: lam     \textcolor{comment}{! The eigenvalue [T2 L-2 ~> s m-1]}
113   \textcolor{keywordtype}{real} :: dlam    \textcolor{comment}{! The change in estimates of the eigenvalue [T2 L-2 ~> s m-1]}
114   \textcolor{keywordtype}{real} :: lam0    \textcolor{comment}{! The first guess of the eigenvalue [T2 L-2 ~> s m-1]}
115   \textcolor{keywordtype}{real} :: min\_h\_frac \textcolor{comment}{! [nondim]}
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]}
117   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
118     htot, hmin, &  \textcolor{comment}{! Thicknesses [Z ~> m]}
119     h\_here, &      \textcolor{comment}{! A thickness [Z ~> m]}
120     hxt\_here, &    \textcolor{comment}{! A layer integrated temperature [degC Z ~> degC m]}
121     hxs\_here, &    \textcolor{comment}{! A layer integrated salinity [ppt Z ~> ppt m]}
122     hxr\_here       \textcolor{comment}{! A layer integrated density [R Z ~> kg m-2]}
123   \textcolor{keywordtype}{real} :: speed2\_tot \textcolor{comment}{! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]}
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]}
125   \textcolor{keywordtype}{real} :: i\_hnew   \textcolor{comment}{! The inverse of a new layer thickness [Z-1 ~> m-1]}
126   \textcolor{keywordtype}{real} :: drxh\_sum \textcolor{comment}{! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2]}
127   \textcolor{keywordtype}{real} :: l2\_to\_z2 \textcolor{comment}{! A scaling factor squared from units of lateral distances to depths [Z2 L-2 ~> 1].}
128   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{pointer}, \textcolor{keywordtype}{dimension(:,:,:)} :: t => null(), s => null()
129   \textcolor{keywordtype}{real} :: g\_rho0   \textcolor{comment}{! G\_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].}
130   \textcolor{keywordtype}{real} :: c2\_scale \textcolor{comment}{! A scaling factor for wave speeds to help control the growth of the determinant}
131                    \textcolor{comment}{! and its derivative with lam between rows of the Thomas algorithm solver.  The}
132                    \textcolor{comment}{! exact value should not matter for the final result if it is an even power of 2.}
133   \textcolor{keywordtype}{real} :: tol\_hfrac \textcolor{comment}{! Layers that together are smaller than this fraction of}
134                     \textcolor{comment}{! the total water column can be merged for efficiency.}
135   \textcolor{keywordtype}{real} :: tol\_solve \textcolor{comment}{! The fractional tolerance with which to solve for the wave speeds [nondim]}
136   \textcolor{keywordtype}{real} :: tol\_merge \textcolor{comment}{! The fractional change in estimated wave speed that is allowed}
137                     \textcolor{comment}{! when deciding to merge layers in the calculation [nondim]}
138   \textcolor{keywordtype}{real} :: rescale, i\_rescale
139   \textcolor{keywordtype}{integer} :: kf(szi\_(g)) \textcolor{comment}{! The number of active layers after filtering.}
140   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: max\_itt = 10
141   \textcolor{keywordtype}{real} :: lam\_it(max\_itt), det\_it(max\_itt), ddet\_it(max\_itt)
142   \textcolor{keywordtype}{logical} :: use\_eos    \textcolor{comment}{! If true, density is calculated from T & S using an equation of state.}
143   \textcolor{keywordtype}{logical} :: better\_est \textcolor{comment}{! If true, use an improved estimate of the first mode internal wave speed.}
144   \textcolor{keywordtype}{logical} :: merge      \textcolor{comment}{! If true, merge the current layer with the one above.}
145   \textcolor{keywordtype}{integer} :: kc         \textcolor{comment}{! The number of layers in the column after merging}
146   \textcolor{keywordtype}{integer} :: i, j, k, k2, itt, is, ie, js, je, nz
147   \textcolor{keywordtype}{real} :: hw, sum\_hc
148   \textcolor{keywordtype}{real} :: gp      \textcolor{comment}{! A limited local copy of gprime [L2 Z-1 T-2 ~> m s-2]}
149   \textcolor{keywordtype}{real} :: n2min   \textcolor{comment}{! A minimum buoyancy frequency [T-2 ~> s-2]}
150   \textcolor{keywordtype}{logical} :: l\_use\_ebt\_mode, calc\_modal\_structure
151   \textcolor{keywordtype}{real} :: l\_mono\_n2\_column\_fraction, l\_mono\_n2\_depth
152   \textcolor{keywordtype}{real} :: mode\_struct(szk\_(g)), ms\_min, ms\_max, ms\_sq
153 
154   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
155 
156   \textcolor{keywordflow}{if} (.not. \textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_wave\_speed: "}// &
157            \textcolor{stringliteral}{"Module must be initialized before it is used."})
158   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(full\_halos)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (full\_halos) \textcolor{keywordflow}{then}
159     is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
160 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
161 
162   l2\_to\_z2 = us%L\_to\_Z**2
163 
164   l\_use\_ebt\_mode = cs%use\_ebt\_mode
165   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(use\_ebt\_mode)) l\_use\_ebt\_mode = use\_ebt\_mode
166   l\_mono\_n2\_column\_fraction = cs%mono\_N2\_column\_fraction
167   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_column\_fraction)) l\_mono\_n2\_column\_fraction = mono\_n2\_column\_fraction
168   l\_mono\_n2\_depth = cs%mono\_N2\_depth
169   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_depth)) l\_mono\_n2\_depth = mono\_n2\_depth
170   calc\_modal\_structure = l\_use\_ebt\_mode
171   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) calc\_modal\_structure = .true.
172   \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}
173     \textcolor{keywordflow}{do} k=1,nz; \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is,ie
174       modal\_structure(i,j,k) = 0.0
175 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
176 \textcolor{keywordflow}{  endif}
177 
178   s => tv%S ; t => tv%T
179   g\_rho0 = gv%g\_Earth / gv%Rho0
180   \textcolor{comment}{! Simplifying the following could change answers at roundoff.}
181   z\_to\_pres = gv%Z\_to\_H * (gv%H\_to\_RZ * gv%g\_Earth)
182   use\_eos = \textcolor{keyword}{associated}(tv%eqn\_of\_state)
183 
184   better\_est = cs%better\_cg1\_est ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(better\_speed\_est)) better\_est = better\_speed\_est
185 
186   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
187     tol\_solve = cs%wave\_speed\_tol ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(wave\_speed\_tol)) tol\_solve = wave\_speed\_tol
188     tol\_hfrac  = 0.1*tol\_solve ; tol\_merge = tol\_solve / \textcolor{keywordtype}{real}(nz)
189   \textcolor{keywordflow}{else}
190     tol\_solve = 0.001 ; tol\_hfrac  = 0.0001 ; tol\_merge = 0.001
191 \textcolor{keywordflow}{  endif}
192 
193   \textcolor{comment}{! The rescaling below can control the growth of the determinant provided that}
194   \textcolor{comment}{! (tol\_merge*cg1\_min2/c2\_scale > I\_rescale).  For default values, this suggests a stable lower}
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.}
196   \textcolor{comment}{! The upper bound on the rate of increase in the determinant is g'H/c2\_scale < rescale or in the}
197   \textcolor{comment}{! worst possible oceanic case of g'H < 0.5*10m/s2*1e4m = 5.e4 m2/s2 < 1024**2*c2\_scale, suggesting}
198   \textcolor{comment}{! that c2\_scale can safely be set to 1/(16*1024**2), which would decrease the stable floor on}
199   \textcolor{comment}{! min\_speed to ~6.9e-8 m/s for 90 layers or 2.33e-7 m/s for 1000 layers.}
200   cg1\_min2 = cs%min\_speed2 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(min\_speed)) cg1\_min2 = min\_speed**2
201   rescale = 1024.0**4 ; i\_rescale = 1.0/rescale
202   c2\_scale = us%m\_s\_to\_L\_T**2 / 4096.0**2 \textcolor{comment}{! Other powers of 2 give identical results.}
203 
204   min\_h\_frac = tol\_hfrac / \textcolor{keywordtype}{real}(nz)
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,&}
206 \textcolor{comment}{!$OMP                                  calc\_modal\_structure,l\_use\_ebt\_mode,modal\_structure, &}
207 \textcolor{comment}{!$OMP                                  l\_mono\_N2\_column\_fraction,l\_mono\_N2\_depth,CS,   &}
208 \textcolor{comment}{!$OMP                                  Z\_to\_pres,cg1,g\_Rho0,rescale,I\_rescale,L2\_to\_Z2, &}
209 \textcolor{comment}{!$OMP                                  better\_est,cg1\_min2,tol\_merge,tol\_solve,c2\_scale) &}
210 \textcolor{comment}{!$OMP                          private(htot,hmin,kf,H\_here,HxT\_here,HxS\_here,HxR\_here, &}
211 \textcolor{comment}{!$OMP                                  Hf,Tf,Sf,Rf,pres,T\_int,S\_int,drho\_dT,drho\_dS,   &}
212 \textcolor{comment}{!$OMP                                  drxh\_sum,kc,Hc,Hc\_H,tC,sc,I\_Hnew,gprime,&}
213 \textcolor{comment}{!$OMP                                  Rc,speed2\_tot,Igl,Igu,lam0,lam,lam\_it,dlam, &}
214 \textcolor{comment}{!$OMP                                  mode\_struct,sum\_hc,N2min,gp,hw,                 &}
215 \textcolor{comment}{!$OMP                                  ms\_min,ms\_max,ms\_sq,H\_top,H\_bot,I\_Htot,merge,   &}
216 \textcolor{comment}{!$OMP                                  det,ddet,detKm1,ddetKm1,detKm2,ddetKm2,det\_it,ddet\_it)}
217   \textcolor{keywordflow}{do} j=js,je
218     \textcolor{comment}{!   First merge very thin layers with the one above (or below if they are}
219     \textcolor{comment}{! at the top).  This also transposes the row order so that columns can}
220     \textcolor{comment}{! be worked upon one at a time.}
221     \textcolor{keywordflow}{do} i=is,ie ; htot(i) = 0.0 ;\textcolor{keywordflow}{ enddo}
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}
223 
224     \textcolor{keywordflow}{do} i=is,ie
225       hmin(i) = htot(i)*min\_h\_frac ; kf(i) = 1 ; h\_here(i) = 0.0
226       hxt\_here(i) = 0.0 ; hxs\_here(i) = 0.0 ; hxr\_here(i) = 0.0
227 \textcolor{keywordflow}{    enddo}
228     \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
229       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
230         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}
231           hf(kf(i),i) = h\_here(i)
232           tf(kf(i),i) = hxt\_here(i) / h\_here(i)
233           sf(kf(i),i) = hxs\_here(i) / h\_here(i)
234           kf(i) = kf(i) + 1
235 
236           \textcolor{comment}{! Start a new layer}
237           h\_here(i) = h(i,j,k)*gv%H\_to\_Z
238           hxt\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*t(i,j,k)
239           hxs\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*s(i,j,k)
240         \textcolor{keywordflow}{else}
241           h\_here(i) = h\_here(i) + h(i,j,k)*gv%H\_to\_Z
242           hxt\_here(i) = hxt\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*t(i,j,k)
243           hxs\_here(i) = hxs\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*s(i,j,k)
244 \textcolor{keywordflow}{        endif}
245 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
246       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}
247         hf(kf(i),i) = h\_here(i)
248         tf(kf(i),i) = hxt\_here(i) / h\_here(i)
249         sf(kf(i),i) = hxs\_here(i) / h\_here(i)
250 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
251     \textcolor{keywordflow}{else}
252       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
253         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}
254           hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)
255           kf(i) = kf(i) + 1
256 
257           \textcolor{comment}{! Start a new layer}
258           h\_here(i) = h(i,j,k)*gv%H\_to\_Z
259           hxr\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*gv%Rlay(k)
260         \textcolor{keywordflow}{else}
261           h\_here(i) = h\_here(i) + h(i,j,k)*gv%H\_to\_Z
262           hxr\_here(i) = hxr\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*gv%Rlay(k)
263 \textcolor{keywordflow}{        endif}
264 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
265       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}
266         hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)
267 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
268 \textcolor{keywordflow}{    endif}
269 
270     \textcolor{comment}{! From this point, we can work on individual columns without causing memory to have page faults.}
271     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then}
272       \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
273         pres(1) = 0.0 ; h\_top(1) = 0.0
274         \textcolor{keywordflow}{do} k=2,kf(i)
275           pres(k) = pres(k-1) + z\_to\_pres*hf(k-1,i)
276           t\_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
277           s\_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
278           h\_top(k) = h\_top(k-1) + hf(k-1,i)
279 \textcolor{keywordflow}{        enddo}
280         \textcolor{keyword}{call }calculate\_density\_derivs(t\_int, s\_int, pres, drho\_dt, drho\_ds, &
281                                       tv%eqn\_of\_state, (/2,kf(i)/) )
282 
283         \textcolor{comment}{! Sum the reduced gravities to find out how small a density difference is negligibly small.}
284         drxh\_sum = 0.0
285         \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
286           \textcolor{comment}{! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for}
287           \textcolor{comment}{! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.}
288           \textcolor{comment}{! For a uniform stratification and a huge number of layers uniformly distributed in}
289           \textcolor{comment}{! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.}
290           \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}
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.}
292             h\_bot(kf(i)+1) = 0.0
293             \textcolor{keywordflow}{do} k=kf(i),2,-1
294               h\_bot(k) = h\_bot(k+1) + hf(k,i)
295               drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * &
296                   max(0.0, drho\_dt(k)*(tf(k,i)-tf(k-1,i)) + drho\_ds(k)*(sf(k,i)-sf(k-1,i)))
297 \textcolor{keywordflow}{            enddo}
298 \textcolor{keywordflow}{          endif}
299         \textcolor{keywordflow}{else}
300           \textcolor{comment}{! This estimate is problematic in that it goes like 1/nz for a large number of layers,}
301           \textcolor{comment}{! but it is an overestimate (as desired) for a small number of layers, by at a factor}
302           \textcolor{comment}{! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.}
303           \textcolor{keywordflow}{do} k=2,kf(i)
304             drxh\_sum = drxh\_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
305                 max(0.0, drho\_dt(k)*(tf(k,i)-tf(k-1,i)) + drho\_ds(k)*(sf(k,i)-sf(k-1,i)))
306 \textcolor{keywordflow}{          enddo}
307 \textcolor{keywordflow}{        endif}
308       \textcolor{keywordflow}{else}
309         drxh\_sum = 0.0
310         \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
311           h\_top(1) = 0.0
312           \textcolor{keywordflow}{do} k=2,kf(i) ; h\_top(k) = h\_top(k-1) + hf(k-1,i) ;\textcolor{keywordflow}{ enddo}
313           \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}
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.}
315             h\_bot(kf(i)+1) = 0.0
316             \textcolor{keywordflow}{do} k=kf(i),2,-1
317               h\_bot(k) = h\_bot(k+1) + hf(k,i)
318               drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * max(0.0,rf(k,i)-rf(k-1,i))
319 \textcolor{keywordflow}{            enddo}
320 \textcolor{keywordflow}{          endif}
321         \textcolor{keywordflow}{else}
322           \textcolor{keywordflow}{do} k=2,kf(i)
323             drxh\_sum = drxh\_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
324 \textcolor{keywordflow}{          enddo}
325 \textcolor{keywordflow}{        endif}
326 \textcolor{keywordflow}{      endif}
327 
328       \textcolor{comment}{!   Find gprime across each internal interface, taking care of convective instabilities by}
329       \textcolor{comment}{! merging layers.  If the estimated wave speed is too small, simply return zero.}
330       \textcolor{keywordflow}{if} (g\_rho0 * drxh\_sum <= cg1\_min2) \textcolor{keywordflow}{then}
331         cg1(i,j) = 0.0
332         \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) modal\_structure(i,j,:) = 0.
333       \textcolor{keywordflow}{else}
334         \textcolor{comment}{! Merge layers to eliminate convective instabilities or exceedingly}
335         \textcolor{comment}{! small reduced gravities.  Merging layers reduces the estimated wave speed by}
336         \textcolor{comment}{! (rho(2)-rho(1))*h(1)*h(2) / H\_tot.}
337         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
338           kc = 1
339           hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
340           \textcolor{keywordflow}{do} k=2,kf(i)
341             \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
342               merge = ((drho\_dt(k)*(tf(k,i)-tc(kc)) + drho\_ds(k)*(sf(k,i)-sc(kc))) * &
343                        ((hc(kc) * hf(k,i))*i\_htot) < 2.0 * tol\_merge*drxh\_sum)
344             \textcolor{keywordflow}{else}
345               merge = ((drho\_dt(k)*(tf(k,i)-tc(kc)) + drho\_ds(k)*(sf(k,i)-sc(kc))) * &
346                        (hc(kc) + hf(k,i)) < 2.0 * tol\_merge*drxh\_sum)
347 \textcolor{keywordflow}{            endif}
348             \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
349               \textcolor{comment}{! Merge this layer with the one above and backtrack.}
350               i\_hnew = 1.0 / (hc(kc) + hf(k,i))
351               tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i\_hnew
352               sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i\_hnew
353               hc(kc) = (hc(kc) + hf(k,i))
354               \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}
355               \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}
356               \textcolor{comment}{! far back we go.}
357               \textcolor{keywordflow}{do} k2=kc,2,-1
358                 \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
359                   merge = ((drho\_dt(k2)*(tc(k2)-tc(k2-1)) + drho\_ds(k2)*(sc(k2)-sc(k2-1))) * &
360                            ((hc(k2) * hc(k2-1))*i\_htot) < tol\_merge*drxh\_sum)
361                 \textcolor{keywordflow}{else}
362                   merge = ((drho\_dt(k2)*(tc(k2)-tc(k2-1)) + drho\_ds(k2)*(sc(k2)-sc(k2-1))) * &
363                            (hc(k2) + hc(k2-1)) < tol\_merge*drxh\_sum)
364 \textcolor{keywordflow}{                endif}
365                 \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
366                   \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}
367                   i\_hnew = 1.0 / (hc(kc) + hc(kc-1))
368                   tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i\_hnew
369                   sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i\_hnew
370                   hc(kc-1) = (hc(kc) + hc(kc-1))
371                   kc = kc - 1
372                 \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}
373 \textcolor{keywordflow}{              enddo}
374             \textcolor{keywordflow}{else}
375               \textcolor{comment}{! Add a new layer to the column.}
376               kc = kc + 1
377               drho\_ds(kc) = drho\_ds(k) ; drho\_dt(kc) = drho\_dt(k)
378               tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
379 \textcolor{keywordflow}{            endif}
380 \textcolor{keywordflow}{          enddo}
381           \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}
382           \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-Boussinesq.}
383             gprime(k) = g\_rho0 * (drho\_dt(k)*(tc(k)-tc(k-1)) + drho\_ds(k)*(sc(k)-sc(k-1)))
384 \textcolor{keywordflow}{          enddo}
385         \textcolor{keywordflow}{else}  \textcolor{comment}{! .not.use\_EOS}
386           \textcolor{comment}{! Do the same with density directly...}
387           kc = 1
388           hc(1) = hf(1,i) ; rc(1) = rf(1,i)
389           \textcolor{keywordflow}{do} k=2,kf(i)
390             \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
391               merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i\_htot) < 2.0*tol\_merge*drxh\_sum)
392             \textcolor{keywordflow}{else}
393               merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol\_merge*drxh\_sum)
394 \textcolor{keywordflow}{            endif}
395             \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
396               \textcolor{comment}{! Merge this layer with the one above and backtrack.}
397               rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
398               hc(kc) = (hc(kc) + hf(k,i))
399               \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}
400               \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}
401               \textcolor{comment}{! far back we go.}
402               \textcolor{keywordflow}{do} k2=kc,2,-1
403                 \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
404                   merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i\_htot) < tol\_merge*drxh\_sum)
405                 \textcolor{keywordflow}{else}
406                   merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol\_merge*drxh\_sum)
407 \textcolor{keywordflow}{                endif}
408                 \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
409                   \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}
410                   rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
411                   hc(kc-1) = (hc(kc) + hc(kc-1))
412                   kc = kc - 1
413                 \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}
414 \textcolor{keywordflow}{              enddo}
415             \textcolor{keywordflow}{else}
416               \textcolor{comment}{! Add a new layer to the column.}
417               kc = kc + 1
418               rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
419 \textcolor{keywordflow}{            endif}
420 \textcolor{keywordflow}{          enddo}
421           \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}
422           \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-Boussinesq.}
423             gprime(k) = g\_rho0 * (rc(k)-rc(k-1))
424 \textcolor{keywordflow}{          enddo}
425 \textcolor{keywordflow}{        endif}  \textcolor{comment}{! use\_EOS}
426 
427         \textcolor{comment}{! Sum the contributions from all of the interfaces to give an over-estimate}
428         \textcolor{comment}{! of the first-mode wave speed.  Also populate Igl and Igu which are the}
429         \textcolor{comment}{! non-leading diagonals of the tridiagonal matrix.}
430         \textcolor{keywordflow}{if} (kc >= 2) \textcolor{keywordflow}{then}
431           speed2\_tot = 0.0
432           \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
433             h\_top(1) = 0.0 ; h\_bot(kc+1) = 0.0
434             \textcolor{keywordflow}{do} k=2,kc+1 ; h\_top(k) = h\_top(k-1) + hc(k-1) ;\textcolor{keywordflow}{ enddo}
435             \textcolor{keywordflow}{do} k=kc,2,-1 ; h\_bot(k) = h\_bot(k+1) + hc(k) ;\textcolor{keywordflow}{ enddo}
436             i\_htot = 0.0 ; \textcolor{keywordflow}{if} (h\_top(kc+1) > 0.0) i\_htot = 1.0 / h\_top(kc+1)
437 \textcolor{keywordflow}{          endif}
438 
439           \textcolor{keywordflow}{if} (l\_use\_ebt\_mode) \textcolor{keywordflow}{then}
440             igu(1) = 0. \textcolor{comment}{! Neumann condition for pressure modes}
441             sum\_hc = hc(1)
442             n2min = l2\_to\_z2*gprime(2)/hc(1)
443             \textcolor{keywordflow}{do} k=2,kc
444               hw = 0.5*(hc(k-1)+hc(k))
445               gp = gprime(k)
446               \textcolor{keywordflow}{if} (l\_mono\_n2\_column\_fraction>0. .or. l\_mono\_n2\_depth>=0.) \textcolor{keywordflow}{then}
447                 \textcolor{keywordflow}{if} ( ((g%bathyT(i,j)-sum\_hc < l\_mono\_n2\_column\_fraction*g%bathyT(i,j)) .or. &
448                       ((l\_mono\_n2\_depth >= 0.) .and. (sum\_hc > l\_mono\_n2\_depth))) .and. &
449                      (l2\_to\_z2*gp > n2min*hw) ) \textcolor{keywordflow}{then}
450                   \textcolor{comment}{! Filters out regions where N2 increases with depth but only in a lower fraction}
451                   \textcolor{comment}{! of the water column or below a certain depth.}
452                   gp = us%Z\_to\_L**2 * (n2min*hw)
453                 \textcolor{keywordflow}{else}
454                   n2min = l2\_to\_z2 * gp/hw
455 \textcolor{keywordflow}{                endif}
456 \textcolor{keywordflow}{              endif}
457               igu(k) = 1.0/(gp*hc(k))
458               igl(k-1) = 1.0/(gp*hc(k-1))
459               sum\_hc = sum\_hc + hc(k)
460               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
461                 \textcolor{comment}{! Estimate that the ebt\_mode is sqrt(2) times the speed of the flat bottom modes.}
462                 speed2\_tot = speed2\_tot + 2.0 * gprime(k)*((h\_top(k) * h\_bot(k)) * i\_htot)
463               \textcolor{keywordflow}{else} \textcolor{comment}{! The ebt\_mode wave should be faster than the flat-bottom mode, so 0.707 should be > 1?}
464                 speed2\_tot = speed2\_tot + gprime(k)*(hc(k-1)+hc(k))*0.707
465 \textcolor{keywordflow}{              endif}
466 \textcolor{keywordflow}{            enddo}
467            \textcolor{comment}{!Igl(kc) = 0. ! Neumann condition for pressure modes}
468             igl(kc) = 2.*igu(kc) \textcolor{comment}{! Dirichlet condition for pressure modes}
469           \textcolor{keywordflow}{else} \textcolor{comment}{! .not. l\_use\_ebt\_mode}
470             \textcolor{keywordflow}{do} k=2,kc
471               igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
472               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
473                 speed2\_tot = speed2\_tot + gprime(k)*((h\_top(k) * h\_bot(k)) * i\_htot)
474               \textcolor{keywordflow}{else}
475                 speed2\_tot = speed2\_tot + gprime(k)*(hc(k-1)+hc(k))
476 \textcolor{keywordflow}{              endif}
477 \textcolor{keywordflow}{            enddo}
478 \textcolor{keywordflow}{          endif}
479 
480           \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}
481             mode\_struct(:) = 0.
482             mode\_struct(1:kc) = 1. \textcolor{comment}{! Uniform flow, first guess}
483 \textcolor{keywordflow}{          endif}
484 
485           \textcolor{comment}{! Under estimate the first eigenvalue (overestimate the speed) to start with.}
486           \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}
487             lam0 = 0.5 / speed2\_tot ; lam = lam0
488           \textcolor{keywordflow}{else}
489             lam0 = 1.0 / speed2\_tot ; lam = lam0
490 \textcolor{keywordflow}{          endif}
491           \textcolor{comment}{! Find the determinant and its derivative with lam.}
492           \textcolor{keywordflow}{do} itt=1,max\_itt
493             lam\_it(itt) = lam
494             \textcolor{keywordflow}{if} (l\_use\_ebt\_mode) \textcolor{keywordflow}{then}
495               \textcolor{comment}{! This initialization of det,ddet imply Neumann boundary conditions for horizontal}
496               \textcolor{comment}{! velocity or pressure modes, so that first 3 rows of the matrix are}
497               \textcolor{comment}{!    /   b(1)-lam  igl(1)      0        0     0  ...  \(\backslash\)}
498               \textcolor{comment}{!    |  igu(2)    b(2)-lam   igl(2)     0     0  ...  |}
499               \textcolor{comment}{!    |    0        igu(3)   b(3)-lam  igl(3)  0  ...  |}
500               \textcolor{comment}{! The last two rows of the pressure equation matrix are}
501               \textcolor{comment}{!    |    ...  0  igu(kc-1)  b(kc-1)-lam  igl(kc-1)  |}
502               \textcolor{comment}{!    \(\backslash\)    ...  0     0        igu(kc)     b(kc)-lam  /}
503               \textcolor{keyword}{call }tridiag\_det(igu, igl, 1, kc, lam, det, ddet, row\_scale=c2\_scale)
504             \textcolor{keywordflow}{else}
505               \textcolor{comment}{! This initialization of det,ddet imply Dirichlet boundary conditions for vertical}
506               \textcolor{comment}{! velocity modes, so that first 3 rows of the matrix are}
507               \textcolor{comment}{!    /  b(2)-lam  igl(2)      0       0     0  ...  |}
508               \textcolor{comment}{!    |  igu(3)  b(3)-lam   igl(3)     0     0  ...  |}
509               \textcolor{comment}{!    |    0       igu(4)  b(4)-lam  igl(4)  0  ...  |}
510               \textcolor{comment}{! The last three rows of the w equation matrix are}
511               \textcolor{comment}{!    |    ...   0  igu(kc-2)  b(kc-2)-lam  igl(kc-2)     0       |}
512               \textcolor{comment}{!    |    ...   0     0        igu(kc-1)  b(kc-1)-lam  igl(kc-1) |}
513               \textcolor{comment}{!    \(\backslash\)    ...   0     0           0        igu(kc)    b(kc)-lam  /}
514               \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lam, det, ddet, row\_scale=c2\_scale)
515 \textcolor{keywordflow}{            endif}
516             \textcolor{comment}{! Use Newton's method iteration to find a new estimate of lam.}
517             det\_it(itt) = det ; ddet\_it(itt) = ddet
518 
519             \textcolor{keywordflow}{if} ((ddet >= 0.0) .or. (-det > -0.5*lam*ddet)) \textcolor{keywordflow}{then}
520               \textcolor{comment}{! lam was not an under-estimate, as intended, so Newton's method}
521               \textcolor{comment}{! may not be reliable; lam must be reduced, but not by more}
522               \textcolor{comment}{! than half.}
523               lam = 0.5 * lam
524               dlam = -lam
525             \textcolor{keywordflow}{else}  \textcolor{comment}{! Newton's method is OK.}
526               dlam = - det / ddet
527               lam = lam + dlam
528 \textcolor{keywordflow}{            endif}
529 
530             \textcolor{keywordflow}{if} (calc\_modal\_structure) \textcolor{keywordflow}{then}
531               \textcolor{keyword}{call }tdma6(kc, igu, igl, lam, mode\_struct)
532               ms\_min = mode\_struct(1)
533               ms\_max = mode\_struct(1)
534               ms\_sq = mode\_struct(1)**2
535               \textcolor{keywordflow}{do} k = 2,kc
536                 ms\_min = min(ms\_min, mode\_struct(k))
537                 ms\_max = max(ms\_max, mode\_struct(k))
538                 ms\_sq = ms\_sq + mode\_struct(k)**2
539 \textcolor{keywordflow}{              enddo}
540               \textcolor{keywordflow}{if} (ms\_min<0. .and. ms\_max>0.) \textcolor{keywordflow}{then} \textcolor{comment}{! Any zero crossings => lam is too high}
541                 lam = 0.5 * ( lam - dlam )
542                 dlam = -lam
543                 mode\_struct(1:kc) = abs(mode\_struct(1:kc)) / sqrt( ms\_sq )
544               \textcolor{keywordflow}{else}
545                 mode\_struct(1:kc) = mode\_struct(1:kc) / sqrt( ms\_sq )
546 \textcolor{keywordflow}{              endif}
547 \textcolor{keywordflow}{            endif}
548 
549             \textcolor{keywordflow}{if} (abs(dlam) < tol\_solve*lam) \textcolor{keywordflow}{exit}
550 \textcolor{keywordflow}{          enddo}
551 
552           cg1(i,j) = 0.0
553           \textcolor{keywordflow}{if} (lam > 0.0) cg1(i,j) = 1.0 / sqrt(lam)
554 
555           \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) \textcolor{keywordflow}{then}
556             \textcolor{keywordflow}{if} (mode\_struct(1)/=0.) \textcolor{keywordflow}{then} \textcolor{comment}{! Normalize}
557               mode\_struct(1:kc) = mode\_struct(1:kc) / mode\_struct(1)
558             \textcolor{keywordflow}{else}
559               mode\_struct(1:kc)=0.
560 \textcolor{keywordflow}{            endif}
561             \textcolor{comment}{! Note that remapping\_core\_h requires that the same units be used}
562             \textcolor{comment}{! for both the source and target grid thicknesses, here [H ~> m or kg m-2].}
563             \textcolor{keywordflow}{do} k = 1,kc
564               hc\_h(k) = gv%Z\_to\_H * hc(k)
565 \textcolor{keywordflow}{            enddo}
566             \textcolor{keywordflow}{if} (cs%remap\_answers\_2018) \textcolor{keywordflow}{then}
567               \textcolor{keyword}{call }remapping\_core\_h(cs%remapping\_CS, kc, hc\_h(:), mode\_struct, &
568                                     nz, h(i,j,:), modal\_structure(i,j,:), &
569                                     1.0e-30*gv%m\_to\_H, 1.0e-10*gv%m\_to\_H)
570             \textcolor{keywordflow}{else}
571               \textcolor{keyword}{call }remapping\_core\_h(cs%remapping\_CS, kc, hc\_h(:), mode\_struct, &
572                                     nz, h(i,j,:), modal\_structure(i,j,:), &
573                                     gv%H\_subroundoff, gv%H\_subroundoff)
574 \textcolor{keywordflow}{            endif}
575 \textcolor{keywordflow}{          endif}
576         \textcolor{keywordflow}{else}
577           cg1(i,j) = 0.0
578           \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) modal\_structure(i,j,:) = 0.
579 \textcolor{keywordflow}{        endif}
580 \textcolor{keywordflow}{      endif} \textcolor{comment}{! cg1 /= 0.0}
581     \textcolor{keywordflow}{else}
582       cg1(i,j) = 0.0 \textcolor{comment}{! This is a land point.}
583       \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(modal\_structure)) modal\_structure(i,j,:) = 0.
584 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i-loop}
585 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-loop}
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}}
\subsubsection{\texorpdfstring{wave\+\_\+speed\+\_\+init()}{wave\_speed\_init()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+speed\+::wave\+\_\+speed\+\_\+init (\begin{DoxyParamCaption}\item[{type(\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{\tt in}  & {\em use\+\_\+ebt\+\_\+mode} & If true, use the equivalent barotropic mode instead of the first baroclinic mode.\\
\hline
\mbox{\tt 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{\tt 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{\tt 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{\tt 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{\tt 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{\tt 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}
1191   \textcolor{keywordtype}{type}(wave\_speed\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}
1192   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: use\_ebt\_mode\textcolor{comment}{  !< If true, use the equivalent}
1193 \textcolor{comment}{                                     !! barotropic mode instead of the first baroclinic mode.}
1194   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_n2\_column\_fraction\textcolor{comment}{ !< The lower fraction of water column over}
1195 \textcolor{comment}{                                     !! which N2 is limited as monotonic for the purposes of}
1196 \textcolor{comment}{                                     !! calculating the vertical modal structure.}
1197   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_n2\_depth\textcolor{comment}{ !< The depth below which N2 is limited}
1198 \textcolor{comment}{                                     !! as monotonic for the purposes of calculating the}
1199 \textcolor{comment}{                                     !! vertical modal structure [Z ~> m].}
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}
1201 \textcolor{comment}{                                     !! that recover the remapping answers from 2018.  Otherwise}
1202 \textcolor{comment}{                                     !! use more robust but mathematically equivalent expressions.}
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}
1204 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}
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}
1206 \textcolor{comment}{                                     !! below which 0 is returned [L T-1 ~> m s-1].}
1207   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}
1208 \textcolor{comment}{                                     !! wave speeds [nondim]}
1209 
1210   \textcolor{comment}{! This include declares and sets the variable "version".}
1211 \textcolor{preprocessor}{# include "version\_variable.h"}
1212 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_wave\_speed"}  \textcolor{comment}{! This module's name.}
1213 
1214   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}
1215     \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_speed\_init called with an "}// &
1216                             \textcolor{stringliteral}{"associated control structure."})
1217     \textcolor{keywordflow}{return}
1218   \textcolor{keywordflow}{else} ; \textcolor{keyword}{allocate}(cs) ;\textcolor{keywordflow}{ endif}
1219 
1220   \textcolor{comment}{! Write all relevant parameters to the model log.}
1221   \textcolor{keyword}{call }log\_version(mdl, version)
1222 
1223   \textcolor{keyword}{call }wave\_speed\_set\_param(cs, use\_ebt\_mode=use\_ebt\_mode, mono\_n2\_column\_fraction=mono\_n2\_column\_fraction,
       &
1224                             better\_speed\_est=better\_speed\_est, min\_speed=min\_speed, wave\_speed\_tol=
      wave\_speed\_tol)
1225 
1226   \textcolor{keyword}{call }initialize\_remapping(cs%remapping\_CS, \textcolor{stringliteral}{'PLM'}, boundary\_extrapolation=.false., &
1227                             answers\_2018=cs%remap\_answers\_2018)
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}}
\subsubsection{\texorpdfstring{wave\+\_\+speed\+\_\+set\+\_\+param()}{wave\_speed\_set\_param()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+wave\+\_\+speed\+::wave\+\_\+speed\+\_\+set\+\_\+param (\begin{DoxyParamCaption}\item[{type(\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{\tt in}  & {\em use\+\_\+ebt\+\_\+mode} & If true, use the equivalent barotropic mode instead of the first baroclinic mode.\\
\hline
\mbox{\tt 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{\tt 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{\tt 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{\tt 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{\tt 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{\tt 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}
1234   \textcolor{keywordtype}{type}(wave\_speed\_cs), \textcolor{keywordtype}{pointer}  :: cs\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}
1235   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: use\_ebt\_mode\textcolor{comment}{  !< If true, use the equivalent}
1236 \textcolor{comment}{                                      !! barotropic mode instead of the first baroclinic mode.}
1237   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_n2\_column\_fraction\textcolor{comment}{ !< The lower fraction of water column over}
1238 \textcolor{comment}{                                      !! which N2 is limited as monotonic for the purposes of}
1239 \textcolor{comment}{                                      !! calculating the vertical modal structure.}
1240   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: mono\_n2\_depth\textcolor{comment}{ !< The depth below which N2 is limited}
1241 \textcolor{comment}{                                      !! as monotonic for the purposes of calculating the}
1242 \textcolor{comment}{                                      !! vertical modal structure [Z ~> m].}
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}
1244 \textcolor{comment}{                                      !! that recover the remapping answers from 2018.  Otherwise}
1245 \textcolor{comment}{                                      !! use more robust but mathematically equivalent expressions.}
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}
1247 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}
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}
1249 \textcolor{comment}{                                     !! below which 0 is returned [L T-1 ~> m s-1].}
1250   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}
1251 \textcolor{comment}{                                     !! wave speeds [nondim]}
1252 
1253   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, &
1254      \textcolor{stringliteral}{"wave\_speed\_set\_param called with an associated control structure."})
1255 
1256   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(use\_ebt\_mode)) cs%use\_ebt\_mode = use\_ebt\_mode
1257   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_column\_fraction)) cs%mono\_N2\_column\_fraction = mono\_n2\_column\_fraction
1258   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(mono\_n2\_depth)) cs%mono\_N2\_depth = mono\_n2\_depth
1259   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(remap\_answers\_2018)) cs%remap\_answers\_2018 = remap\_answers\_2018
1260   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(better\_speed\_est)) cs%better\_cg1\_est = better\_speed\_est
1261   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(min\_speed)) cs%min\_speed2 = min\_speed**2
1262   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(wave\_speed\_tol)) cs%wave\_speed\_tol = wave\_speed\_tol
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}}
\subsubsection{\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(\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{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & Thermodynamic variables\\
\hline
\mbox{\tt in}  & {\em nmodes} & Number of modes\\
\hline
\mbox{\tt 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{\tt in}  & {\em full\+\_\+halos} & If true, do the calculation over the entire computational domain.\\
\hline
\mbox{\tt 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{\tt 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{\tt 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}
642   \textcolor{keywordtype}{type}(ocean\_grid\_type),                    \textcolor{keywordtype}{intent(in)}  :: g\textcolor{comment}{ !< Ocean grid structure}
643   \textcolor{keywordtype}{type}(verticalgrid\_type),                  \textcolor{keywordtype}{intent(in)}  :: gv\textcolor{comment}{ !< Vertical grid structure}
644   \textcolor{keywordtype}{type}(unit\_scale\_type),                    \textcolor{keywordtype}{intent(in)}  :: us\textcolor{comment}{ !< A dimensional unit scaling type}
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]}
646   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),                    \textcolor{keywordtype}{intent(in)}  :: tv\textcolor{comment}{ !< Thermodynamic variables}
647   \textcolor{keywordtype}{integer},                                  \textcolor{keywordtype}{intent(in)}  :: nmodes\textcolor{comment}{ !< Number of modes}
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]}
649   \textcolor{keywordtype}{type}(wave\_speed\_cs), \textcolor{keywordtype}{optional},            \textcolor{keywordtype}{pointer}     :: cs\textcolor{comment}{ !< Control structure for MOM\_wave\_speed}
650   \textcolor{keywordtype}{logical},             \textcolor{keywordtype}{optional},            \textcolor{keywordtype}{intent(in)}  :: full\_halos\textcolor{comment}{ !< If true, do the calculation}
651 \textcolor{comment}{                                                                      !! over the entire computational
       domain.}
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}
653 \textcolor{comment}{                                     !! mode speed as the starting point for iterations.}
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}
655 \textcolor{comment}{                                     !! below which 0 is returned [L T-1 ~> m s-1].}
656   \textcolor{keywordtype}{real},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: wave\_speed\_tol\textcolor{comment}{ !< The fractional tolerance for finding the}
657 \textcolor{comment}{                                     !! wave speeds [nondim]}
658 
659   \textcolor{comment}{! Local variables}
660   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G)+1)} :: &
661     drho\_dt, &    \textcolor{comment}{! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]}
662     drho\_ds, &    \textcolor{comment}{! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]}
663     pres, &       \textcolor{comment}{! Interface pressure [R L2 T-2 ~> Pa]}
664     t\_int, &      \textcolor{comment}{! Temperature interpolated to interfaces [degC]}
665     s\_int, &      \textcolor{comment}{! Salinity interpolated to interfaces [ppt]}
666     h\_top, &      \textcolor{comment}{! The distance of each filtered interface from the ocean surface [Z ~> m]}
667     h\_bot, &      \textcolor{comment}{! The distance of each filtered interface from the bottom [Z ~> m]}
668     gprime        \textcolor{comment}{! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].}
669   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G),SZI\_(G))} :: &
670     hf, &         \textcolor{comment}{! Layer thicknesses after very thin layers are combined [Z ~> m]}
671     tf, &         \textcolor{comment}{! Layer temperatures after very thin layers are combined [degC]}
672     sf, &         \textcolor{comment}{! Layer salinities after very thin layers are combined [ppt]}
673     rf            \textcolor{comment}{! Layer densities after very thin layers are combined [R ~> kg m-3]}
674   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZK\_(G))} :: &
675     igl, igu, &   \textcolor{comment}{! The inverse of the reduced gravity across an interface times}
676                   \textcolor{comment}{! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].}
677     hc, &         \textcolor{comment}{! A column of layer thicknesses after convective istabilities are removed [Z ~> m]}
678     tc, &         \textcolor{comment}{! A column of layer temperatures after convective istabilities are removed [degC]}
679     sc, &         \textcolor{comment}{! A column of layer salinites after convective istabilities are removed [ppt]}
680     rc            \textcolor{comment}{! A column of layer densities after convective istabilities are removed [R ~> kg m-3]}
681   \textcolor{keywordtype}{real} :: i\_htot  \textcolor{comment}{! The inverse of the total filtered thicknesses [Z ~> m]}
682   \textcolor{keywordtype}{real} :: c1\_thresh  \textcolor{comment}{! if c1 is below this value, don't bother calculating}
683                      \textcolor{comment}{! cn values for higher modes [L T-1 ~> m s-1]}
684   \textcolor{keywordtype}{real} :: c2\_scale \textcolor{comment}{! A scaling factor for wave speeds to help control the growth of the determinant}
685                    \textcolor{comment}{! and its derivative with lam between rows of the Thomas algorithm solver.  The}
686                    \textcolor{comment}{! exact value should not matter for the final result if it is an even power of 2.}
687   \textcolor{keywordtype}{real} :: det, ddet       \textcolor{comment}{! determinant & its derivative of eigen system}
688   \textcolor{keywordtype}{real} :: lam\_1           \textcolor{comment}{! approximate mode-1 eigenvalue [T2 L-2 ~> s2 m-2]}
689   \textcolor{keywordtype}{real} :: lam\_n           \textcolor{comment}{! approximate mode-n eigenvalue [T2 L-2 ~> s2 m-2]}
690   \textcolor{keywordtype}{real} :: dlam            \textcolor{comment}{! The change in estimates of the eigenvalue [T2 L-2 ~> s m-1]}
691   \textcolor{keywordtype}{real} :: lammin          \textcolor{comment}{! minimum lam value for root searching range [T2 L-2 ~> s2 m-2]}
692   \textcolor{keywordtype}{real} :: lammax          \textcolor{comment}{! maximum lam value for root searching range [T2 L-2 ~> s2 m-2]}
693   \textcolor{keywordtype}{real} :: laminc          \textcolor{comment}{! width of moving window for root searching [T2 L-2 ~> s2 m-2]}
694   \textcolor{keywordtype}{real} :: det\_l,det\_r     \textcolor{comment}{! determinant value at left and right of window}
695   \textcolor{keywordtype}{real} :: ddet\_l,ddet\_r   \textcolor{comment}{! derivative of determinant at left and right of window}
696   \textcolor{keywordtype}{real} :: det\_sub,ddet\_sub\textcolor{comment}{! derivative of determinant at subinterval endpoint}
697   \textcolor{keywordtype}{real} :: xl,xr           \textcolor{comment}{! lam guesses at left and right of window [T2 L-2 ~> s2 m-2]}
698   \textcolor{keywordtype}{real} :: xl\_sub          \textcolor{comment}{! lam guess at left of subinterval window [T2 L-2 ~> s2 m-2]}
699   \textcolor{keywordtype}{real},\textcolor{keywordtype}{dimension(nmodes)} :: &
700           xbl,xbr         \textcolor{comment}{! lam guesses bracketing a zero-crossing (root) [T2 L-2 ~> s2 m-2]}
701   \textcolor{keywordtype}{integer} :: numint       \textcolor{comment}{! number of widows (intervals) in root searching range}
702   \textcolor{keywordtype}{integer} :: nrootsfound  \textcolor{comment}{! number of extra roots found (not including 1st root)}
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]}
704   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
705     htot, hmin, &  \textcolor{comment}{! Thicknesses [Z ~> m]}
706     h\_here, &      \textcolor{comment}{! A thickness [Z ~> m]}
707     hxt\_here, &    \textcolor{comment}{! A layer integrated temperature [degC Z ~> degC m]}
708     hxs\_here, &    \textcolor{comment}{! A layer integrated salinity [ppt Z ~> ppt m]}
709     hxr\_here       \textcolor{comment}{! A layer integrated density [R Z ~> kg m-2]}
710   \textcolor{keywordtype}{real} :: speed2\_tot \textcolor{comment}{! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]}
711   \textcolor{keywordtype}{real} :: speed2\_min \textcolor{comment}{! minimum mode speed (squared) to consider in root searching [L2 T-2 ~> m2 s-2]}
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]}
713   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: reduct\_factor = 0.5
714                      \textcolor{comment}{! A factor used in setting speed2\_min [nondim]}
715   \textcolor{keywordtype}{real} :: i\_hnew   \textcolor{comment}{! The inverse of a new layer thickness [Z-1 ~> m-1]}
716   \textcolor{keywordtype}{real} :: drxh\_sum \textcolor{comment}{! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2]}
717   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{pointer}, \textcolor{keywordtype}{dimension(:,:,:)} :: t => null(), s => null()
718   \textcolor{keywordtype}{real} :: g\_rho0   \textcolor{comment}{! G\_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].}
719   \textcolor{keywordtype}{real} :: tol\_hfrac  \textcolor{comment}{! Layers that together are smaller than this fraction of}
720                      \textcolor{comment}{! the total water column can be merged for efficiency.}
721   \textcolor{keywordtype}{real} :: min\_h\_frac \textcolor{comment}{! tol\_Hfrac divided by the total number of layers [nondim].}
722   \textcolor{keywordtype}{real} :: tol\_solve  \textcolor{comment}{! The fractional tolerance with which to solve for the wave speeds [nondim].}
723   \textcolor{keywordtype}{real} :: tol\_merge  \textcolor{comment}{! The fractional change in estimated wave speed that is allowed}
724                      \textcolor{comment}{! when deciding to merge layers in the calculation [nondim]}
725   \textcolor{keywordtype}{integer} :: kf(szi\_(g)) \textcolor{comment}{! The number of active layers after filtering.}
726   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: max\_itt = 10
727   \textcolor{keywordtype}{logical} :: use\_eos    \textcolor{comment}{! If true, density is calculated from T & S using the equation of state.}
728   \textcolor{keywordtype}{logical} :: better\_est \textcolor{comment}{! If true, use an improved estimate of the first mode internal wave speed.}
729   \textcolor{keywordtype}{logical} :: merge      \textcolor{comment}{! If true, merge the current layer with the one above.}
730   \textcolor{keywordtype}{integer} :: nsub       \textcolor{comment}{! number of subintervals used for root finding}
731   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{parameter} :: sub\_it\_max = 4
732                         \textcolor{comment}{! maximum number of times to subdivide interval}
733                         \textcolor{comment}{! for root finding (# intervals = 2**sub\_it\_max)}
734   \textcolor{keywordtype}{logical} :: sub\_rootfound \textcolor{comment}{! if true, subdivision has located root}
735   \textcolor{keywordtype}{integer} :: kc         \textcolor{comment}{! The number of layers in the column after merging}
736   \textcolor{keywordtype}{integer} :: sub, sub\_it
737   \textcolor{keywordtype}{integer} :: i, j, k, k2, itt, is, ie, js, je, nz, row, iint, m, ig, jg
738 
739   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
740 
741   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) \textcolor{keywordflow}{then}
742     \textcolor{keywordflow}{if} (.not. \textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_wave\_speed: "}// &
743            \textcolor{stringliteral}{"Module must be initialized before it is used."})
744 \textcolor{keywordflow}{  endif}
745 
746   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(full\_halos)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (full\_halos) \textcolor{keywordflow}{then}
747     is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
748 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
749 
750   s => tv%S ; t => tv%T
751   g\_rho0 = gv%g\_Earth / gv%Rho0
752   \textcolor{comment}{! Simplifying the following could change answers at roundoff.}
753   z\_to\_pres = gv%Z\_to\_H * (gv%H\_to\_RZ * gv%g\_Earth)
754   use\_eos = \textcolor{keyword}{associated}(tv%eqn\_of\_state)
755   c1\_thresh = 0.01*us%m\_s\_to\_L\_T
756   c2\_scale = us%m\_s\_to\_L\_T**2 / 4096.0**2 \textcolor{comment}{! Other powers of 2 give identical results.}
757 
758   better\_est = .false. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) better\_est = cs%better\_cg1\_est
759   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(better\_speed\_est)) better\_est = better\_speed\_est
760   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
761     tol\_solve = 0.001 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) tol\_solve = cs%wave\_speed\_tol
762     \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(wave\_speed\_tol)) tol\_solve = wave\_speed\_tol
763     tol\_hfrac  = 0.1*tol\_solve ; tol\_merge = tol\_solve / \textcolor{keywordtype}{real}(nz)
764   \textcolor{keywordflow}{else}
765     tol\_solve = 0.001 ; tol\_hfrac  = 0.0001 ; tol\_merge = 0.001
766 \textcolor{keywordflow}{  endif}
767   cg1\_min2 = 0.0 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(cs)) cg1\_min2 = cs%min\_speed2
768   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(min\_speed)) cg1\_min2 = min\_speed**2
769 
770   \textcolor{comment}{! Zero out all wave speeds.  Values over land or for columns that are too weakly stratified}
771   \textcolor{comment}{! are not changed from this zero value.}
772   cn(:,:,:) = 0.0
773 
774   min\_h\_frac = tol\_hfrac / \textcolor{keywordtype}{real}(nz)
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, &}
776   \textcolor{comment}{!$OMP                                     Z\_to\_pres,tv,cn,g\_Rho0,nmodes,cg1\_min2,better\_est, &}
777   \textcolor{comment}{!$OMP                                     c1\_thresh,tol\_solve,tol\_merge)}
778   \textcolor{keywordflow}{do} j=js,je
779     \textcolor{comment}{!   First merge very thin layers with the one above (or below if they are}
780     \textcolor{comment}{! at the top).  This also transposes the row order so that columns can}
781     \textcolor{comment}{! be worked upon one at a time.}
782     \textcolor{keywordflow}{do} i=is,ie ; htot(i) = 0.0 ;\textcolor{keywordflow}{ enddo}
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}
784 
785     \textcolor{keywordflow}{do} i=is,ie
786       hmin(i) = htot(i)*min\_h\_frac ; kf(i) = 1 ; h\_here(i) = 0.0
787       hxt\_here(i) = 0.0 ; hxs\_here(i) = 0.0 ; hxr\_here(i) = 0.0
788 \textcolor{keywordflow}{    enddo}
789     \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
790       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
791         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}
792           hf(kf(i),i) = h\_here(i)
793           tf(kf(i),i) = hxt\_here(i) / h\_here(i)
794           sf(kf(i),i) = hxs\_here(i) / h\_here(i)
795           kf(i) = kf(i) + 1
796 
797           \textcolor{comment}{! Start a new layer}
798           h\_here(i) = h(i,j,k)*gv%H\_to\_Z
799           hxt\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*t(i,j,k)
800           hxs\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*s(i,j,k)
801         \textcolor{keywordflow}{else}
802           h\_here(i) = h\_here(i) + h(i,j,k)*gv%H\_to\_Z
803           hxt\_here(i) = hxt\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*t(i,j,k)
804           hxs\_here(i) = hxs\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*s(i,j,k)
805 \textcolor{keywordflow}{        endif}
806 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
807       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}
808         hf(kf(i),i) = h\_here(i)
809         tf(kf(i),i) = hxt\_here(i) / h\_here(i)
810         sf(kf(i),i) = hxs\_here(i) / h\_here(i)
811 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
812     \textcolor{keywordflow}{else}
813       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
814         \textcolor{keywordflow}{if} ((h\_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H\_to\_Z > hmin(i))) \textcolor{keywordflow}{then}
815           hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)
816           kf(i) = kf(i) + 1
817 
818           \textcolor{comment}{! Start a new layer}
819           h\_here(i) = h(i,j,k)*gv%H\_to\_Z
820           hxr\_here(i) = (h(i,j,k)*gv%H\_to\_Z)*gv%Rlay(k)
821         \textcolor{keywordflow}{else}
822           h\_here(i) = h\_here(i) + h(i,j,k)*gv%H\_to\_Z
823           hxr\_here(i) = hxr\_here(i) + (h(i,j,k)*gv%H\_to\_Z)*gv%Rlay(k)
824 \textcolor{keywordflow}{        endif}
825 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
826       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (h\_here(i) > 0.0) \textcolor{keywordflow}{then}
827         hf(kf(i),i) = h\_here(i) ; rf(kf(i),i) = hxr\_here(i) / h\_here(i)
828 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
829 \textcolor{keywordflow}{    endif}
830 
831     \textcolor{comment}{! From this point, we can work on individual columns without causing memory to have page faults.}
832     \textcolor{keywordflow}{do} i=is,ie
833       \textcolor{keywordflow}{if} (g%mask2dT(i,j) > 0.5) \textcolor{keywordflow}{then}
834         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
835           pres(1) = 0.0 ; h\_top(1) = 0.0
836           \textcolor{keywordflow}{do} k=2,kf(i)
837             pres(k) = pres(k-1) + z\_to\_pres*hf(k-1,i)
838             t\_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
839             s\_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
840             h\_top(k) = h\_top(k-1) + hf(k-1,i)
841 \textcolor{keywordflow}{          enddo}
842           \textcolor{keyword}{call }calculate\_density\_derivs(t\_int, s\_int, pres, drho\_dt, drho\_ds, &
843                                         tv%eqn\_of\_state, (/2,kf(i)/) )
844 
845           \textcolor{comment}{! Sum the reduced gravities to find out how small a density difference is negligibly small.}
846           drxh\_sum = 0.0
847           \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
848             \textcolor{comment}{! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for}
849             \textcolor{comment}{! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.}
850             \textcolor{comment}{! For a uniform stratification and a huge number of layers uniformly distributed in}
851             \textcolor{comment}{! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.}
852             \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}
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.}
854               h\_bot(kf(i)+1) = 0.0
855               \textcolor{keywordflow}{do} k=kf(i),2,-1
856                 h\_bot(k) = h\_bot(k+1) + hf(k,i)
857                 drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * &
858                     max(0.0, drho\_dt(k)*(tf(k,i)-tf(k-1,i)) + drho\_ds(k)*(sf(k,i)-sf(k-1,i)))
859 \textcolor{keywordflow}{              enddo}
860 \textcolor{keywordflow}{            endif}
861           \textcolor{keywordflow}{else}
862             \textcolor{comment}{! This estimate is problematic in that it goes like 1/nz for a large number of layers,}
863             \textcolor{comment}{! but it is an overestimate (as desired) for a small number of layers, by at a factor}
864             \textcolor{comment}{! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.}
865             \textcolor{keywordflow}{do} k=2,kf(i)
866               drxh\_sum = drxh\_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
867                   max(0.0, drho\_dt(k)*(tf(k,i)-tf(k-1,i)) + drho\_ds(k)*(sf(k,i)-sf(k-1,i)))
868 \textcolor{keywordflow}{            enddo}
869 \textcolor{keywordflow}{          endif}
870         \textcolor{keywordflow}{else}
871           drxh\_sum = 0.0
872           \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
873             h\_top(1) = 0.0
874             \textcolor{keywordflow}{do} k=2,kf(i) ; h\_top(k) = h\_top(k-1) + hf(k-1,i) ;\textcolor{keywordflow}{ enddo}
875               \textcolor{keywordflow}{if} (h\_top(kf(i)) > 0.0) \textcolor{keywordflow}{then}
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.}
877               h\_bot(kf(i)+1) = 0.0
878               \textcolor{keywordflow}{do} k=kf(i),2,-1
879                 h\_bot(k) = h\_bot(k+1) + hf(k,i)
880                 drxh\_sum = drxh\_sum + ((h\_top(k) * h\_bot(k)) * i\_htot) * max(0.0,rf(k,i)-rf(k-1,i))
881 \textcolor{keywordflow}{              enddo}
882 \textcolor{keywordflow}{            endif}
883           \textcolor{keywordflow}{else}
884             \textcolor{keywordflow}{do} k=2,kf(i)
885               drxh\_sum = drxh\_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
886 \textcolor{keywordflow}{            enddo}
887 \textcolor{keywordflow}{          endif}
888 \textcolor{keywordflow}{        endif}
889 
890         \textcolor{comment}{!   Find gprime across each internal interface, taking care of convective}
891         \textcolor{comment}{! instabilities by merging layers.}
892         \textcolor{keywordflow}{if} (g\_rho0 * drxh\_sum > cg1\_min2) \textcolor{keywordflow}{then}
893           \textcolor{comment}{! Merge layers to eliminate convective instabilities or exceedingly}
894           \textcolor{comment}{! small reduced gravities.  Merging layers reduces the estimated wave speed by}
895           \textcolor{comment}{! (rho(2)-rho(1))*h(1)*h(2) / H\_tot.}
896           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
897             kc = 1
898             hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
899             \textcolor{keywordflow}{do} k=2,kf(i)
900               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
901                 merge = ((drho\_dt(k)*(tf(k,i)-tc(kc)) + drho\_ds(k)*(sf(k,i)-sc(kc))) * &
902                          ((hc(kc) * hf(k,i))*i\_htot) < 2.0 * tol\_merge*drxh\_sum)
903               \textcolor{keywordflow}{else}
904                 merge = ((drho\_dt(k)*(tf(k,i)-tc(kc)) + drho\_ds(k)*(sf(k,i)-sc(kc))) * &
905                          (hc(kc) + hf(k,i)) < 2.0 * tol\_merge*drxh\_sum)
906 \textcolor{keywordflow}{              endif}
907               \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
908                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}
909                 i\_hnew = 1.0 / (hc(kc) + hf(k,i))
910                 tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i\_hnew
911                 sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i\_hnew
912                 hc(kc) = (hc(kc) + hf(k,i))
913                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}
914                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}
915                 \textcolor{comment}{! far back we go.}
916                 \textcolor{keywordflow}{do} k2=kc,2,-1
917                   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
918                     merge = ((drho\_dt(k2)*(tc(k2)-tc(k2-1)) + drho\_ds(k2)*(sc(k2)-sc(k2-1))) * &
919                              ((hc(k2) * hc(k2-1))*i\_htot) < tol\_merge*drxh\_sum)
920                   \textcolor{keywordflow}{else}
921                     merge = ((drho\_dt(k2)*(tc(k2)-tc(k2-1)) + drho\_ds(k2)*(sc(k2)-sc(k2-1))) * &
922                              (hc(k2) + hc(k2-1)) < tol\_merge*drxh\_sum)
923 \textcolor{keywordflow}{                  endif}
924                   \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
925                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}
926                     i\_hnew = 1.0 / (hc(kc) + hc(kc-1))
927                     tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i\_hnew
928                     sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i\_hnew
929                     hc(kc-1) = (hc(kc) + hc(kc-1))
930                     kc = kc - 1
931                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}
932 \textcolor{keywordflow}{                enddo}
933               \textcolor{keywordflow}{else}
934                 \textcolor{comment}{! Add a new layer to the column.}
935                 kc = kc + 1
936                 drho\_ds(kc) = drho\_ds(k) ; drho\_dt(kc) = drho\_dt(k)
937                 tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
938 \textcolor{keywordflow}{              endif}
939 \textcolor{keywordflow}{            enddo}
940             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}
941             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-Boussinesq.}
942               gprime(k) = g\_rho0 * (drho\_dt(k)*(tc(k)-tc(k-1)) + drho\_ds(k)*(sc(k)-sc(k-1)))
943 \textcolor{keywordflow}{            enddo}
944           \textcolor{keywordflow}{else}  \textcolor{comment}{! .not.use\_EOS}
945             \textcolor{comment}{! Do the same with density directly...}
946             kc = 1
947             hc(1) = hf(1,i) ; rc(1) = rf(1,i)
948             \textcolor{keywordflow}{do} k=2,kf(i)
949               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
950                 merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i\_htot) < 2.0 * tol\_merge*drxh\_sum)
951               \textcolor{keywordflow}{else}
952                 merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol\_merge*drxh\_sum)
953 \textcolor{keywordflow}{              endif}
954               \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
955                 \textcolor{comment}{! Merge this layer with the one above and backtrack.}
956                 rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
957                 hc(kc) = (hc(kc) + hf(k,i))
958                 \textcolor{comment}{! Backtrack to remove any convective instabilities above...  Note}
959                 \textcolor{comment}{! that the tolerance is a factor of two larger, to avoid limit how}
960                 \textcolor{comment}{! far back we go.}
961                 \textcolor{keywordflow}{do} k2=kc,2,-1
962                   \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
963                     merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i\_htot) < tol\_merge*drxh\_sum)
964                   \textcolor{keywordflow}{else}
965                     merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol\_merge*drxh\_sum)
966 \textcolor{keywordflow}{                  endif}
967                   \textcolor{keywordflow}{if} (merge) \textcolor{keywordflow}{then}
968                     \textcolor{comment}{! Merge the two bottommost layers.  At this point kc = k2.}
969                     rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
970                     hc(kc-1) = (hc(kc) + hc(kc-1))
971                     kc = kc - 1
972                   \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ endif}
973 \textcolor{keywordflow}{                enddo}
974               \textcolor{keywordflow}{else}
975                 \textcolor{comment}{! Add a new layer to the column.}
976                 kc = kc + 1
977                 rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
978 \textcolor{keywordflow}{              endif}
979 \textcolor{keywordflow}{            enddo}
980             \textcolor{comment}{! At this point there are kc layers and the gprimes should be positive.}
981             \textcolor{keywordflow}{do} k=2,kc \textcolor{comment}{! Revisit this if non-Boussinesq.}
982               gprime(k) = g\_rho0 * (rc(k)-rc(k-1))
983 \textcolor{keywordflow}{            enddo}
984 \textcolor{keywordflow}{          endif}  \textcolor{comment}{! use\_EOS}
985 
986           \textcolor{comment}{!-----------------NOW FIND WAVE SPEEDS---------------------------------------}
987           \textcolor{comment}{! ig = i + G%idg\_offset ; jg = j + G%jdg\_offset}
988           \textcolor{comment}{!   Sum the contributions from all of the interfaces to give an over-estimate}
989           \textcolor{comment}{! of the first-mode wave speed.  Also populate Igl and Igu which are the}
990           \textcolor{comment}{! non-leading diagonals of the tridiagonal matrix.}
991           \textcolor{keywordflow}{if} (kc >= 2) \textcolor{keywordflow}{then}
992             \textcolor{comment}{! initialize speed2\_tot}
993             speed2\_tot = 0.0
994             \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
995               h\_top(1) = 0.0 ; h\_bot(kc+1) = 0.0
996               \textcolor{keywordflow}{do} k=2,kc+1 ; h\_top(k) = h\_top(k-1) + hc(k-1) ;\textcolor{keywordflow}{ enddo}
997               \textcolor{keywordflow}{do} k=kc,2,-1 ; h\_bot(k) = h\_bot(k+1) + hc(k) ;\textcolor{keywordflow}{ enddo}
998               i\_htot = 0.0 ; \textcolor{keywordflow}{if} (h\_top(kc+1) > 0.0) i\_htot = 1.0 / h\_top(kc+1)
999 \textcolor{keywordflow}{            endif}
1000 
1001             \textcolor{comment}{! Calculate Igu, Igl, depth, and N2 at each interior interface}
1002             \textcolor{comment}{! [excludes surface (K=1) and bottom (K=kc+1)]}
1003             \textcolor{keywordflow}{do} k=2,kc
1004               igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
1005               \textcolor{keywordflow}{if} (better\_est) \textcolor{keywordflow}{then}
1006                 speed2\_tot = speed2\_tot + gprime(k)*((h\_top(k) * h\_bot(k)) * i\_htot)
1007               \textcolor{keywordflow}{else}
1008                 speed2\_tot = speed2\_tot + gprime(k)*(hc(k-1)+hc(k))
1009 \textcolor{keywordflow}{              endif}
1010 \textcolor{keywordflow}{            enddo}
1011 
1012             \textcolor{comment}{! Under estimate the first eigenvalue (overestimate the speed) to start with.}
1013             lam\_1 = 1.0 / speed2\_tot
1014 
1015             \textcolor{comment}{! Find the first eigen value}
1016             \textcolor{keywordflow}{do} itt=1,max\_itt
1017               \textcolor{comment}{! calculate the determinant of (A-lam\_1*I)}
1018               \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lam\_1, det, ddet, row\_scale=c2\_scale)
1019 
1020               \textcolor{comment}{! If possible, use Newton's method iteration to find a new estimate of lam\_1}
1021               \textcolor{comment}{!det = det\_it(itt) ; ddet = ddet\_it(itt)}
1022               \textcolor{keywordflow}{if} ((ddet >= 0.0) .or. (-det > -0.5*lam\_1*ddet)) \textcolor{keywordflow}{then}
1023                 \textcolor{comment}{! lam\_1 was not an under-estimate, as intended, so Newton's method}
1024                 \textcolor{comment}{! may not be reliable; lam\_1 must be reduced, but not by more than half.}
1025                 lam\_1 = 0.5 * lam\_1
1026                 dlam = -lam\_1
1027               \textcolor{keywordflow}{else}  \textcolor{comment}{! Newton's method is OK.}
1028                 dlam = - det / ddet
1029                 lam\_1 = lam\_1 + dlam
1030 \textcolor{keywordflow}{              endif}
1031 
1032               \textcolor{keywordflow}{if} (abs(dlam) < tol\_solve*lam\_1) \textcolor{keywordflow}{exit}
1033 \textcolor{keywordflow}{            enddo}
1034 
1035             \textcolor{keywordflow}{if} (lam\_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam\_1)
1036 
1037             \textcolor{comment}{! Find other eigen values if c1 is of significant magnitude, > cn\_thresh}
1038             nrootsfound = 0    \textcolor{comment}{! number of extra roots found (not including 1st root)}
1039             \textcolor{keywordflow}{if} (nmodes>1 .and. kc>=nmodes+1 .and. cn(i,j,1)>c1\_thresh) \textcolor{keywordflow}{then}
1040               \textcolor{comment}{! Set the the range to look for the other desired eigen values}
1041               \textcolor{comment}{! set min value just greater than the 1st root (found above)}
1042               lammin = lam\_1*(1.0 + tol\_solve)
1043               \textcolor{comment}{! set max value based on a low guess at wavespeed for highest mode}
1044               speed2\_min = (reduct\_factor*cn(i,j,1)/\textcolor{keywordtype}{real}(nmodes))**2
1045               lammax = 1.0 / speed2\_min
1046               \textcolor{comment}{! set width of interval (not sure about this - BDM)}
1047               laminc = 0.5*lam\_1
1048               \textcolor{comment}{! set number of intervals within search range}
1049               numint = nint((lammax - lammin)/laminc)
1050 
1051               \textcolor{comment}{!   Find intervals containing zero-crossings (roots) of the determinant}
1052               \textcolor{comment}{! that are beyond the first root}
1053 
1054               \textcolor{comment}{! find det\_l of first interval (det at left endpoint)}
1055               \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lammin, det\_l, ddet\_l, row\_scale=c2\_scale)
1056               \textcolor{comment}{! move interval window looking for zero-crossings************************}
1057               \textcolor{keywordflow}{do} iint=1,numint
1058                 xr = lammin + laminc * iint
1059                 xl = xr - laminc
1060                 \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, xr, det\_r, ddet\_r, row\_scale=c2\_scale)
1061                 \textcolor{keywordflow}{if} (det\_l*det\_r < 0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! if function changes sign}
1062                   \textcolor{keywordflow}{if} (det\_l*ddet\_l < 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! if function at left is headed to zero}
1063                     nrootsfound = nrootsfound + 1
1064                     xbl(nrootsfound) = xl
1065                     xbr(nrootsfound) = xr
1066                   \textcolor{keywordflow}{else}
1067                     \textcolor{comment}{!   function changes sign but has a local max/min in interval,}
1068                     \textcolor{comment}{! try subdividing interval as many times as necessary (or sub\_it\_max).}
1069                     \textcolor{comment}{! loop that increases number of subintervals:}
1070                     \textcolor{comment}{!call MOM\_error(WARNING, "determinant changes sign"// &}
1071                     \textcolor{comment}{!            "but has a local max/min in interval;"//&}
1072                     \textcolor{comment}{!            " reduce increment in lam.")}
1073                     \textcolor{comment}{! begin subdivision loop -------------------------------------------}
1074                     sub\_rootfound = .false. \textcolor{comment}{! initialize}
1075                     \textcolor{keywordflow}{do} sub\_it=1,sub\_it\_max
1076                       nsub = 2**sub\_it \textcolor{comment}{! number of subintervals; nsub=2,4,8,...}
1077                       \textcolor{comment}{! loop over each subinterval:}
1078                       \textcolor{keywordflow}{do} sub=1,nsub-1,2 \textcolor{comment}{! only check odds; sub = 1; 1,3; 1,3,5,7;...}
1079                         xl\_sub = xl + laminc/(nsub)*sub
1080                         \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, xl\_sub, det\_sub, ddet\_sub, &
1081                                          row\_scale=c2\_scale)
1082                         \textcolor{keywordflow}{if} (det\_sub*det\_r < 0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! if function changes sign}
1083                           \textcolor{keywordflow}{if} (det\_sub*ddet\_sub < 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! if function at left is headed to zero}
1084                             sub\_rootfound = .true.
1085                             nrootsfound = nrootsfound + 1
1086                             xbl(nrootsfound) = xl\_sub
1087                             xbr(nrootsfound) = xr
1088                             \textcolor{keywordflow}{exit} \textcolor{comment}{! exit sub loop}
1089 \textcolor{keywordflow}{                          endif} \textcolor{comment}{! headed toward zero}
1090 \textcolor{keywordflow}{                        endif} \textcolor{comment}{! sign change}
1091 \textcolor{keywordflow}{                      enddo} \textcolor{comment}{! sub-loop}
1092                       \textcolor{keywordflow}{if} (sub\_rootfound) \textcolor{keywordflow}{exit} \textcolor{comment}{! root has been found, exit sub\_it loop}
1093                       \textcolor{comment}{!   Otherwise, function changes sign but has a local max/min in one of the}
1094                       \textcolor{comment}{! sub intervals, try subdividing again unless sub\_it\_max has been reached.}
1095                       \textcolor{keywordflow}{if} (sub\_it == sub\_it\_max) \textcolor{keywordflow}{then}
1096                         \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"wave\_speed: root not found "}// &
1097                                        \textcolor{stringliteral}{" after sub\_it\_max subdivisions of original"}// &
1098                                        \textcolor{stringliteral}{" interval."})
1099 \textcolor{keywordflow}{                      endif} \textcolor{comment}{! sub\_it == sub\_it\_max}
1100 \textcolor{keywordflow}{                    enddo} \textcolor{comment}{! sub\_it-loop-------------------------------------------------}
1101 \textcolor{keywordflow}{                  endif} \textcolor{comment}{! det\_l*ddet\_l < 0.0}
1102 \textcolor{keywordflow}{                endif} \textcolor{comment}{! det\_l*det\_r < 0.0}
1103                 \textcolor{comment}{! exit iint-loop if all desired roots have been found}
1104                 \textcolor{keywordflow}{if} (nrootsfound >= nmodes-1) \textcolor{keywordflow}{then}
1105                   \textcolor{comment}{! exit if all additional roots found}
1106                   \textcolor{keywordflow}{exit}
1107                 \textcolor{keywordflow}{elseif} (iint == numint) \textcolor{keywordflow}{then}
1108                   \textcolor{comment}{! oops, lamMax not large enough - could add code to increase (BDM)}
1109                   \textcolor{comment}{! set unfound modes to zero for now (BDM)}
1110                   \textcolor{comment}{!   cn(i,j,nrootsfound+2:nmodes) = 0.0}
1111                 \textcolor{keywordflow}{else}
1112                   \textcolor{comment}{! else shift interval and keep looking until nmodes or numint is reached}
1113                   det\_l = det\_r
1114                   ddet\_l = ddet\_r
1115 \textcolor{keywordflow}{                endif}
1116 \textcolor{keywordflow}{              enddo} \textcolor{comment}{! iint-loop}
1117 
1118               \textcolor{comment}{! Use Newton's method to find the roots within the identified windows}
1119               \textcolor{keywordflow}{do} m=1,nrootsfound \textcolor{comment}{! loop over the root-containing widows (excluding 1st mode)}
1120                 lam\_n = xbl(m) \textcolor{comment}{! first guess is left edge of window}
1121                 \textcolor{keywordflow}{do} itt=1,max\_itt
1122                   \textcolor{comment}{! calculate the determinant of (A-lam\_n*I)}
1123                   \textcolor{keyword}{call }tridiag\_det(igu, igl, 2, kc, lam\_n, det, ddet, row\_scale=c2\_scale)
1124                   \textcolor{comment}{! Use Newton's method to find a new estimate of lam\_n}
1125                   dlam = - det / ddet
1126                   lam\_n = lam\_n + dlam
1127                   \textcolor{keywordflow}{if} (abs(dlam) < tol\_solve*lam\_1)  \textcolor{keywordflow}{exit}
1128 \textcolor{keywordflow}{                enddo} \textcolor{comment}{! itt-loop}
1129                 \textcolor{comment}{! calculate nth mode speed}
1130                 \textcolor{keywordflow}{if} (lam\_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam\_n)
1131 \textcolor{keywordflow}{              enddo} \textcolor{comment}{! n-loop}
1132 \textcolor{keywordflow}{            endif} \textcolor{comment}{! if nmodes>1 .and. kc>nmodes .and. c1>c1\_thresh}
1133 \textcolor{keywordflow}{          endif} \textcolor{comment}{! if more than 2 layers}
1134 \textcolor{keywordflow}{        endif} \textcolor{comment}{! if drxh\_sum < 0}
1135 \textcolor{keywordflow}{      endif} \textcolor{comment}{! if not land}
1136 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! i-loop}
1137 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-loop}
1138 
\end{DoxyCode}
