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