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