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