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