\hypertarget{namespacemom__density__integrals}{}\section{mom\+\_\+density\+\_\+integrals Module Reference} \label{namespacemom__density__integrals}\index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsection{Detailed Description} Provides integrals of density. \subsection*{Functions/\+Subroutines} \begin{DoxyCompactItemize} \item subroutine, public \hyperlink{namespacemom__density__integrals_ac36ae5f4af2d02df0a1adf41b762e017}{int\+\_\+density\+\_\+dz} (T, S, z\+\_\+t, z\+\_\+b, rho\+\_\+ref, rho\+\_\+0, G\+\_\+e, HI, E\+OS, US, dpa, intz\+\_\+dpa, intx\+\_\+dpa, inty\+\_\+dpa, bathyT, dz\+\_\+neglect, use\+Mass\+Wght\+Interp) \begin{DoxyCompactList}\small\item\em Calls the appropriate subroutine to calculate analytical and nearly-\/analytical integrals in z across layers of pressure anomalies, which are required for calculating the finite-\/volume form pressure accelerations in a Boussinesq model. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__density__integrals_ad0ff54518bfacf7beebbe4dba687a914}{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+pcm} (T, S, z\+\_\+t, z\+\_\+b, rho\+\_\+ref, rho\+\_\+0, G\+\_\+e, HI, E\+OS, US, dpa, intz\+\_\+dpa, intx\+\_\+dpa, inty\+\_\+dpa, bathyT, dz\+\_\+neglect, use\+Mass\+Wght\+Interp) \begin{DoxyCompactList}\small\item\em Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which are required for calculating the finite-\/volume form pressure accelerations in a Boussinesq model. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__density__integrals_a454d4d62eb599716cc3389fb8aa90b4b}{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+plm} (k, tv, T\+\_\+t, T\+\_\+b, S\+\_\+t, S\+\_\+b, e, rho\+\_\+ref, rho\+\_\+0, G\+\_\+e, dz\+\_\+subroundoff, bathyT, HI, GV, E\+OS, US, dpa, intz\+\_\+dpa, intx\+\_\+dpa, inty\+\_\+dpa, use\+Mass\+Wght\+Interp) \begin{DoxyCompactList}\small\item\em Compute pressure gradient force integrals by quadrature for the case where T and S are linear profiles. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__density__integrals_a3bf090e5cbb58811b3dd0ee4de80f999}{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+ppm} (k, tv, T\+\_\+t, T\+\_\+b, S\+\_\+t, S\+\_\+b, e, rho\+\_\+ref, rho\+\_\+0, G\+\_\+e, dz\+\_\+subroundoff, bathyT, HI, GV, E\+OS, US, dpa, intz\+\_\+dpa, intx\+\_\+dpa, inty\+\_\+dpa, use\+Mass\+Wght\+Interp) \begin{DoxyCompactList}\small\item\em Compute pressure gradient force integrals for layer \char`\"{}k\char`\"{} and the case where T and S are parabolic profiles. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__density__integrals_a759c2ae7aec17c59d532050f68a1e518}{int\+\_\+specific\+\_\+vol\+\_\+dp} (T, S, p\+\_\+t, p\+\_\+b, alpha\+\_\+ref, HI, E\+OS, US, dza, intp\+\_\+dza, intx\+\_\+dza, inty\+\_\+dza, halo\+\_\+size, bathyP, d\+P\+\_\+tiny, use\+Mass\+Wght\+Interp) \begin{DoxyCompactList}\small\item\em Calls the appropriate subroutine to calculate analytical and nearly-\/analytical integrals in pressure across layers of geopotential anomalies, which are required for calculating the finite-\/volume form pressure accelerations in a non-\/\+Boussinesq model. There are essentially no free assumptions, apart from the use of Boole\textquotesingle{}s rule to do the horizontal integrals, and from a truncation in the series for log(1-\/eps/1+eps) that assumes that $\vert$eps$\vert$ $<$ 0.\+34. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__density__integrals_a86ebdfebaaea2ac0339dcabed482dd11}{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+pcm} (T, S, p\+\_\+t, p\+\_\+b, alpha\+\_\+ref, HI, E\+OS, US, dza, intp\+\_\+dza, intx\+\_\+dza, inty\+\_\+dza, halo\+\_\+size, bathyP, d\+P\+\_\+neglect, use\+Mass\+Wght\+Interp) \begin{DoxyCompactList}\small\item\em This subroutine calculates integrals of specific volume anomalies in pressure across layers, which are required for calculating the finite-\/volume form pressure accelerations in a non-\/\+Boussinesq model. There are essentially no free assumptions, apart from the use of Boole\textquotesingle{}s rule quadrature to do the integrals. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__density__integrals_a71a518f80b5f2ecf83fd9c638cad140d}{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+plm} (T\+\_\+t, T\+\_\+b, S\+\_\+t, S\+\_\+b, p\+\_\+t, p\+\_\+b, alpha\+\_\+ref, d\+P\+\_\+neglect, bathyP, HI, E\+OS, US, dza, intp\+\_\+dza, intx\+\_\+dza, inty\+\_\+dza, use\+Mass\+Wght\+Interp) \begin{DoxyCompactList}\small\item\em This subroutine calculates integrals of specific volume anomalies in pressure across layers, which are required for calculating the finite-\/volume form pressure accelerations in a non-\/\+Boussinesq model. There are essentially no free assumptions, apart from the use of Boole\textquotesingle{}s rule quadrature to do the integrals. \end{DoxyCompactList}\item subroutine, public \hyperlink{namespacemom__density__integrals_af9946b0e15d53d8f1fd1ab1c1f81f3a2}{find\+\_\+depth\+\_\+of\+\_\+pressure\+\_\+in\+\_\+cell} (T\+\_\+t, T\+\_\+b, S\+\_\+t, S\+\_\+b, z\+\_\+t, z\+\_\+b, P\+\_\+t, P\+\_\+tgt, rho\+\_\+ref, G\+\_\+e, E\+OS, US, P\+\_\+b, z\+\_\+out, z\+\_\+tol) \begin{DoxyCompactList}\small\item\em Find the depth at which the reconstructed pressure matches P\+\_\+tgt. \end{DoxyCompactList}\item real function \hyperlink{namespacemom__density__integrals_abb8d1f2c24def176a27ef2d30fa373df}{frac\+\_\+dp\+\_\+at\+\_\+pos} (T\+\_\+t, T\+\_\+b, S\+\_\+t, S\+\_\+b, z\+\_\+t, z\+\_\+b, rho\+\_\+ref, G\+\_\+e, pos, E\+OS) \begin{DoxyCompactList}\small\item\em Returns change in anomalous pressure change from top to non-\/dimensional position pos between z\+\_\+t and z\+\_\+b. \end{DoxyCompactList}\end{DoxyCompactItemize} \subsection{Function/\+Subroutine Documentation} \mbox{\Hypertarget{namespacemom__density__integrals_af9946b0e15d53d8f1fd1ab1c1f81f3a2}\label{namespacemom__density__integrals_af9946b0e15d53d8f1fd1ab1c1f81f3a2}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!find\+\_\+depth\+\_\+of\+\_\+pressure\+\_\+in\+\_\+cell@{find\+\_\+depth\+\_\+of\+\_\+pressure\+\_\+in\+\_\+cell}} \index{find\+\_\+depth\+\_\+of\+\_\+pressure\+\_\+in\+\_\+cell@{find\+\_\+depth\+\_\+of\+\_\+pressure\+\_\+in\+\_\+cell}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{find\+\_\+depth\+\_\+of\+\_\+pressure\+\_\+in\+\_\+cell()}{find\_depth\_of\_pressure\_in\_cell()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::find\+\_\+depth\+\_\+of\+\_\+pressure\+\_\+in\+\_\+cell (\begin{DoxyParamCaption}\item[{real, intent(in)}]{T\+\_\+t, }\item[{real, intent(in)}]{T\+\_\+b, }\item[{real, intent(in)}]{S\+\_\+t, }\item[{real, intent(in)}]{S\+\_\+b, }\item[{real, intent(in)}]{z\+\_\+t, }\item[{real, intent(in)}]{z\+\_\+b, }\item[{real, intent(in)}]{P\+\_\+t, }\item[{real, intent(in)}]{P\+\_\+tgt, }\item[{real, intent(in)}]{rho\+\_\+ref, }\item[{real, intent(in)}]{G\+\_\+e, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, intent(out)}]{P\+\_\+b, }\item[{real, intent(out)}]{z\+\_\+out, }\item[{real, intent(in), optional}]{z\+\_\+tol }\end{DoxyParamCaption})} Find the depth at which the reconstructed pressure matches P\+\_\+tgt. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em t\+\_\+t} & Potential temperature at the cell top \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em t\+\_\+b} & Potential temperature at the cell bottom \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+t} & Salinity at the cell top \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+b} & Salinity at the cell bottom \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+t} & Absolute height of top of cell \mbox{[}Z $\sim$$>$ m\mbox{]} (Boussinesq ????)\\ \hline \mbox{\tt in} & {\em z\+\_\+b} & Absolute height of bottom of cell \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+t} & Anomalous pressure of top of cell, relative to g$\ast$rho\+\_\+ref$\ast$z\+\_\+t \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+tgt} & Target pressure at height z\+\_\+out, relative to g$\ast$rho\+\_\+ref$\ast$z\+\_\+out \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em rho\+\_\+ref} & Reference density with which calculation are anomalous to \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]}\\ \hline \mbox{\tt in} & {\em g\+\_\+e} & Gravitational acceleration \mbox{[}L2 Z-\/1 T-\/2 $\sim$$>$ m s-\/2\mbox{]}\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt out} & {\em p\+\_\+b} & Pressure at the bottom of the cell \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]}\\ \hline \mbox{\tt out} & {\em z\+\_\+out} & Absolute depth at which anomalous pressure = p\+\_\+tgt \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+tol} & The tolerance in finding z\+\_\+out \mbox{[}Z $\sim$$>$ m\mbox{]} \\ \hline \end{DoxyParams} Definition at line 1527 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 1527 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: t\_t\textcolor{comment}{ !< Potential temperature at the cell top [degC]} 1528 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: t\_b\textcolor{comment}{ !< Potential temperature at the cell bottom [degC]} 1529 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: s\_t\textcolor{comment}{ !< Salinity at the cell top [ppt]} 1530 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: s\_b\textcolor{comment}{ !< Salinity at the cell bottom [ppt]} 1531 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: z\_t\textcolor{comment}{ !< Absolute height of top of cell [Z ~> m] (Boussinesq ????)} 1532 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: z\_b\textcolor{comment}{ !< Absolute height of bottom of cell [Z ~> m]} 1533 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: p\_t\textcolor{comment}{ !< Anomalous pressure of top of cell, relative} 1534 \textcolor{comment}{ !! to g*rho\_ref*z\_t [R L2 T-2 ~> Pa]} 1535 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: p\_tgt\textcolor{comment}{ !< Target pressure at height z\_out, relative} 1536 \textcolor{comment}{ !! to g*rho\_ref*z\_out [R L2 T-2 ~> Pa]} 1537 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_ref\textcolor{comment}{ !< Reference density with which calculation} 1538 \textcolor{comment}{ !! are anomalous to [R ~> kg m-3]} 1539 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: g\_e\textcolor{comment}{ !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]} 1540 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 1541 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 1542 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(out)} :: p\_b\textcolor{comment}{ !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa]} 1543 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(out)} :: z\_out\textcolor{comment}{ !< Absolute depth at which anomalous pressure = p\_tgt [Z ~> m]} 1544 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: z\_tol\textcolor{comment}{ !< The tolerance in finding z\_out [Z ~> m]} 1545 1546 \textcolor{comment}{! Local variables} 1547 \textcolor{keywordtype}{real} :: dp \textcolor{comment}{! Pressure thickness of the layer [R L2 T-2 ~> Pa]} 1548 \textcolor{keywordtype}{real} :: f\_guess, f\_l, f\_r \textcolor{comment}{! Fractional positions [nondim]} 1549 \textcolor{keywordtype}{real} :: gxrho \textcolor{comment}{! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1]} 1550 \textcolor{keywordtype}{real} :: pa, pa\_left, pa\_right, pa\_tol \textcolor{comment}{! Pressure anomalies, P = integral of g*(rho-rho\_ref) dz [R L2 T-2 ~> Pa]} 1551 \textcolor{keywordtype}{character(len=240)} :: msg 1552 1553 gxrho = g\_e * rho\_ref 1554 1555 \textcolor{comment}{! Anomalous pressure difference across whole cell} 1556 dp = frac\_dp\_at\_pos(t\_t, t\_b, s\_t, s\_b, z\_t, z\_b, rho\_ref, g\_e, 1.0, eos) 1557 1558 p\_b = p\_t + dp \textcolor{comment}{! Anomalous pressure at bottom of cell} 1559 1560 \textcolor{keywordflow}{if} (p\_tgt <= p\_t ) \textcolor{keywordflow}{then} 1561 z\_out = z\_t 1562 \textcolor{keywordflow}{return} 1563 \textcolor{keywordflow}{ endif} 1564 1565 \textcolor{keywordflow}{if} (p\_tgt >= p\_b) \textcolor{keywordflow}{then} 1566 z\_out = z\_b 1567 \textcolor{keywordflow}{return} 1568 \textcolor{keywordflow}{ endif} 1569 1570 f\_l = 0. 1571 pa\_left = p\_t - p\_tgt \textcolor{comment}{! Pa\_left < 0} 1572 f\_r = 1. 1573 pa\_right = p\_b - p\_tgt \textcolor{comment}{! Pa\_right > 0} 1574 pa\_tol = gxrho * 1.0e-5*us%m\_to\_Z 1575 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(z\_tol)) pa\_tol = gxrho * z\_tol 1576 1577 f\_guess = f\_l - pa\_left / (pa\_right - pa\_left) * (f\_r - f\_l) 1578 pa = pa\_right - pa\_left \textcolor{comment}{! To get into iterative loop} 1579 \textcolor{keywordflow}{do} \textcolor{keywordflow}{while} ( abs(pa) > pa\_tol ) 1580 1581 z\_out = z\_t + ( z\_b - z\_t ) * f\_guess 1582 pa = frac\_dp\_at\_pos(t\_t, t\_b, s\_t, s\_b, z\_t, z\_b, rho\_ref, g\_e, f\_guess, eos) - ( p\_tgt - p\_t ) 1583 1584 \textcolor{keywordflow}{if} (papa\_right) \textcolor{keywordflow}{then} 1591 \textcolor{keyword}{write}(msg,*) pa\_left,pa,pa\_right,p\_t-p\_tgt,p\_b-p\_tgt 1592 \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{'find\_depth\_of\_pressure\_in\_cell out of bounds positive: /n'}//msg) 1593 \textcolor{keywordflow}{elseif} (pa>0.) \textcolor{keywordflow}{then} 1594 pa\_right = pa 1595 f\_r = f\_guess 1596 \textcolor{keywordflow}{else} \textcolor{comment}{! Pa == 0} 1597 \textcolor{keywordflow}{return} 1598 \textcolor{keywordflow}{ endif} 1599 f\_guess = f\_l - pa\_left / (pa\_right - pa\_left) * (f\_r - f\_l) 1600 1601 \textcolor{keywordflow}{ enddo} 1602 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_abb8d1f2c24def176a27ef2d30fa373df}\label{namespacemom__density__integrals_abb8d1f2c24def176a27ef2d30fa373df}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!frac\+\_\+dp\+\_\+at\+\_\+pos@{frac\+\_\+dp\+\_\+at\+\_\+pos}} \index{frac\+\_\+dp\+\_\+at\+\_\+pos@{frac\+\_\+dp\+\_\+at\+\_\+pos}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{frac\+\_\+dp\+\_\+at\+\_\+pos()}{frac\_dp\_at\_pos()}} {\footnotesize\ttfamily real function mom\+\_\+density\+\_\+integrals\+::frac\+\_\+dp\+\_\+at\+\_\+pos (\begin{DoxyParamCaption}\item[{real, intent(in)}]{T\+\_\+t, }\item[{real, intent(in)}]{T\+\_\+b, }\item[{real, intent(in)}]{S\+\_\+t, }\item[{real, intent(in)}]{S\+\_\+b, }\item[{real, intent(in)}]{z\+\_\+t, }\item[{real, intent(in)}]{z\+\_\+b, }\item[{real, intent(in)}]{rho\+\_\+ref, }\item[{real, intent(in)}]{G\+\_\+e, }\item[{real, intent(in)}]{pos, }\item[{type(eos\+\_\+type), pointer}]{E\+OS }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}} Returns change in anomalous pressure change from top to non-\/dimensional position pos between z\+\_\+t and z\+\_\+b. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em t\+\_\+t} & Potential temperature at the cell top \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em t\+\_\+b} & Potential temperature at the cell bottom \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+t} & Salinity at the cell top \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+b} & Salinity at the cell bottom \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+t} & The geometric height at the top of the layer \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+b} & The geometric height at the bottom of the layer \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em rho\+\_\+ref} & A mean density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]}, that is subtracted out to reduce the magnitude of each of the integrals.\\ \hline \mbox{\tt in} & {\em g\+\_\+e} & The Earth\textquotesingle{}s gravitational acceleration \mbox{[}L2 Z-\/1 T-\/2 $\sim$$>$ m s-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em pos} & The fractional vertical position, 0 to 1 \mbox{[}nondim\mbox{]}\\ \hline & {\em eos} & Equation of state structure \\ \hline \end{DoxyParams} Definition at line 1609 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 1609 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: t\_t\textcolor{comment}{ !< Potential temperature at the cell top [degC]} 1610 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: t\_b\textcolor{comment}{ !< Potential temperature at the cell bottom [degC]} 1611 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: s\_t\textcolor{comment}{ !< Salinity at the cell top [ppt]} 1612 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: s\_b\textcolor{comment}{ !< Salinity at the cell bottom [ppt]} 1613 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: z\_t\textcolor{comment}{ !< The geometric height at the top of the layer [Z ~> m]} 1614 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: z\_b\textcolor{comment}{ !< The geometric height at the bottom of the layer [Z ~> m]} 1615 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_ref\textcolor{comment}{ !< A mean density [R ~> kg m-3], that is subtracted out to} 1616 \textcolor{comment}{ !! reduce the magnitude of each of the integrals.} 1617 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: g\_e\textcolor{comment}{ !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]} 1618 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: pos\textcolor{comment}{ !< The fractional vertical position, 0 to 1 [nondim]} 1619 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 1620 \textcolor{keywordtype}{real} :: fract\_dp\_at\_pos\textcolor{comment}{ !< The change in pressure from the layer top to} 1621 \textcolor{comment}{ !! fractional position pos [R L2 T-2 ~> Pa]} 1622 \textcolor{comment}{! Local variables} 1623 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_90 = 1.0/90.0 \textcolor{comment}{! A rational constant [nondim]} 1624 \textcolor{keywordtype}{real} :: dz \textcolor{comment}{! Distance from the layer top [Z ~> m]} 1625 \textcolor{keywordtype}{real} :: top\_weight, bottom\_weight \textcolor{comment}{! Fractional weights at quadrature points [nondim]} 1626 \textcolor{keywordtype}{real} :: rho\_ave \textcolor{comment}{! Average density [R ~> kg m-3]} 1627 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(5)} :: t5 \textcolor{comment}{! Temperatures at quadrature points [degC]} 1628 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(5)} :: s5 \textcolor{comment}{! Salinities at quadrature points [ppt]} 1629 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(5)} :: p5 \textcolor{comment}{! Pressures at quadrature points [R L2 T-2 ~> Pa]} 1630 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(5)} :: rho5 \textcolor{comment}{! Densities at quadrature points [R ~> kg m-3]} 1631 \textcolor{keywordtype}{integer} :: n 1632 1633 \textcolor{keywordflow}{do} n=1,5 1634 \textcolor{comment}{! Evaluate density at five quadrature points} 1635 bottom\_weight = 0.25*\textcolor{keywordtype}{real(n-1)} * pos 1636 top\_weight = 1.0 - bottom\_weight 1637 \textcolor{comment}{! Salinity and temperature points are linearly interpolated} 1638 s5(n) = top\_weight * s\_t + bottom\_weight * s\_b 1639 t5(n) = top\_weight * t\_t + bottom\_weight * t\_b 1640 p5(n) = ( top\_weight * z\_t + bottom\_weight * z\_b ) * ( g\_e * rho\_ref ) 1641 \textcolor{keywordflow}{ enddo} 1642 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, rho5, eos) 1643 rho5(:) = rho5(:) \textcolor{comment}{!- rho\_ref ! Work with anomalies relative to rho\_ref} 1644 1645 \textcolor{comment}{! Use Boole's rule to estimate the average density} 1646 rho\_ave = c1\_90*(7.0*(rho5(1)+rho5(5)) + 32.0*(rho5(2)+rho5(4)) + 12.0*rho5(3)) 1647 1648 dz = ( z\_t - z\_b ) * pos 1649 frac\_dp\_at\_pos = g\_e * dz * rho\_ave \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_ac36ae5f4af2d02df0a1adf41b762e017}\label{namespacemom__density__integrals_ac36ae5f4af2d02df0a1adf41b762e017}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!int\+\_\+density\+\_\+dz@{int\+\_\+density\+\_\+dz}} \index{int\+\_\+density\+\_\+dz@{int\+\_\+density\+\_\+dz}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{int\+\_\+density\+\_\+dz()}{int\_density\_dz()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::int\+\_\+density\+\_\+dz (\begin{DoxyParamCaption}\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{T, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{S, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{z\+\_\+t, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{z\+\_\+b, }\item[{real, intent(in)}]{rho\+\_\+ref, }\item[{real, intent(in)}]{rho\+\_\+0, }\item[{real, intent(in)}]{G\+\_\+e, }\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(inout)}]{dpa, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(inout), optional}]{intz\+\_\+dpa, }\item[{real, dimension(szib\+\_\+(hi),szj\+\_\+(hi)), intent(inout), optional}]{intx\+\_\+dpa, }\item[{real, dimension(szi\+\_\+(hi),szjb\+\_\+(hi)), intent(inout), optional}]{inty\+\_\+dpa, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in), optional}]{bathyT, }\item[{real, intent(in), optional}]{dz\+\_\+neglect, }\item[{logical, intent(in), optional}]{use\+Mass\+Wght\+Interp }\end{DoxyParamCaption})} Calls the appropriate subroutine to calculate analytical and nearly-\/analytical integrals in z across layers of pressure anomalies, which are required for calculating the finite-\/volume form pressure accelerations in a Boussinesq model. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em hi} & Ocean horizontal index structures for the arrays\\ \hline \mbox{\tt in} & {\em t} & Potential temperature referenced to the surface \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s} & Salinity \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+t} & Height at the top of the layer in depth units \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+b} & Height at the bottom of the layer \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em rho\+\_\+ref} & A mean density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is subtracted out to reduce the magnitude of each of the integrals.\\ \hline \mbox{\tt in} & {\em rho\+\_\+0} & A density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is used to calculate the pressure (as p$\sim$=-\/z$\ast$rho\+\_\+0$\ast$\+G\+\_\+e) used in the equation of state.\\ \hline \mbox{\tt in} & {\em g\+\_\+e} & The Earth\textquotesingle{}s gravitational acceleration \mbox{[}L2 Z-\/1 T-\/2 $\sim$$>$ m s-\/2\mbox{]} or \mbox{[}m2 Z-\/1 s-\/2 $\sim$$>$ m s-\/2\mbox{]}\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in,out} & {\em dpa} & The change in the pressure anomaly\\ \hline \mbox{\tt in,out} & {\em intz\+\_\+dpa} & The integral through the thickness of the\\ \hline \mbox{\tt in,out} & {\em intx\+\_\+dpa} & The integral in x of the difference between\\ \hline \mbox{\tt in,out} & {\em inty\+\_\+dpa} & The integral in y of the difference between\\ \hline \mbox{\tt in} & {\em bathyt} & The depth of the bathymetry \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em dz\+\_\+neglect} & A minuscule thickness change \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em usemasswghtinterp} & If true, uses mass weighting to interpolate T/S for top and bottom integrals. \\ \hline \end{DoxyParams} Definition at line 42 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 42 \textcolor{keywordtype}{type}(hor\_index\_type), \textcolor{keywordtype}{intent(in)} :: hi\textcolor{comment}{ !< Ocean horizontal index structures for the arrays} 43 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 44 \textcolor{keywordtype}{intent(in)} :: t\textcolor{comment}{ !< Potential temperature referenced to the surface [degC]} 45 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 46 \textcolor{keywordtype}{intent(in)} :: s\textcolor{comment}{ !< Salinity [ppt]} 47 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 48 \textcolor{keywordtype}{intent(in)} :: z\_t\textcolor{comment}{ !< Height at the top of the layer in depth units [Z ~> m]} 49 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 50 \textcolor{keywordtype}{intent(in)} :: z\_b\textcolor{comment}{ !< Height at the bottom of the layer [Z ~> m]} 51 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_ref\textcolor{comment}{ !< A mean density [R ~> kg m-3] or [kg m-3], that is} 52 \textcolor{comment}{ !! subtracted out to reduce the magnitude of each of the} 53 \textcolor{comment}{ !! integrals.} 54 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_0\textcolor{comment}{ !< A density [R ~> kg m-3] or [kg m-3], that is used} 55 \textcolor{comment}{ !! to calculate the pressure (as p~=-z*rho\_0*G\_e)} 56 \textcolor{comment}{ !! used in the equation of state.} 57 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: g\_e\textcolor{comment}{ !< The Earth's gravitational acceleration} 58 \textcolor{comment}{ !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2]} 59 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 60 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 61 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 62 \textcolor{keywordtype}{intent(inout)} :: dpa\textcolor{comment}{ !< The change in the pressure anomaly} 63 \textcolor{comment}{ !! across the layer [R L2 T-2 ~> Pa] or [Pa]} 64 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 65 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intz\_dpa\textcolor{comment}{ !< The integral through the thickness of the} 66 \textcolor{comment}{ !! layer of the pressure anomaly relative to the} 67 \textcolor{comment}{ !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]} 68 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(HI),SZJ\_(HI))}, & 69 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intx\_dpa\textcolor{comment}{ !< The integral in x of the difference between} 70 \textcolor{comment}{ !! the pressure anomaly at the top and bottom of the} 71 \textcolor{comment}{ !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]} 72 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJB\_(HI))}, & 73 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: inty\_dpa\textcolor{comment}{ !< The integral in y of the difference between} 74 \textcolor{comment}{ !! the pressure anomaly at the top and bottom of the} 75 \textcolor{comment}{ !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]} 76 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 77 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: bathyt\textcolor{comment}{ !< The depth of the bathymetry [Z ~> m]} 78 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: dz\_neglect\textcolor{comment}{ !< A minuscule thickness change [Z ~> m]} 79 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: usemasswghtinterp\textcolor{comment}{ !< If true, uses mass weighting to} 80 \textcolor{comment}{ !! interpolate T/S for top and bottom integrals.} 81 82 \textcolor{keywordflow}{if} (eos\_quadrature(eos)) \textcolor{keywordflow}{then} 83 \textcolor{keyword}{call }int\_density\_dz\_generic\_pcm(t, s, z\_t, z\_b, rho\_ref, rho\_0, g\_e, hi, eos, us, dpa, & 84 intz\_dpa, intx\_dpa, inty\_dpa, bathyt, dz\_neglect, usemasswghtinterp) 85 \textcolor{keywordflow}{else} 86 \textcolor{keyword}{call }analytic\_int\_density\_dz(t, s, z\_t, z\_b, rho\_ref, rho\_0, g\_e, hi, eos, dpa, & 87 intz\_dpa, intx\_dpa, inty\_dpa, bathyt, dz\_neglect, usemasswghtinterp) 88 \textcolor{keywordflow}{ endif} 89 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_ad0ff54518bfacf7beebbe4dba687a914}\label{namespacemom__density__integrals_ad0ff54518bfacf7beebbe4dba687a914}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+pcm@{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+pcm}} \index{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+pcm@{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+pcm}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+pcm()}{int\_density\_dz\_generic\_pcm()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+pcm (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{T, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{S, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{z\+\_\+t, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{z\+\_\+b, }\item[{real, intent(in)}]{rho\+\_\+ref, }\item[{real, intent(in)}]{rho\+\_\+0, }\item[{real, intent(in)}]{G\+\_\+e, }\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout)}]{dpa, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intz\+\_\+dpa, }\item[{real, dimension( hi \%isdb\+: hi \%iedb, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intx\+\_\+dpa, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsdb\+: hi \%jedb), intent(inout), optional}]{inty\+\_\+dpa, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in), optional}]{bathyT, }\item[{real, intent(in), optional}]{dz\+\_\+neglect, }\item[{logical, intent(in), optional}]{use\+Mass\+Wght\+Interp }\end{DoxyParamCaption})} Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which are required for calculating the finite-\/volume form pressure accelerations in a Boussinesq model. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em hi} & Horizontal index type for input variables.\\ \hline \mbox{\tt in} & {\em t} & Potential temperature of the layer \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s} & Salinity of the layer \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+t} & Height at the top of the layer in depth units \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em z\+\_\+b} & Height at the bottom of the layer \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em rho\+\_\+ref} & A mean density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is subtracted out to reduce the magnitude of each of the integrals.\\ \hline \mbox{\tt in} & {\em rho\+\_\+0} & A density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is used to calculate the pressure (as p$\sim$=-\/z$\ast$rho\+\_\+0$\ast$\+G\+\_\+e) used in the equation of state.\\ \hline \mbox{\tt in} & {\em g\+\_\+e} & The Earth\textquotesingle{}s gravitational acceleration \mbox{[}L2 Z-\/1 T-\/2 $\sim$$>$ m s-\/2\mbox{]} or \mbox{[}m2 Z-\/1 s-\/2 $\sim$$>$ m s-\/2\mbox{]}\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in,out} & {\em dpa} & The change in the pressure anomaly\\ \hline \mbox{\tt in,out} & {\em intz\+\_\+dpa} & The integral through the thickness of the\\ \hline \mbox{\tt in,out} & {\em intx\+\_\+dpa} & The integral in x of the difference between\\ \hline \mbox{\tt in,out} & {\em inty\+\_\+dpa} & The integral in y of the difference between\\ \hline \mbox{\tt in} & {\em bathyt} & The depth of the bathymetry \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em dz\+\_\+neglect} & A minuscule thickness change \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em usemasswghtinterp} & If true, uses mass weighting to interpolate T/S for top and bottom integrals. \\ \hline \end{DoxyParams} Definition at line 98 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 98 \textcolor{keywordtype}{type}(hor\_index\_type), \textcolor{keywordtype}{intent(in)} :: hi\textcolor{comment}{ !< Horizontal index type for input variables.} 99 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 100 \textcolor{keywordtype}{intent(in)} :: t\textcolor{comment}{ !< Potential temperature of the layer [degC]} 101 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 102 \textcolor{keywordtype}{intent(in)} :: s\textcolor{comment}{ !< Salinity of the layer [ppt]} 103 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 104 \textcolor{keywordtype}{intent(in)} :: z\_t\textcolor{comment}{ !< Height at the top of the layer in depth units [Z ~> m]} 105 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 106 \textcolor{keywordtype}{intent(in)} :: z\_b\textcolor{comment}{ !< Height at the bottom of the layer [Z ~> m]} 107 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_ref\textcolor{comment}{ !< A mean density [R ~> kg m-3] or [kg m-3], that is} 108 \textcolor{comment}{ !! subtracted out to reduce the magnitude} 109 \textcolor{comment}{ !! of each of the integrals.} 110 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_0\textcolor{comment}{ !< A density [R ~> kg m-3] or [kg m-3], that is used} 111 \textcolor{comment}{ !! to calculate the pressure (as p~=-z*rho\_0*G\_e)} 112 \textcolor{comment}{ !! used in the equation of state.} 113 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: g\_e\textcolor{comment}{ !< The Earth's gravitational acceleration} 114 \textcolor{comment}{ !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2]} 115 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 116 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 117 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 118 \textcolor{keywordtype}{intent(inout)} :: dpa\textcolor{comment}{ !< The change in the pressure anomaly} 119 \textcolor{comment}{ !! across the layer [R L2 T-2 ~> Pa] or [Pa]} 120 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 121 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intz\_dpa\textcolor{comment}{ !< The integral through the thickness of the} 122 \textcolor{comment}{ !! layer of the pressure anomaly relative to the} 123 \textcolor{comment}{ !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]} 124 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(HI),SZJ\_(HI))}, & 125 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intx\_dpa\textcolor{comment}{ !< The integral in x of the difference between} 126 \textcolor{comment}{ !! the pressure anomaly at the top and bottom of the} 127 \textcolor{comment}{ !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]} 128 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJB\_(HI))}, & 129 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: inty\_dpa\textcolor{comment}{ !< The integral in y of the difference between} 130 \textcolor{comment}{ !! the pressure anomaly at the top and bottom of the} 131 \textcolor{comment}{ !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]} 132 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 133 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: bathyt\textcolor{comment}{ !< The depth of the bathymetry [Z ~> m]} 134 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: dz\_neglect\textcolor{comment}{ !< A minuscule thickness change [Z ~> m]} 135 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: usemasswghtinterp\textcolor{comment}{ !< If true, uses mass weighting to} 136 \textcolor{comment}{ !! interpolate T/S for top and bottom integrals.} 137 \textcolor{comment}{! Local variables} 138 \textcolor{keywordtype}{real} :: t5(5), s5(5) \textcolor{comment}{! Temperatures and salinities at five quadrature points [degC] and [ppt]} 139 \textcolor{keywordtype}{real} :: p5(5) \textcolor{comment}{! Pressures at five quadrature points, never rescaled from Pa [Pa]} 140 \textcolor{keywordtype}{real} :: r5(5) \textcolor{comment}{! Densities at five quadrature points [R ~> kg m-3] or [kg m-3]} 141 \textcolor{keywordtype}{real} :: rho\_anom \textcolor{comment}{! The depth averaged density anomaly [R ~> kg m-3] or [kg m-3]} 142 \textcolor{keywordtype}{real} :: w\_left, w\_right \textcolor{comment}{! Left and right weights [nondim]} 143 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_90 = 1.0/90.0 \textcolor{comment}{! Rational constants.} 144 \textcolor{keywordtype}{real} :: gxrho \textcolor{comment}{! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]} 145 \textcolor{keywordtype}{real} :: i\_rho \textcolor{comment}{! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]} 146 \textcolor{keywordtype}{real} :: rho\_scale \textcolor{comment}{! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]} 147 \textcolor{keywordtype}{real} :: rho\_ref\_mks \textcolor{comment}{! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]} 148 \textcolor{keywordtype}{real} :: dz \textcolor{comment}{! The layer thickness [Z ~> m]} 149 \textcolor{keywordtype}{real} :: hwght \textcolor{comment}{! A pressure-thickness below topography [Z ~> m]} 150 \textcolor{keywordtype}{real} :: hl, hr \textcolor{comment}{! Pressure-thicknesses of the columns to the left and right [Z ~> m]} 151 \textcolor{keywordtype}{real} :: idenom \textcolor{comment}{! The inverse of the denominator in the weights [Z-2 ~> m-2]} 152 \textcolor{keywordtype}{real} :: hwt\_ll, hwt\_lr \textcolor{comment}{! hWt\_LA is the weighted influence of A on the left column [nondim]} 153 \textcolor{keywordtype}{real} :: hwt\_rl, hwt\_rr \textcolor{comment}{! hWt\_RA is the weighted influence of A on the right column [nondim]} 154 \textcolor{keywordtype}{real} :: wt\_l, wt\_r \textcolor{comment}{! The linear weights of the left and right columns [nondim]} 155 \textcolor{keywordtype}{real} :: wtt\_l, wtt\_r \textcolor{comment}{! The weights for tracers from the left and right columns [nondim]} 156 \textcolor{keywordtype}{real} :: intz(5) \textcolor{comment}{! The gravitational acceleration times the integrals of density} 157 \textcolor{comment}{! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]} 158 \textcolor{keywordtype}{logical} :: do\_massweight \textcolor{comment}{! Indicates whether to do mass weighting.} 159 \textcolor{keywordtype}{integer} :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n 160 161 \textcolor{comment}{! These array bounds work for the indexing convention of the input arrays, but} 162 \textcolor{comment}{! on the computational domain defined for the output arrays.} 163 isq = hi%IscB ; ieq = hi%IecB 164 jsq = hi%JscB ; jeq = hi%JecB 165 is = hi%isc ; ie = hi%iec 166 js = hi%jsc ; je = hi%jec 167 168 rho\_scale = us%kg\_m3\_to\_R 169 gxrho = us%RL2\_T2\_to\_Pa * g\_e * rho\_0 170 rho\_ref\_mks = rho\_ref * us%R\_to\_kg\_m3 171 i\_rho = 1.0 / rho\_0 172 173 do\_massweight = .false. 174 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(usemasswghtinterp)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (usemasswghtinterp) \textcolor{keywordflow}{then} 175 do\_massweight = .true. 176 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{present}(bathyt)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"int\_density\_dz\_generic: "}//& 177 \textcolor{stringliteral}{"bathyT must be present if useMassWghtInterp is present and true."}) 178 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{present}(dz\_neglect)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"int\_density\_dz\_generic: "}//& 179 \textcolor{stringliteral}{"dz\_neglect must be present if useMassWghtInterp is present and true."}) 180 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ endif} 181 182 \textcolor{keywordflow}{do} j=jsq,jeq+1 ; \textcolor{keywordflow}{do} i=isq,ieq+1 183 dz = z\_t(i,j) - z\_b(i,j) 184 \textcolor{keywordflow}{do} n=1,5 185 t5(n) = t(i,j) ; s5(n) = s(i,j) 186 p5(n) = -gxrho*(z\_t(i,j) - 0.25*\textcolor{keywordtype}{real}(n-1)*dz) 187 \textcolor{keywordflow}{ enddo} 188 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 189 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 190 \textcolor{keywordflow}{else} 191 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks) 192 \textcolor{keywordflow}{ endif} 193 194 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 195 rho\_anom = c1\_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) 196 dpa(i,j) = g\_e*dz*rho\_anom 197 \textcolor{comment}{! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of} 198 \textcolor{comment}{! the pressure anomaly.} 199 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intz\_dpa)) intz\_dpa(i,j) = 0.5*g\_e*dz**2 * & 200 (rho\_anom - c1\_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) ) 201 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 202 203 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intx\_dpa)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=isq,ieq 204 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 205 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation of} 206 \textcolor{comment}{! T & S along the top and bottom integrals, akin to thickness weighting.} 207 hwght = 0.0 208 \textcolor{keywordflow}{if} (do\_massweight) & 209 hwght = max(0., -bathyt(i,j)-z\_t(i+1,j), -bathyt(i+1,j)-z\_t(i,j)) 210 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 211 hl = (z\_t(i,j) - z\_b(i,j)) + dz\_neglect 212 hr = (z\_t(i+1,j) - z\_b(i+1,j)) + dz\_neglect 213 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 214 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr ) 215 hwt\_ll = (hwght*hl + hr*hl) * idenom ; hwt\_lr = (hwght*hr) * idenom 216 hwt\_rr = (hwght*hr + hr*hl) * idenom ; hwt\_rl = (hwght*hl) * idenom 217 \textcolor{keywordflow}{else} 218 hwt\_ll = 1.0 ; hwt\_lr = 0.0 ; hwt\_rr = 1.0 ; hwt\_rl = 0.0 219 \textcolor{keywordflow}{ endif} 220 221 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) 222 \textcolor{keywordflow}{do} m=2,4 223 \textcolor{comment}{! T, S, and z are interpolated in the horizontal. The z interpolation} 224 \textcolor{comment}{! is linear, but for T and S it may be thickness weighted.} 225 wt\_l = 0.25*\textcolor{keywordtype}{real(5-m)} ; wt\_r = 1.0-wt\_l 226 wtt\_l = wt\_l*hwt\_ll + wt\_r*hwt\_rl ; wtt\_r = wt\_l*hwt\_lr + wt\_r*hwt\_rr 227 dz = wt\_l*(z\_t(i,j) - z\_b(i,j)) + wt\_r*(z\_t(i+1,j) - z\_b(i+1,j)) 228 t5(1) = wtt\_l*t(i,j) + wtt\_r*t(i+1,j) 229 s5(1) = wtt\_l*s(i,j) + wtt\_r*s(i+1,j) 230 p5(1) = -gxrho*(wt\_l*z\_t(i,j) + wt\_r*z\_t(i+1,j)) 231 \textcolor{keywordflow}{do} n=2,5 232 t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) + gxrho*0.25*dz 233 \textcolor{keywordflow}{ enddo} 234 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 235 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 236 \textcolor{keywordflow}{else} 237 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks) 238 \textcolor{keywordflow}{ endif} 239 240 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 241 intz(m) = g\_e*dz*( c1\_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) 242 \textcolor{keywordflow}{ enddo} 243 \textcolor{comment}{! Use Boole's rule to integrate the bottom pressure anomaly values in x.} 244 intx\_dpa(i,j) = c1\_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & 245 12.0*intz(3)) 246 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 247 248 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(inty\_dpa)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} i=is,ie 249 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 250 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation of} 251 \textcolor{comment}{! T & S along the top and bottom integrals, akin to thickness weighting.} 252 hwght = 0.0 253 \textcolor{keywordflow}{if} (do\_massweight) & 254 hwght = max(0., -bathyt(i,j)-z\_t(i,j+1), -bathyt(i,j+1)-z\_t(i,j)) 255 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 256 hl = (z\_t(i,j) - z\_b(i,j)) + dz\_neglect 257 hr = (z\_t(i,j+1) - z\_b(i,j+1)) + dz\_neglect 258 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 259 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr ) 260 hwt\_ll = (hwght*hl + hr*hl) * idenom ; hwt\_lr = (hwght*hr) * idenom 261 hwt\_rr = (hwght*hr + hr*hl) * idenom ; hwt\_rl = (hwght*hl) * idenom 262 \textcolor{keywordflow}{else} 263 hwt\_ll = 1.0 ; hwt\_lr = 0.0 ; hwt\_rr = 1.0 ; hwt\_rl = 0.0 264 \textcolor{keywordflow}{ endif} 265 266 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) 267 \textcolor{keywordflow}{do} m=2,4 268 \textcolor{comment}{! T, S, and z are interpolated in the horizontal. The z interpolation} 269 \textcolor{comment}{! is linear, but for T and S it may be thickness weighted.} 270 wt\_l = 0.25*\textcolor{keywordtype}{real(5-m)} ; wt\_r = 1.0-wt\_l 271 wtt\_l = wt\_l*hwt\_ll + wt\_r*hwt\_rl ; wtt\_r = wt\_l*hwt\_lr + wt\_r*hwt\_rr 272 dz = wt\_l*(z\_t(i,j) - z\_b(i,j)) + wt\_r*(z\_t(i,j+1) - z\_b(i,j+1)) 273 t5(1) = wtt\_l*t(i,j) + wtt\_r*t(i,j+1) 274 s5(1) = wtt\_l*s(i,j) + wtt\_r*s(i,j+1) 275 p5(1) = -gxrho*(wt\_l*z\_t(i,j) + wt\_r*z\_t(i,j+1)) 276 \textcolor{keywordflow}{do} n=2,5 277 t5(n) = t5(1) ; s5(n) = s5(1) 278 p5(n) = p5(n-1) + gxrho*0.25*dz 279 \textcolor{keywordflow}{ enddo} 280 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 281 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 282 \textcolor{keywordflow}{else} 283 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks) 284 \textcolor{keywordflow}{ endif} 285 286 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 287 intz(m) = g\_e*dz*( c1\_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) 288 \textcolor{keywordflow}{ enddo} 289 \textcolor{comment}{! Use Boole's rule to integrate the values.} 290 inty\_dpa(i,j) = c1\_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & 291 12.0*intz(3)) 292 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_a454d4d62eb599716cc3389fb8aa90b4b}\label{namespacemom__density__integrals_a454d4d62eb599716cc3389fb8aa90b4b}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+plm@{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+plm}} \index{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+plm@{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+plm}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+plm()}{int\_density\_dz\_generic\_plm()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+plm (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{k, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed, gv \%ke), intent(in)}]{T\+\_\+t, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed, gv \%ke), intent(in)}]{T\+\_\+b, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed, gv \%ke), intent(in)}]{S\+\_\+t, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed, gv \%ke), intent(in)}]{S\+\_\+b, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed, gv \%ke+1), intent(in)}]{e, }\item[{real, intent(in)}]{rho\+\_\+ref, }\item[{real, intent(in)}]{rho\+\_\+0, }\item[{real, intent(in)}]{G\+\_\+e, }\item[{real, intent(in)}]{dz\+\_\+subroundoff, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{bathyT, }\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(inout)}]{dpa, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(inout), optional}]{intz\+\_\+dpa, }\item[{real, dimension(szib\+\_\+(hi),szj\+\_\+(hi)), intent(inout), optional}]{intx\+\_\+dpa, }\item[{real, dimension(szi\+\_\+(hi),szjb\+\_\+(hi)), intent(inout), optional}]{inty\+\_\+dpa, }\item[{logical, intent(in), optional}]{use\+Mass\+Wght\+Interp }\end{DoxyParamCaption})} Compute pressure gradient force integrals by quadrature for the case where T and S are linear profiles. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em k} & Layer index to calculate integrals for\\ \hline \mbox{\tt in} & {\em hi} & Ocean horizontal index structures for the input arrays\\ \hline \mbox{\tt in} & {\em gv} & Vertical grid structure\\ \hline \mbox{\tt in} & {\em tv} & Thermodynamic variables\\ \hline \mbox{\tt in} & {\em t\+\_\+t} & Potential temperature at the cell top \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em t\+\_\+b} & Potential temperature at the cell bottom \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+t} & Salinity at the cell top \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+b} & Salinity at the cell bottom \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em e} & Height of interfaces \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em rho\+\_\+ref} & A mean density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is subtracted out to reduce the magnitude of each of the integrals.\\ \hline \mbox{\tt in} & {\em rho\+\_\+0} & A density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is used to calculate the pressure (as p$\sim$=-\/z$\ast$rho\+\_\+0$\ast$\+G\+\_\+e) used in the equation of state.\\ \hline \mbox{\tt in} & {\em g\+\_\+e} & The Earth\textquotesingle{}s gravitational acceleration \mbox{[}L2 Z-\/1 T-\/2 $\sim$$>$ m s-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em dz\+\_\+subroundoff} & A minuscule thickness change \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em bathyt} & The depth of the bathymetry \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in,out} & {\em dpa} & The change in the pressure anomaly across the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]}\\ \hline \mbox{\tt in,out} & {\em intz\+\_\+dpa} & The integral through the thickness of the layer of\\ \hline \mbox{\tt in,out} & {\em intx\+\_\+dpa} & The integral in x of the difference between the\\ \hline \mbox{\tt in,out} & {\em inty\+\_\+dpa} & The integral in y of the difference between the\\ \hline \mbox{\tt in} & {\em usemasswghtinterp} & If true, uses mass weighting to interpolate T/S for top and bottom integrals. \\ \hline \end{DoxyParams} Definition at line 301 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 301 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: k\textcolor{comment}{ !< Layer index to calculate integrals for} 302 \textcolor{keywordtype}{type}(hor\_index\_type), \textcolor{keywordtype}{intent(in)} :: hi\textcolor{comment}{ !< Ocean horizontal index structures for the input arrays} 303 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< Vertical grid structure} 304 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{ !< Thermodynamic variables} 305 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 306 \textcolor{keywordtype}{intent(in)} :: t\_t\textcolor{comment}{ !< Potential temperature at the cell top [degC]} 307 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 308 \textcolor{keywordtype}{intent(in)} :: t\_b\textcolor{comment}{ !< Potential temperature at the cell bottom [degC]} 309 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 310 \textcolor{keywordtype}{intent(in)} :: s\_t\textcolor{comment}{ !< Salinity at the cell top [ppt]} 311 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 312 \textcolor{keywordtype}{intent(in)} :: s\_b\textcolor{comment}{ !< Salinity at the cell bottom [ppt]} 313 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV)+1)}, & 314 \textcolor{keywordtype}{intent(in)} :: e\textcolor{comment}{ !< Height of interfaces [Z ~> m]} 315 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_ref\textcolor{comment}{ !< A mean density [R ~> kg m-3] or [kg m-3], that is subtracted} 316 \textcolor{comment}{ !! out to reduce the magnitude of each of the integrals.} 317 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_0\textcolor{comment}{ !< A density [R ~> kg m-3] or [kg m-3], that is used to calculate} 318 \textcolor{comment}{ !! the pressure (as p~=-z*rho\_0*G\_e) used in the equation of state.} 319 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: g\_e\textcolor{comment}{ !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]} 320 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: dz\_subroundoff\textcolor{comment}{ !< A minuscule thickness change [Z ~> m]} 321 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 322 \textcolor{keywordtype}{intent(in)} :: bathyt\textcolor{comment}{ !< The depth of the bathymetry [Z ~> m]} 323 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 324 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 325 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 326 \textcolor{keywordtype}{intent(inout)} :: dpa\textcolor{comment}{ !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]} 327 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 328 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intz\_dpa\textcolor{comment}{ !< The integral through the thickness of the layer of} 329 \textcolor{comment}{ !! the pressure anomaly relative to the anomaly at the} 330 \textcolor{comment}{ !! top of the layer [R L2 Z T-2 ~> Pa Z]} 331 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(HI),SZJ\_(HI))}, & 332 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intx\_dpa\textcolor{comment}{ !< The integral in x of the difference between the} 333 \textcolor{comment}{ !! pressure anomaly at the top and bottom of the layer} 334 \textcolor{comment}{ !! divided by the x grid spacing [R L2 T-2 ~> Pa]} 335 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJB\_(HI))}, & 336 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: inty\_dpa\textcolor{comment}{ !< The integral in y of the difference between the} 337 \textcolor{comment}{ !! pressure anomaly at the top and bottom of the layer} 338 \textcolor{comment}{ !! divided by the y grid spacing [R L2 T-2 ~> Pa]} 339 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: usemasswghtinterp\textcolor{comment}{ !< If true, uses mass weighting to} 340 \textcolor{comment}{ !! interpolate T/S for top and bottom integrals.} 341 342 \textcolor{comment}{! This subroutine calculates (by numerical quadrature) integrals of} 343 \textcolor{comment}{! pressure anomalies across layers, which are required for calculating the} 344 \textcolor{comment}{! finite-volume form pressure accelerations in a Boussinesq model. The one} 345 \textcolor{comment}{! potentially dodgy assumption here is that rho\_0 is used both in the denominator} 346 \textcolor{comment}{! of the accelerations, and in the pressure used to calculated density (the} 347 \textcolor{comment}{! latter being -z*rho\_0*G\_e). These two uses could be separated if need be.} 348 \textcolor{comment}{!} 349 \textcolor{comment}{! It is assumed that the salinity and temperature profiles are linear in the} 350 \textcolor{comment}{! vertical. The top and bottom values within each layer are provided and} 351 \textcolor{comment}{! a linear interpolation is used to compute intermediate values.} 352 353 \textcolor{comment}{! Local variables} 354 \textcolor{keywordtype}{real} :: t5((5*hi%iscb+1):(5*(hi%iecb+2))) \textcolor{comment}{! Temperatures along a line of subgrid locations [degC]} 355 \textcolor{keywordtype}{real} :: s5((5*hi%iscb+1):(5*(hi%iecb+2))) \textcolor{comment}{! Salinities along a line of subgrid locations [ppt]} 356 \textcolor{keywordtype}{real} :: t25((5*hi%iscb+1):(5*(hi%iecb+2))) \textcolor{comment}{! SGS temperature variance along a line of subgrid locations [degC2]} 357 \textcolor{keywordtype}{real} :: ts5((5*hi%iscb+1):(5*(hi%iecb+2))) \textcolor{comment}{! SGS temp-salt covariance along a line of subgrid locations [degC ppt]} 358 \textcolor{keywordtype}{real} :: s25((5*hi%iscb+1):(5*(hi%iecb+2))) \textcolor{comment}{! SGS salinity variance along a line of subgrid locations [ppt2]} 359 \textcolor{keywordtype}{real} :: p5((5*hi%iscb+1):(5*(hi%iecb+2))) \textcolor{comment}{! Pressures along a line of subgrid locations, never} 360 \textcolor{comment}{! rescaled from Pa [Pa]} 361 \textcolor{keywordtype}{real} :: r5((5*hi%iscb+1):(5*(hi%iecb+2))) \textcolor{comment}{! Densities anomalies along a line of subgrid} 362 \textcolor{comment}{! locations [R ~> kg m-3] or [kg m-3]} 363 \textcolor{keywordtype}{real} :: t15((15*hi%iscb+1):(15*(hi%iecb+1))) \textcolor{comment}{! Temperatures at an array of subgrid locations [degC]} 364 \textcolor{keywordtype}{real} :: s15((15*hi%iscb+1):(15*(hi%iecb+1))) \textcolor{comment}{! Salinities at an array of subgrid locations [ppt]} 365 \textcolor{keywordtype}{real} :: t215((15*hi%iscb+1):(15*(hi%iecb+1))) \textcolor{comment}{! SGS temperature variance along a line of subgrid locations [degC2]} 366 \textcolor{keywordtype}{real} :: ts15((15*hi%iscb+1):(15*(hi%iecb+1))) \textcolor{comment}{! SGS temp-salt covariance along a line of subgrid locations [degC ppt]} 367 \textcolor{keywordtype}{real} :: s215((15*hi%iscb+1):(15*(hi%iecb+1))) \textcolor{comment}{! SGS salinity variance along a line of subgrid locations [ppt2]} 368 \textcolor{keywordtype}{real} :: p15((15*hi%iscb+1):(15*(hi%iecb+1))) \textcolor{comment}{! Pressures at an array of subgrid locations [Pa]} 369 \textcolor{keywordtype}{real} :: r15((15*hi%iscb+1):(15*(hi%iecb+1))) \textcolor{comment}{! Densities at an array of subgrid locations} 370 \textcolor{comment}{! [R ~> kg m-3] or [kg m-3]} 371 \textcolor{keywordtype}{real} :: wt\_t(5), wt\_b(5) \textcolor{comment}{! Top and bottom weights [nondim]} 372 \textcolor{keywordtype}{real} :: rho\_anom \textcolor{comment}{! A density anomaly [R ~> kg m-3] or [kg m-3]} 373 \textcolor{keywordtype}{real} :: w\_left, w\_right \textcolor{comment}{! Left and right weights [nondim]} 374 \textcolor{keywordtype}{real} :: intz(5) \textcolor{comment}{! The gravitational acceleration times the integrals of density} 375 \textcolor{comment}{! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]} 376 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_90 = 1.0/90.0 \textcolor{comment}{! A rational constant [nondim]} 377 \textcolor{keywordtype}{real} :: gxrho \textcolor{comment}{! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]} 378 \textcolor{keywordtype}{real} :: i\_rho \textcolor{comment}{! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]} 379 \textcolor{keywordtype}{real} :: rho\_scale \textcolor{comment}{! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]} 380 \textcolor{keywordtype}{real} :: rho\_ref\_mks \textcolor{comment}{! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]} 381 \textcolor{keywordtype}{real} :: dz(hi%iscb:hi%iecb+1) \textcolor{comment}{! Layer thicknesses at tracer points [Z ~> m]} 382 \textcolor{keywordtype}{real} :: dz\_x(5,hi%iscb:hi%iecb) \textcolor{comment}{! Layer thicknesses along an x-line of subrid locations [Z ~> m]} 383 \textcolor{keywordtype}{real} :: dz\_y(5,hi%isc:hi%iec) \textcolor{comment}{! Layer thicknesses along a y-line of subrid locations [Z ~> m]} 384 \textcolor{keywordtype}{real} :: massweighttoggle \textcolor{comment}{! A non-dimensional toggle factor (0 or 1) [nondim]} 385 \textcolor{keywordtype}{real} :: ttl, tbl, ttr, tbr \textcolor{comment}{! Temperatures at the velocity cell corners [degC]} 386 \textcolor{keywordtype}{real} :: stl, sbl, str, sbr \textcolor{comment}{! Salinities at the velocity cell corners [ppt]} 387 \textcolor{keywordtype}{real} :: hwght \textcolor{comment}{! A topographically limited thicknes weight [Z ~> m]} 388 \textcolor{keywordtype}{real} :: hl, hr \textcolor{comment}{! Thicknesses to the left and right [Z ~> m]} 389 \textcolor{keywordtype}{real} :: idenom \textcolor{comment}{! The denominator of the thickness weight expressions [Z-2 ~> m-2]} 390 \textcolor{keywordtype}{logical} :: use\_stanley\_eos \textcolor{comment}{! True is SGS variance fields exist in tv.} 391 \textcolor{keywordtype}{logical} :: use\_vart, use\_vars, use\_covarts 392 \textcolor{keywordtype}{integer} :: isq, ieq, jsq, jeq, i, j, m, n 393 \textcolor{keywordtype}{integer} :: pos 394 395 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB 396 397 rho\_scale = us%kg\_m3\_to\_R 398 gxrho = us%RL2\_T2\_to\_Pa * g\_e * rho\_0 399 rho\_ref\_mks = rho\_ref * us%R\_to\_kg\_m3 400 i\_rho = 1.0 / rho\_0 401 massweighttoggle = 0. 402 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(usemasswghtinterp)) \textcolor{keywordflow}{then} 403 \textcolor{keywordflow}{if} (usemasswghtinterp) massweighttoggle = 1. 404 \textcolor{keywordflow}{ endif} 405 406 use\_vart = \textcolor{keyword}{associated}(tv%varT) 407 use\_covarts = \textcolor{keyword}{associated}(tv%covarTS) 408 use\_vars = \textcolor{keyword}{associated}(tv%varS) 409 use\_stanley\_eos = use\_vart .or. use\_covarts .or. use\_vars 410 t25(:) = 0. 411 ts5(:) = 0. 412 s25(:) = 0. 413 t215(:) = 0. 414 ts15(:) = 0. 415 s215(:) = 0. 416 417 \textcolor{keywordflow}{do} n = 1, 5 418 wt\_t(n) = 0.25 * \textcolor{keywordtype}{real}(5-n) 419 wt\_b(n) = 1.0 - wt\_t(n) 420 \textcolor{keywordflow}{ enddo} 421 422 \textcolor{comment}{! 1. Compute vertical integrals} 423 \textcolor{keywordflow}{do} j=jsq,jeq+1 424 \textcolor{keywordflow}{do} i = isq,ieq+1 425 dz(i) = e(i,j,k) - e(i,j,k+1) 426 \textcolor{keywordflow}{do} n=1,5 427 p5(i*5+n) = -gxrho*(e(i,j,k) - 0.25*\textcolor{keywordtype}{real}(n-1)*dz(i)) 428 \textcolor{comment}{! Salinity and temperature points are linearly interpolated} 429 s5(i*5+n) = wt\_t(n) * s\_t(i,j,k) + wt\_b(n) * s\_b(i,j,k) 430 t5(i*5+n) = wt\_t(n) * t\_t(i,j,k) + wt\_b(n) * t\_b(i,j,k) 431 \textcolor{keywordflow}{ enddo} 432 \textcolor{keywordflow}{if} (use\_vart) t25(i*5+1:i*5+5) = tv%varT(i,j,k) 433 \textcolor{keywordflow}{if} (use\_covarts) ts5(i*5+1:i*5+5) = tv%covarTS(i,j,k) 434 \textcolor{keywordflow}{if} (use\_vars) s25(i*5+1:i*5+5) = tv%varS(i,j,k) 435 \textcolor{keywordflow}{ enddo} 436 \textcolor{keywordflow}{if} (use\_stanley\_eos) \textcolor{keywordflow}{then} 437 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 438 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, t25, ts5, s25, r5, 1, (ieq-isq+2)*5, eos, & 439 rho\_ref=rho\_ref\_mks, scale=rho\_scale) 440 \textcolor{keywordflow}{else} 441 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, t25, ts5, s25, r5, 1, (ieq-isq+2)*5, eos, & 442 rho\_ref=rho\_ref\_mks) 443 \textcolor{keywordflow}{ endif} 444 \textcolor{keywordflow}{else} 445 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 446 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho\_ref=rho\_ref\_mks, & 447 scale=rho\_scale) 448 \textcolor{keywordflow}{else} 449 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho\_ref=rho\_ref\_mks) 450 \textcolor{keywordflow}{ endif} 451 \textcolor{keywordflow}{ endif} 452 453 \textcolor{keywordflow}{do} i=isq,ieq+1 454 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 455 rho\_anom = c1\_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3)) 456 dpa(i,j) = g\_e*dz(i)*rho\_anom 457 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intz\_dpa)) \textcolor{keywordflow}{then} 458 \textcolor{comment}{! Use a Boole's-rule-like fifth-order accurate estimate of} 459 \textcolor{comment}{! the double integral of the pressure anomaly.} 460 intz\_dpa(i,j) = 0.5*g\_e*dz(i)**2 * & 461 (rho\_anom - c1\_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) ) 462 \textcolor{keywordflow}{ endif} 463 \textcolor{keywordflow}{ enddo} 464 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loops on j} 465 466 \textcolor{comment}{! 2. Compute horizontal integrals in the x direction} 467 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intx\_dpa)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=hi%jsc,hi%jec 468 \textcolor{keywordflow}{do} i=isq,ieq 469 \textcolor{comment}{! Corner values of T and S} 470 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 471 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation} 472 \textcolor{comment}{! of T,S along the top and bottom integrals, almost like thickness} 473 \textcolor{comment}{! weighting.} 474 \textcolor{comment}{! Note: To work in terrain following coordinates we could offset} 475 \textcolor{comment}{! this distance by the layer thickness to replicate other models.} 476 hwght = massweighttoggle * & 477 max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k)) 478 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 479 hl = (e(i,j,k) - e(i,j,k+1)) + dz\_subroundoff 480 hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz\_subroundoff 481 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 482 idenom = 1./( hwght*(hr + hl) + hl*hr ) 483 ttl = ( (hwght*hr)*t\_t(i+1,j,k) + (hwght*hl + hr*hl)*t\_t(i,j,k) ) * idenom 484 ttr = ( (hwght*hl)*t\_t(i,j,k) + (hwght*hr + hr*hl)*t\_t(i+1,j,k) ) * idenom 485 tbl = ( (hwght*hr)*t\_b(i+1,j,k) + (hwght*hl + hr*hl)*t\_b(i,j,k) ) * idenom 486 tbr = ( (hwght*hl)*t\_b(i,j,k) + (hwght*hr + hr*hl)*t\_b(i+1,j,k) ) * idenom 487 stl = ( (hwght*hr)*s\_t(i+1,j,k) + (hwght*hl + hr*hl)*s\_t(i,j,k) ) * idenom 488 str = ( (hwght*hl)*s\_t(i,j,k) + (hwght*hr + hr*hl)*s\_t(i+1,j,k) ) * idenom 489 sbl = ( (hwght*hr)*s\_b(i+1,j,k) + (hwght*hl + hr*hl)*s\_b(i,j,k) ) * idenom 490 sbr = ( (hwght*hl)*s\_b(i,j,k) + (hwght*hr + hr*hl)*s\_b(i+1,j,k) ) * idenom 491 \textcolor{keywordflow}{else} 492 ttl = t\_t(i,j,k); tbl = t\_b(i,j,k); ttr = t\_t(i+1,j,k); tbr = t\_b(i+1,j,k) 493 stl = s\_t(i,j,k); sbl = s\_b(i,j,k); str = s\_t(i+1,j,k); sbr = s\_b(i+1,j,k) 494 \textcolor{keywordflow}{ endif} 495 496 \textcolor{keywordflow}{do} m=2,4 497 w\_left = wt\_t(m) ; w\_right = wt\_b(m) 498 dz\_x(m,i) = w\_left*(e(i,j,k) - e(i,j,k+1)) + w\_right*(e(i+1,j,k) - e(i+1,j,k+1)) 499 500 \textcolor{comment}{! Salinity and temperature points are linearly interpolated in} 501 \textcolor{comment}{! the horizontal. The subscript (1) refers to the top value in} 502 \textcolor{comment}{! the vertical profile while subscript (5) refers to the bottom} 503 \textcolor{comment}{! value in the vertical profile.} 504 pos = i*15+(m-2)*5 505 t15(pos+1) = w\_left*ttl + w\_right*ttr 506 t15(pos+5) = w\_left*tbl + w\_right*tbr 507 508 s15(pos+1) = w\_left*stl + w\_right*str 509 s15(pos+5) = w\_left*sbl + w\_right*sbr 510 511 p15(pos+1) = -gxrho*(w\_left*e(i,j,k) + w\_right*e(i+1,j,k)) 512 513 \textcolor{comment}{! Pressure} 514 \textcolor{keywordflow}{do} n=2,5 515 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz\_x(m,i) 516 \textcolor{keywordflow}{ enddo} 517 518 \textcolor{comment}{! Salinity and temperature (linear interpolation in the vertical)} 519 \textcolor{keywordflow}{do} n=2,4 520 s15(pos+n) = wt\_t(n) * s15(pos+1) + wt\_b(n) * s15(pos+5) 521 t15(pos+n) = wt\_t(n) * t15(pos+1) + wt\_b(n) * t15(pos+5) 522 \textcolor{keywordflow}{ enddo} 523 \textcolor{keywordflow}{if} (use\_vart) t215(pos+1:pos+5) = w\_left*tv%varT(i,j,k) + w\_right*tv%varT(i+1,j,k) 524 \textcolor{keywordflow}{if} (use\_covarts) ts15(pos+1:pos+5) = w\_left*tv%covarTS(i,j,k) + w\_right*tv%covarTS(i+1,j,k) 525 \textcolor{keywordflow}{if} (use\_vars) s215(pos+1:pos+5) = w\_left*tv%varS(i,j,k) + w\_right*tv%varS(i+1,j,k) 526 \textcolor{keywordflow}{ enddo} 527 \textcolor{keywordflow}{ enddo} 528 529 \textcolor{keywordflow}{if} (use\_stanley\_eos) \textcolor{keywordflow}{then} 530 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 531 \textcolor{keyword}{call }calculate\_density(t15, s15, p15, t215, ts15, s215, r15, 1, 15*(ieq-isq+1), eos, & 532 rho\_ref=rho\_ref\_mks, scale=rho\_scale) 533 \textcolor{keywordflow}{else} 534 \textcolor{keyword}{call }calculate\_density(t15, s15, p15, t215, ts15, s215, r15, 1, 15*(ieq-isq+1), eos, & 535 rho\_ref=rho\_ref\_mks) 536 \textcolor{keywordflow}{ endif} 537 \textcolor{keywordflow}{else} 538 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 539 \textcolor{keyword}{call }calculate\_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho\_ref=rho\_ref\_mks, & 540 scale=rho\_scale) 541 \textcolor{keywordflow}{else} 542 \textcolor{keyword}{call }calculate\_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho\_ref=rho\_ref\_mks) 543 \textcolor{keywordflow}{ endif} 544 \textcolor{keywordflow}{ endif} 545 546 \textcolor{keywordflow}{do} i=isq,ieq 547 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) 548 549 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 550 \textcolor{keywordflow}{do} m = 2,4 551 pos = i*15+(m-2)*5 552 intz(m) = g\_e*dz\_x(m,i)*( c1\_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + & 553 12.0*r15(pos+3))) 554 \textcolor{keywordflow}{ enddo} 555 \textcolor{comment}{! Use Boole's rule to integrate the bottom pressure anomaly values in x.} 556 intx\_dpa(i,j) = c1\_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & 557 12.0*intz(3)) 558 \textcolor{keywordflow}{ enddo} 559 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 560 561 \textcolor{comment}{! 3. Compute horizontal integrals in the y direction} 562 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(inty\_dpa)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=jsq,jeq 563 \textcolor{keywordflow}{do} i=hi%isc,hi%iec 564 \textcolor{comment}{! Corner values of T and S} 565 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 566 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation} 567 \textcolor{comment}{! of T,S along the top and bottom integrals, almost like thickness} 568 \textcolor{comment}{! weighting.} 569 \textcolor{comment}{! Note: To work in terrain following coordinates we could offset} 570 \textcolor{comment}{! this distance by the layer thickness to replicate other models.} 571 hwght = massweighttoggle * & 572 max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k)) 573 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 574 hl = (e(i,j,k) - e(i,j,k+1)) + dz\_subroundoff 575 hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz\_subroundoff 576 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 577 idenom = 1./( hwght*(hr + hl) + hl*hr ) 578 ttl = ( (hwght*hr)*t\_t(i,j+1,k) + (hwght*hl + hr*hl)*t\_t(i,j,k) ) * idenom 579 ttr = ( (hwght*hl)*t\_t(i,j,k) + (hwght*hr + hr*hl)*t\_t(i,j+1,k) ) * idenom 580 tbl = ( (hwght*hr)*t\_b(i,j+1,k) + (hwght*hl + hr*hl)*t\_b(i,j,k) ) * idenom 581 tbr = ( (hwght*hl)*t\_b(i,j,k) + (hwght*hr + hr*hl)*t\_b(i,j+1,k) ) * idenom 582 stl = ( (hwght*hr)*s\_t(i,j+1,k) + (hwght*hl + hr*hl)*s\_t(i,j,k) ) * idenom 583 str = ( (hwght*hl)*s\_t(i,j,k) + (hwght*hr + hr*hl)*s\_t(i,j+1,k) ) * idenom 584 sbl = ( (hwght*hr)*s\_b(i,j+1,k) + (hwght*hl + hr*hl)*s\_b(i,j,k) ) * idenom 585 sbr = ( (hwght*hl)*s\_b(i,j,k) + (hwght*hr + hr*hl)*s\_b(i,j+1,k) ) * idenom 586 \textcolor{keywordflow}{else} 587 ttl = t\_t(i,j,k); tbl = t\_b(i,j,k); ttr = t\_t(i,j+1,k); tbr = t\_b(i,j+1,k) 588 stl = s\_t(i,j,k); sbl = s\_b(i,j,k); str = s\_t(i,j+1,k); sbr = s\_b(i,j+1,k) 589 \textcolor{keywordflow}{ endif} 590 591 \textcolor{keywordflow}{do} m=2,4 592 w\_left = wt\_t(m) ; w\_right = wt\_b(m) 593 dz\_y(m,i) = w\_left*(e(i,j,k) - e(i,j,k+1)) + w\_right*(e(i,j+1,k) - e(i,j+1,k+1)) 594 595 \textcolor{comment}{! Salinity and temperature points are linearly interpolated in} 596 \textcolor{comment}{! the horizontal. The subscript (1) refers to the top value in} 597 \textcolor{comment}{! the vertical profile while subscript (5) refers to the bottom} 598 \textcolor{comment}{! value in the vertical profile.} 599 pos = i*15+(m-2)*5 600 t15(pos+1) = w\_left*ttl + w\_right*ttr 601 t15(pos+5) = w\_left*tbl + w\_right*tbr 602 603 s15(pos+1) = w\_left*stl + w\_right*str 604 s15(pos+5) = w\_left*sbl + w\_right*sbr 605 606 p15(pos+1) = -gxrho*(w\_left*e(i,j,k) + w\_right*e(i,j+1,k)) 607 608 \textcolor{comment}{! Pressure} 609 \textcolor{keywordflow}{do} n=2,5 610 p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz\_y(m,i) 611 \textcolor{keywordflow}{ enddo} 612 613 \textcolor{comment}{! Salinity and temperature (linear interpolation in the vertical)} 614 \textcolor{keywordflow}{do} n=2,4 615 s15(pos+n) = wt\_t(n) * s15(pos+1) + wt\_b(n) * s15(pos+5) 616 t15(pos+n) = wt\_t(n) * t15(pos+1) + wt\_b(n) * t15(pos+5) 617 \textcolor{keywordflow}{ enddo} 618 \textcolor{keywordflow}{if} (use\_vart) t215(pos+1:pos+5) = w\_left*tv%varT(i,j,k) + w\_right*tv%varT(i,j+1,k) 619 \textcolor{keywordflow}{if} (use\_covarts) ts15(pos+1:pos+5) = w\_left*tv%covarTS(i,j,k) + w\_right*tv%covarTS(i,j+1,k) 620 \textcolor{keywordflow}{if} (use\_vars) s215(pos+1:pos+5) = w\_left*tv%varS(i,j,k) + w\_right*tv%varS(i,j+1,k) 621 \textcolor{keywordflow}{ enddo} 622 \textcolor{keywordflow}{ enddo} 623 624 \textcolor{keywordflow}{if} (use\_stanley\_eos) \textcolor{keywordflow}{then} 625 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 626 \textcolor{keyword}{call }calculate\_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), & 627 t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), & 628 r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, & 629 rho\_ref=rho\_ref\_mks, scale=rho\_scale) 630 \textcolor{keywordflow}{else} 631 \textcolor{keyword}{call }calculate\_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), & 632 t215(15*hi%isc+1:), ts15(15*hi%isc+1:), s215(15*hi%isc+1:), & 633 r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, rho\_ref=rho\_ref\_mks) 634 \textcolor{keywordflow}{ endif} 635 \textcolor{keywordflow}{else} 636 \textcolor{keywordflow}{if} (rho\_scale /= 1.0) \textcolor{keywordflow}{then} 637 \textcolor{keyword}{call }calculate\_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), & 638 r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, & 639 rho\_ref=rho\_ref\_mks, scale=rho\_scale) 640 \textcolor{keywordflow}{else} 641 \textcolor{keyword}{call }calculate\_density(t15(15*hi%isc+1:), s15(15*hi%isc+1:), p15(15*hi%isc+1:), & 642 r15(15*hi%isc+1:), 1, 15*(hi%iec-hi%isc+1), eos, rho\_ref=rho\_ref\_mks) 643 \textcolor{keywordflow}{ endif} 644 \textcolor{keywordflow}{ endif} 645 646 \textcolor{keywordflow}{do} i=hi%isc,hi%iec 647 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) 648 649 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 650 \textcolor{keywordflow}{do} m = 2,4 651 pos = i*15+(m-2)*5 652 intz(m) = g\_e*dz\_y(m,i)*( c1\_90*(7.0*(r15(pos+1)+r15(pos+5)) + & 653 32.0*(r15(pos+2)+r15(pos+4)) + & 654 12.0*r15(pos+3))) 655 \textcolor{keywordflow}{ enddo} 656 \textcolor{comment}{! Use Boole's rule to integrate the values.} 657 inty\_dpa(i,j) = c1\_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + & 658 12.0*intz(3)) 659 \textcolor{keywordflow}{ enddo} 660 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 661 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_a3bf090e5cbb58811b3dd0ee4de80f999}\label{namespacemom__density__integrals_a3bf090e5cbb58811b3dd0ee4de80f999}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+ppm@{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+ppm}} \index{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+ppm@{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+ppm}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+ppm()}{int\_density\_dz\_generic\_ppm()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::int\+\_\+density\+\_\+dz\+\_\+generic\+\_\+ppm (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{k, }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi),szk\+\_\+(gv)), intent(in)}]{T\+\_\+t, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi),szk\+\_\+(gv)), intent(in)}]{T\+\_\+b, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi),szk\+\_\+(gv)), intent(in)}]{S\+\_\+t, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi),szk\+\_\+(gv)), intent(in)}]{S\+\_\+b, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi),szk\+\_\+(gv)+1), intent(in)}]{e, }\item[{real, intent(in)}]{rho\+\_\+ref, }\item[{real, intent(in)}]{rho\+\_\+0, }\item[{real, intent(in)}]{G\+\_\+e, }\item[{real, intent(in)}]{dz\+\_\+subroundoff, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{bathyT, }\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI, }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout)}]{dpa, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intz\+\_\+dpa, }\item[{real, dimension( hi \%isdb\+: hi \%iedb, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intx\+\_\+dpa, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsdb\+: hi \%jedb), intent(inout), optional}]{inty\+\_\+dpa, }\item[{logical, intent(in), optional}]{use\+Mass\+Wght\+Interp }\end{DoxyParamCaption})} Compute pressure gradient force integrals for layer \char`\"{}k\char`\"{} and the case where T and S are parabolic profiles. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em k} & Layer index to calculate integrals for\\ \hline \mbox{\tt in} & {\em hi} & Ocean horizontal index structures for the input arrays\\ \hline \mbox{\tt in} & {\em gv} & Vertical grid structure\\ \hline \mbox{\tt in} & {\em tv} & Thermodynamic variables\\ \hline \mbox{\tt in} & {\em t\+\_\+t} & Potential temperature at the cell top \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em t\+\_\+b} & Potential temperature at the cell bottom \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+t} & Salinity at the cell top \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+b} & Salinity at the cell bottom \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em e} & Height of interfaces \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em rho\+\_\+ref} & A mean density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is subtracted out to reduce the magnitude of each of the integrals.\\ \hline \mbox{\tt in} & {\em rho\+\_\+0} & A density \mbox{[}R $\sim$$>$ kg m-\/3\mbox{]} or \mbox{[}kg m-\/3\mbox{]}, that is used to calculate the pressure (as p$\sim$=-\/z$\ast$rho\+\_\+0$\ast$\+G\+\_\+e) used in the equation of state.\\ \hline \mbox{\tt in} & {\em g\+\_\+e} & The Earth\textquotesingle{}s gravitational acceleration \mbox{[}m s-\/2\mbox{]}\\ \hline \mbox{\tt in} & {\em dz\+\_\+subroundoff} & A minuscule thickness change \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline \mbox{\tt in} & {\em bathyt} & The depth of the bathymetry \mbox{[}Z $\sim$$>$ m\mbox{]}\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in,out} & {\em dpa} & The change in the pressure anomaly across the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]}\\ \hline \mbox{\tt in,out} & {\em intz\+\_\+dpa} & The integral through the thickness of the layer of\\ \hline \mbox{\tt in,out} & {\em intx\+\_\+dpa} & The integral in x of the difference between the\\ \hline \mbox{\tt in,out} & {\em inty\+\_\+dpa} & The integral in y of the difference between the\\ \hline \mbox{\tt in} & {\em usemasswghtinterp} & If true, uses mass weighting to interpolate T/S for top and bottom integrals. \\ \hline \end{DoxyParams} Definition at line 670 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 670 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{intent(in)} :: k\textcolor{comment}{ !< Layer index to calculate integrals for} 671 \textcolor{keywordtype}{type}(hor\_index\_type), \textcolor{keywordtype}{intent(in)} :: hi\textcolor{comment}{ !< Ocean horizontal index structures for the input arrays} 672 \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{ !< Vertical grid structure} 673 \textcolor{keywordtype}{type}(thermo\_var\_ptrs), \textcolor{keywordtype}{intent(in)} :: tv\textcolor{comment}{ !< Thermodynamic variables} 674 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 675 \textcolor{keywordtype}{intent(in)} :: t\_t\textcolor{comment}{ !< Potential temperature at the cell top [degC]} 676 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 677 \textcolor{keywordtype}{intent(in)} :: t\_b\textcolor{comment}{ !< Potential temperature at the cell bottom [degC]} 678 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 679 \textcolor{keywordtype}{intent(in)} :: s\_t\textcolor{comment}{ !< Salinity at the cell top [ppt]} 680 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV))}, & 681 \textcolor{keywordtype}{intent(in)} :: s\_b\textcolor{comment}{ !< Salinity at the cell bottom [ppt]} 682 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI),SZK\_(GV)+1)}, & 683 \textcolor{keywordtype}{intent(in)} :: e\textcolor{comment}{ !< Height of interfaces [Z ~> m]} 684 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_ref\textcolor{comment}{ !< A mean density [R ~> kg m-3] or [kg m-3], that is} 685 \textcolor{comment}{ !! subtracted out to reduce the magnitude of each of the integrals.} 686 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: rho\_0\textcolor{comment}{ !< A density [R ~> kg m-3] or [kg m-3], that is used to calculate} 687 \textcolor{comment}{ !! the pressure (as p~=-z*rho\_0*G\_e) used in the equation of state.} 688 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: g\_e\textcolor{comment}{ !< The Earth's gravitational acceleration [m s-2]} 689 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: dz\_subroundoff\textcolor{comment}{ !< A minuscule thickness change [Z ~> m]} 690 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 691 \textcolor{keywordtype}{intent(in)} :: bathyt\textcolor{comment}{ !< The depth of the bathymetry [Z ~> m]} 692 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 693 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 694 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 695 \textcolor{keywordtype}{intent(inout)} :: dpa\textcolor{comment}{ !< The change in the pressure anomaly across the layer [R L2 T-2 ~> Pa]} 696 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 697 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intz\_dpa\textcolor{comment}{ !< The integral through the thickness of the layer of} 698 \textcolor{comment}{ !! the pressure anomaly relative to the anomaly at the} 699 \textcolor{comment}{ !! top of the layer [R L2 Z T-2 ~> Pa m]} 700 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(HI),SZJ\_(HI))}, & 701 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intx\_dpa\textcolor{comment}{ !< The integral in x of the difference between the} 702 \textcolor{comment}{ !! pressure anomaly at the top and bottom of the layer} 703 \textcolor{comment}{ !! divided by the x grid spacing [R L2 T-2 ~> Pa]} 704 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJB\_(HI))}, & 705 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: inty\_dpa\textcolor{comment}{ !< The integral in y of the difference between the} 706 \textcolor{comment}{ !! pressure anomaly at the top and bottom of the layer} 707 \textcolor{comment}{ !! divided by the y grid spacing [R L2 T-2 ~> Pa]} 708 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: usemasswghtinterp\textcolor{comment}{ !< If true, uses mass weighting to} 709 \textcolor{comment}{ !! interpolate T/S for top and bottom integrals.} 710 711 \textcolor{comment}{! This subroutine calculates (by numerical quadrature) integrals of} 712 \textcolor{comment}{! pressure anomalies across layers, which are required for calculating the} 713 \textcolor{comment}{! finite-volume form pressure accelerations in a Boussinesq model. The one} 714 \textcolor{comment}{! potentially dodgy assumption here is that rho\_0 is used both in the denominator} 715 \textcolor{comment}{! of the accelerations, and in the pressure used to calculated density (the} 716 \textcolor{comment}{! latter being -z*rho\_0*G\_e). These two uses could be separated if need be.} 717 \textcolor{comment}{!} 718 \textcolor{comment}{! It is assumed that the salinity and temperature profiles are parabolic in the} 719 \textcolor{comment}{! vertical. The top and bottom values within each layer are provided and} 720 \textcolor{comment}{! a parabolic interpolation is used to compute intermediate values.} 721 722 \textcolor{comment}{! Local variables} 723 \textcolor{keywordtype}{real} :: t5(5) \textcolor{comment}{! Temperatures along a line of subgrid locations [degC]} 724 \textcolor{keywordtype}{real} :: s5(5) \textcolor{comment}{! Salinities along a line of subgrid locations [ppt]} 725 \textcolor{keywordtype}{real} :: t25(5) \textcolor{comment}{! SGS temperature variance along a line of subgrid locations [degC2]} 726 \textcolor{keywordtype}{real} :: ts5(5) \textcolor{comment}{! SGS temperature-salinity covariance along a line of subgrid locations [degC ppt]} 727 \textcolor{keywordtype}{real} :: s25(5) \textcolor{comment}{! SGS salinity variance along a line of subgrid locations [ppt2]} 728 \textcolor{keywordtype}{real} :: p5(5) \textcolor{comment}{! Pressures at five quadrature points, never rescaled from Pa [Pa]} 729 \textcolor{keywordtype}{real} :: r5(5) \textcolor{comment}{! Density anomalies from rho\_ref at quadrature points [R ~> kg m-3] or [kg m-3]} 730 \textcolor{keywordtype}{real} :: wt\_t(5), wt\_b(5) \textcolor{comment}{! Top and bottom weights [nondim]} 731 \textcolor{keywordtype}{real} :: rho\_anom \textcolor{comment}{! The integrated density anomaly [R ~> kg m-3] or [kg m-3]} 732 \textcolor{keywordtype}{real} :: w\_left, w\_right \textcolor{comment}{! Left and right weights [nondim]} 733 \textcolor{keywordtype}{real} :: intz(5) \textcolor{comment}{! The gravitational acceleration times the integrals of density} 734 \textcolor{comment}{! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]} 735 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_90 = 1.0/90.0 \textcolor{comment}{! Rational constants.} 736 \textcolor{keywordtype}{real} :: gxrho \textcolor{comment}{! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]} 737 \textcolor{keywordtype}{real} :: i\_rho \textcolor{comment}{! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1]} 738 \textcolor{keywordtype}{real} :: rho\_scale \textcolor{comment}{! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1]} 739 \textcolor{keywordtype}{real} :: rho\_ref\_mks \textcolor{comment}{! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]} 740 \textcolor{keywordtype}{real} :: dz \textcolor{comment}{! Layer thicknesses at tracer points [Z ~> m]} 741 \textcolor{keywordtype}{real} :: massweighttoggle \textcolor{comment}{! A non-dimensional toggle factor (0 or 1) [nondim]} 742 \textcolor{keywordtype}{real} :: ttl, tbl, tml, ttr, tbr, tmr \textcolor{comment}{! Temperatures at the velocity cell corners [degC]} 743 \textcolor{keywordtype}{real} :: stl, sbl, sml, str, sbr, smr \textcolor{comment}{! Salinities at the velocity cell corners [ppt]} 744 \textcolor{keywordtype}{real} :: s6 \textcolor{comment}{! PPM curvature coefficient for S [ppt]} 745 \textcolor{keywordtype}{real} :: t6 \textcolor{comment}{! PPM curvature coefficient for T [degC]} 746 \textcolor{keywordtype}{real} :: t\_top, t\_mn, t\_bot \textcolor{comment}{! Left edge, cell mean and right edge values used in PPM reconstructions of T} 747 \textcolor{keywordtype}{real} :: s\_top, s\_mn, s\_bot \textcolor{comment}{! Left edge, cell mean and right edge values used in PPM reconstructions of S} 748 \textcolor{keywordtype}{real} :: hwght \textcolor{comment}{! A topographically limited thicknes weight [Z ~> m]} 749 \textcolor{keywordtype}{real} :: hl, hr \textcolor{comment}{! Thicknesses to the left and right [Z ~> m]} 750 \textcolor{keywordtype}{real} :: idenom \textcolor{comment}{! The denominator of the thickness weight expressions [Z-2 ~> m-2]} 751 \textcolor{keywordtype}{integer} :: isq, ieq, jsq, jeq, i, j, m, n 752 \textcolor{keywordtype}{logical} :: use\_ppm \textcolor{comment}{! If false, assume zero curvature in reconstruction, i.e. PLM} 753 \textcolor{keywordtype}{logical} :: use\_stanley\_eos \textcolor{comment}{! True is SGS variance fields exist in tv.} 754 \textcolor{keywordtype}{logical} :: use\_vart, use\_vars, use\_covarts 755 756 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB 757 758 rho\_scale = us%kg\_m3\_to\_R 759 gxrho = us%RL2\_T2\_to\_Pa * g\_e * rho\_0 760 rho\_ref\_mks = rho\_ref * us%R\_to\_kg\_m3 761 i\_rho = 1.0 / rho\_0 762 massweighttoggle = 0. 763 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(usemasswghtinterp)) \textcolor{keywordflow}{then} 764 \textcolor{keywordflow}{if} (usemasswghtinterp) massweighttoggle = 1. 765 \textcolor{keywordflow}{ endif} 766 767 \textcolor{comment}{! In event PPM calculation is bypassed with use\_PPM=False} 768 s6 = 0. 769 t6 = 0. 770 use\_ppm = .true. \textcolor{comment}{! This is a place-holder to allow later re-use of this function} 771 772 use\_vart = \textcolor{keyword}{associated}(tv%varT) 773 use\_covarts = \textcolor{keyword}{associated}(tv%covarTS) 774 use\_vars = \textcolor{keyword}{associated}(tv%varS) 775 use\_stanley\_eos = use\_vart .or. use\_covarts .or. use\_vars 776 t25(:) = 0. 777 ts5(:) = 0. 778 s25(:) = 0. 779 780 \textcolor{keywordflow}{do} n = 1, 5 781 wt\_t(n) = 0.25 * \textcolor{keywordtype}{real}(5-n) 782 wt\_b(n) = 1.0 - wt\_t(n) 783 \textcolor{keywordflow}{ enddo} 784 785 \textcolor{comment}{! 1. Compute vertical integrals} 786 \textcolor{keywordflow}{do} j=jsq,jeq+1 ; \textcolor{keywordflow}{do} i=isq,ieq+1 787 \textcolor{keywordflow}{if} (use\_ppm) \textcolor{keywordflow}{then} 788 \textcolor{comment}{! Curvature coefficient of the parabolas} 789 s6 = 3.0 * ( 2.0*tv%S(i,j,k) - ( s\_t(i,j,k) + s\_b(i,j,k) ) ) 790 t6 = 3.0 * ( 2.0*tv%T(i,j,k) - ( t\_t(i,j,k) + t\_b(i,j,k) ) ) 791 \textcolor{keywordflow}{ endif} 792 dz = e(i,j,k) - e(i,j,k+1) 793 \textcolor{keywordflow}{do} n=1,5 794 p5(n) = -gxrho*(e(i,j,k) - 0.25*\textcolor{keywordtype}{real}(n-1)*dz) 795 \textcolor{comment}{! Salinity and temperature points are reconstructed with PPM} 796 s5(n) = wt\_t(n) * s\_t(i,j,k) + wt\_b(n) * ( s\_b(i,j,k) + s6 * wt\_t(n) ) 797 t5(n) = wt\_t(n) * t\_t(i,j,k) + wt\_b(n) * ( t\_b(i,j,k) + t6 * wt\_t(n) ) 798 \textcolor{keywordflow}{ enddo} 799 \textcolor{keywordflow}{if} (use\_stanley\_eos) \textcolor{keywordflow}{then} 800 \textcolor{keywordflow}{if} (use\_vart) t25(:) = tv%varT(i,j,k) 801 \textcolor{keywordflow}{if} (use\_covarts) ts5(:) = tv%covarTS(i,j,k) 802 \textcolor{keywordflow}{if} (use\_vars) s25(:) = tv%varS(i,j,k) 803 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, t25, ts5, s25, r5, & 804 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 805 \textcolor{keywordflow}{else} 806 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 807 \textcolor{keywordflow}{ endif} 808 809 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 810 rho\_anom = c1\_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) 811 dpa(i,j) = g\_e*dz*rho\_anom 812 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intz\_dpa)) \textcolor{keywordflow}{then} 813 \textcolor{comment}{! Use a Boole's-rule-like fifth-order accurate estimate of} 814 \textcolor{comment}{! the double integral of the pressure anomaly.} 815 intz\_dpa(i,j) = 0.5*g\_e*dz**2 * & 816 (rho\_anom - c1\_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) ) 817 \textcolor{keywordflow}{ endif} 818 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loops on j and i} 819 820 \textcolor{comment}{! 2. Compute horizontal integrals in the x direction} 821 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intx\_dpa)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=hi%jsc,hi%jec ; \textcolor{keywordflow}{do} i=isq,ieq 822 \textcolor{comment}{! Corner values of T and S} 823 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 824 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation} 825 \textcolor{comment}{! of T,S along the top and bottom integrals, almost like thickness} 826 \textcolor{comment}{! weighting.} 827 \textcolor{comment}{! Note: To work in terrain following coordinates we could offset} 828 \textcolor{comment}{! this distance by the layer thickness to replicate other models.} 829 hwght = massweighttoggle * & 830 max(0., -bathyt(i,j)-e(i+1,j,k), -bathyt(i+1,j)-e(i,j,k)) 831 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 832 hl = (e(i,j,k) - e(i,j,k+1)) + dz\_subroundoff 833 hr = (e(i+1,j,k) - e(i+1,j,k+1)) + dz\_subroundoff 834 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 835 idenom = 1./( hwght*(hr + hl) + hl*hr ) 836 ttl = ( (hwght*hr)*t\_t(i+1,j,k) + (hwght*hl + hr*hl)*t\_t(i,j,k) ) * idenom 837 tbl = ( (hwght*hr)*t\_b(i+1,j,k) + (hwght*hl + hr*hl)*t\_b(i,j,k) ) * idenom 838 tml = ( (hwght*hr)*tv%T(i+1,j,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom 839 ttr = ( (hwght*hl)*t\_t(i,j,k) + (hwght*hr + hr*hl)*t\_t(i+1,j,k) ) * idenom 840 tbr = ( (hwght*hl)*t\_b(i,j,k) + (hwght*hr + hr*hl)*t\_b(i+1,j,k) ) * idenom 841 tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i+1,j,k) ) * idenom 842 stl = ( (hwght*hr)*s\_t(i+1,j,k) + (hwght*hl + hr*hl)*s\_t(i,j,k) ) * idenom 843 sbl = ( (hwght*hr)*s\_b(i+1,j,k) + (hwght*hl + hr*hl)*s\_b(i,j,k) ) * idenom 844 sml = ( (hwght*hr)*tv%S(i+1,j,k) + (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom 845 str = ( (hwght*hl)*s\_t(i,j,k) + (hwght*hr + hr*hl)*s\_t(i+1,j,k) ) * idenom 846 sbr = ( (hwght*hl)*s\_b(i,j,k) + (hwght*hr + hr*hl)*s\_b(i+1,j,k) ) * idenom 847 smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i+1,j,k) ) * idenom 848 \textcolor{keywordflow}{else} 849 ttl = t\_t(i,j,k); tbl = t\_b(i,j,k); ttr = t\_t(i+1,j,k); tbr = t\_b(i+1,j,k) 850 tml = tv%T(i,j,k); tmr = tv%T(i+1,j,k) 851 stl = s\_t(i,j,k); sbl = s\_b(i,j,k); str = s\_t(i+1,j,k); sbr = s\_b(i+1,j,k) 852 sml = tv%S(i,j,k); smr = tv%S(i+1,j,k) 853 \textcolor{keywordflow}{ endif} 854 855 \textcolor{keywordflow}{do} m=2,4 856 w\_left = wt\_t(m) ; w\_right = wt\_b(m) 857 858 \textcolor{comment}{! Salinity and temperature points are linearly interpolated in} 859 \textcolor{comment}{! the horizontal. The subscript (1) refers to the top value in} 860 \textcolor{comment}{! the vertical profile while subscript (5) refers to the bottom} 861 \textcolor{comment}{! value in the vertical profile.} 862 t\_top = w\_left*ttl + w\_right*ttr 863 t\_mn = w\_left*tml + w\_right*tmr 864 t\_bot = w\_left*tbl + w\_right*tbr 865 866 s\_top = w\_left*stl + w\_right*str 867 s\_mn = w\_left*sml + w\_right*smr 868 s\_bot = w\_left*sbl + w\_right*sbr 869 870 \textcolor{comment}{! Pressure} 871 dz = w\_left*(e(i,j,k) - e(i,j,k+1)) + w\_right*(e(i+1,j,k) - e(i+1,j,k+1)) 872 p5(1) = -gxrho*(w\_left*e(i,j,k) + w\_right*e(i+1,j,k)) 873 \textcolor{keywordflow}{do} n=2,5 874 p5(n) = p5(n-1) + gxrho*0.25*dz 875 \textcolor{keywordflow}{ enddo} 876 877 \textcolor{comment}{! Parabolic reconstructions in the vertical for T and S} 878 \textcolor{keywordflow}{if} (use\_ppm) \textcolor{keywordflow}{then} 879 \textcolor{comment}{! Coefficients of the parabolas} 880 s6 = 3.0 * ( 2.0*s\_mn - ( s\_top + s\_bot ) ) 881 t6 = 3.0 * ( 2.0*t\_mn - ( t\_top + t\_bot ) ) 882 \textcolor{keywordflow}{ endif} 883 \textcolor{keywordflow}{do} n=1,5 884 s5(n) = wt\_t(n) * s\_top + wt\_b(n) * ( s\_bot + s6 * wt\_t(n) ) 885 t5(n) = wt\_t(n) * t\_top + wt\_b(n) * ( t\_bot + t6 * wt\_t(n) ) 886 \textcolor{keywordflow}{ enddo} 887 888 \textcolor{keywordflow}{if} (use\_stanley\_eos) \textcolor{keywordflow}{then} 889 \textcolor{keywordflow}{if} (use\_vart) t25(:) = w\_left*tv%varT(i,j,k) + w\_right*tv%varT(i+1,j,k) 890 \textcolor{keywordflow}{if} (use\_covarts) ts5(:) = w\_left*tv%covarTS(i,j,k) + w\_right*tv%covarTS(i+1,j,k) 891 \textcolor{keywordflow}{if} (use\_vars) s25(:) = w\_left*tv%varS(i,j,k) + w\_right*tv%varS(i+1,j,k) 892 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, t25, ts5, s25, r5, & 893 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 894 \textcolor{keywordflow}{else} 895 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 896 \textcolor{keywordflow}{ endif} 897 898 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 899 intz(m) = g\_e*dz*( c1\_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) ) 900 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! m} 901 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) 902 903 \textcolor{comment}{! Use Boole's rule to integrate the bottom pressure anomaly values in x.} 904 intx\_dpa(i,j) = c1\_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3)) 905 906 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 907 908 \textcolor{comment}{! 3. Compute horizontal integrals in the y direction} 909 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(inty\_dpa)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} i=hi%isc,hi%iec 910 \textcolor{comment}{! Corner values of T and S} 911 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 912 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation} 913 \textcolor{comment}{! of T,S along the top and bottom integrals, almost like thickness} 914 \textcolor{comment}{! weighting.} 915 \textcolor{comment}{! Note: To work in terrain following coordinates we could offset} 916 \textcolor{comment}{! this distance by the layer thickness to replicate other models.} 917 hwght = massweighttoggle * & 918 max(0., -bathyt(i,j)-e(i,j+1,k), -bathyt(i,j+1)-e(i,j,k)) 919 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 920 hl = (e(i,j,k) - e(i,j,k+1)) + dz\_subroundoff 921 hr = (e(i,j+1,k) - e(i,j+1,k+1)) + dz\_subroundoff 922 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 923 idenom = 1./( hwght*(hr + hl) + hl*hr ) 924 ttl = ( (hwght*hr)*t\_t(i,j+1,k) + (hwght*hl + hr*hl)*t\_t(i,j,k) ) * idenom 925 tbl = ( (hwght*hr)*t\_b(i,j+1,k) + (hwght*hl + hr*hl)*t\_b(i,j,k) ) * idenom 926 tml = ( (hwght*hr)*tv%T(i,j+1,k)+ (hwght*hl + hr*hl)*tv%T(i,j,k) ) * idenom 927 ttr = ( (hwght*hl)*t\_t(i,j,k) + (hwght*hr + hr*hl)*t\_t(i,j+1,k) ) * idenom 928 tbr = ( (hwght*hl)*t\_b(i,j,k) + (hwght*hr + hr*hl)*t\_b(i,j+1,k) ) * idenom 929 tmr = ( (hwght*hl)*tv%T(i,j,k) + (hwght*hr + hr*hl)*tv%T(i,j+1,k) ) * idenom 930 stl = ( (hwght*hr)*s\_t(i,j+1,k) + (hwght*hl + hr*hl)*s\_t(i,j,k) ) * idenom 931 sbl = ( (hwght*hr)*s\_b(i,j+1,k) + (hwght*hl + hr*hl)*s\_b(i,j,k) ) * idenom 932 sml = ( (hwght*hr)*tv%S(i,j+1,k)+ (hwght*hl + hr*hl)*tv%S(i,j,k) ) * idenom 933 str = ( (hwght*hl)*s\_t(i,j,k) + (hwght*hr + hr*hl)*s\_t(i,j+1,k) ) * idenom 934 sbr = ( (hwght*hl)*s\_b(i,j,k) + (hwght*hr + hr*hl)*s\_b(i,j+1,k) ) * idenom 935 smr = ( (hwght*hl)*tv%S(i,j,k) + (hwght*hr + hr*hl)*tv%S(i,j+1,k) ) * idenom 936 \textcolor{keywordflow}{else} 937 ttl = t\_t(i,j,k); tbl = t\_b(i,j,k); ttr = t\_t(i,j+1,k); tbr = t\_b(i,j+1,k) 938 tml = tv%T(i,j,k); tmr = tv%T(i,j+1,k) 939 stl = s\_t(i,j,k); sbl = s\_b(i,j,k); str = s\_t(i,j+1,k); sbr = s\_b(i,j+1,k) 940 sml = tv%S(i,j,k); smr = tv%S(i,j+1,k) 941 \textcolor{keywordflow}{ endif} 942 943 \textcolor{keywordflow}{do} m=2,4 944 w\_left = wt\_t(m) ; w\_right = wt\_b(m) 945 946 \textcolor{comment}{! Salinity and temperature points are linearly interpolated in} 947 \textcolor{comment}{! the horizontal. The subscript (1) refers to the top value in} 948 \textcolor{comment}{! the vertical profile while subscript (5) refers to the bottom} 949 \textcolor{comment}{! value in the vertical profile.} 950 t\_top = w\_left*ttl + w\_right*ttr 951 t\_mn = w\_left*tml + w\_right*tmr 952 t\_bot = w\_left*tbl + w\_right*tbr 953 954 s\_top = w\_left*stl + w\_right*str 955 s\_mn = w\_left*sml + w\_right*smr 956 s\_bot = w\_left*sbl + w\_right*sbr 957 958 \textcolor{comment}{! Pressure} 959 dz = w\_left*(e(i,j,k) - e(i,j,k+1)) + w\_right*(e(i,j+1,k) - e(i,j+1,k+1)) 960 p5(1) = -gxrho*(w\_left*e(i,j,k) + w\_right*e(i,j+1,k)) 961 \textcolor{keywordflow}{do} n=2,5 962 p5(n) = p5(n-1) + gxrho*0.25*dz 963 \textcolor{keywordflow}{ enddo} 964 965 \textcolor{comment}{! Parabolic reconstructions in the vertical for T and S} 966 \textcolor{keywordflow}{if} (use\_ppm) \textcolor{keywordflow}{then} 967 \textcolor{comment}{! Coefficients of the parabolas} 968 s6 = 3.0 * ( 2.0*s\_mn - ( s\_top + s\_bot ) ) 969 t6 = 3.0 * ( 2.0*t\_mn - ( t\_top + t\_bot ) ) 970 \textcolor{keywordflow}{ endif} 971 \textcolor{keywordflow}{do} n=1,5 972 s5(n) = wt\_t(n) * s\_top + wt\_b(n) * ( s\_bot + s6 * wt\_t(n) ) 973 t5(n) = wt\_t(n) * t\_top + wt\_b(n) * ( t\_bot + t6 * wt\_t(n) ) 974 \textcolor{keywordflow}{ enddo} 975 976 \textcolor{keywordflow}{if} (use\_stanley\_eos) \textcolor{keywordflow}{then} 977 \textcolor{keywordflow}{if} (use\_vart) t25(:) = w\_left*tv%varT(i,j,k) + w\_right*tv%varT(i,j+1,k) 978 \textcolor{keywordflow}{if} (use\_covarts) ts5(:) = w\_left*tv%covarTS(i,j,k) + w\_right*tv%covarTS(i,j+1,k) 979 \textcolor{keywordflow}{if} (use\_vars) s25(:) = w\_left*tv%varS(i,j,k) + w\_right*tv%varS(i,j+1,k) 980 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, t25, ts5, s25, r5, & 981 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 982 \textcolor{keywordflow}{else} 983 \textcolor{keyword}{call }calculate\_density(t5, s5, p5, r5, 1, 5, eos, rho\_ref=rho\_ref\_mks, scale=rho\_scale) 984 \textcolor{keywordflow}{ endif} 985 986 \textcolor{comment}{! Use Boole's rule to estimate the pressure anomaly change.} 987 intz(m) = g\_e*dz*( c1\_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) ) 988 \textcolor{keywordflow}{ enddo} \textcolor{comment}{! m} 989 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) 990 991 \textcolor{comment}{! Use Boole's rule to integrate the bottom pressure anomaly values in y.} 992 inty\_dpa(i,j) = c1\_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3)) 993 994 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 995 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_a86ebdfebaaea2ac0339dcabed482dd11}\label{namespacemom__density__integrals_a86ebdfebaaea2ac0339dcabed482dd11}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+pcm@{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+pcm}} \index{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+pcm@{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+pcm}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+pcm()}{int\_spec\_vol\_dp\_generic\_pcm()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+pcm (\begin{DoxyParamCaption}\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{T, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{S, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{p\+\_\+t, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{p\+\_\+b, }\item[{real, intent(in)}]{alpha\+\_\+ref, }\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout)}]{dza, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intp\+\_\+dza, }\item[{real, dimension( hi \%isdb\+: hi \%iedb, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intx\+\_\+dza, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsdb\+: hi \%jedb), intent(inout), optional}]{inty\+\_\+dza, }\item[{integer, intent(in), optional}]{halo\+\_\+size, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in), optional}]{bathyP, }\item[{real, intent(in), optional}]{d\+P\+\_\+neglect, }\item[{logical, intent(in), optional}]{use\+Mass\+Wght\+Interp }\end{DoxyParamCaption})} This subroutine calculates integrals of specific volume anomalies in pressure across layers, which are required for calculating the finite-\/volume form pressure accelerations in a non-\/\+Boussinesq model. There are essentially no free assumptions, apart from the use of Boole\textquotesingle{}s rule quadrature to do the integrals. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em hi} & A horizontal index type structure.\\ \hline \mbox{\tt in} & {\em t} & Potential temperature of the layer \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s} & Salinity of the layer \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+t} & Pressure atop the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+b} & Pressure below the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em alpha\+\_\+ref} & A mean specific volume that is subtracted out to reduce the magnitude of each of the integrals \mbox{[}R-\/1 $\sim$$>$ m3 kg-\/1\mbox{]} The calculation is mathematically identical with different values of alpha\+\_\+ref, but alpha\+\_\+ref alters the effects of roundoff, and answers do change.\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in,out} & {\em dza} & The change in the geopotential anomaly\\ \hline \mbox{\tt in,out} & {\em intp\+\_\+dza} & The integral in pressure through the layer of\\ \hline \mbox{\tt in,out} & {\em intx\+\_\+dza} & The integral in x of the difference between\\ \hline \mbox{\tt in,out} & {\em inty\+\_\+dza} & The integral in y of the difference between\\ \hline \mbox{\tt in} & {\em halo\+\_\+size} & The width of halo points on which to calculate dza.\\ \hline \mbox{\tt in} & {\em bathyp} & The pressure at the bathymetry \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em dp\+\_\+neglect} & A minuscule pressure change with the same units as p\+\_\+t \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em usemasswghtinterp} & If true, uses mass weighting to interpolate T/S for top and bottom integrals. \\ \hline \end{DoxyParams} Definition at line 1065 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 1065 \textcolor{keywordtype}{type}(hor\_index\_type), \textcolor{keywordtype}{intent(in)} :: hi\textcolor{comment}{ !< A horizontal index type structure.} 1066 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1067 \textcolor{keywordtype}{intent(in)} :: t\textcolor{comment}{ !< Potential temperature of the layer [degC]} 1068 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1069 \textcolor{keywordtype}{intent(in)} :: s\textcolor{comment}{ !< Salinity of the layer [ppt]} 1070 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1071 \textcolor{keywordtype}{intent(in)} :: p\_t\textcolor{comment}{ !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]} 1072 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1073 \textcolor{keywordtype}{intent(in)} :: p\_b\textcolor{comment}{ !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]} 1074 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: alpha\_ref\textcolor{comment}{ !< A mean specific volume that is subtracted out} 1075 \textcolor{comment}{ !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]} 1076 \textcolor{comment}{ !! The calculation is mathematically identical with different values of} 1077 \textcolor{comment}{ !! alpha\_ref, but alpha\_ref alters the effects of roundoff, and} 1078 \textcolor{comment}{ !! answers do change.} 1079 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 1080 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 1081 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1082 \textcolor{keywordtype}{intent(inout)} :: dza\textcolor{comment}{ !< The change in the geopotential anomaly} 1083 \textcolor{comment}{ !! across the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1084 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1085 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intp\_dza\textcolor{comment}{ !< The integral in pressure through the layer of} 1086 \textcolor{comment}{ !! the geopotential anomaly relative to the anomaly at the bottom of the} 1087 \textcolor{comment}{ !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]} 1088 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(HI),SZJ\_(HI))}, & 1089 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intx\_dza\textcolor{comment}{ !< The integral in x of the difference between} 1090 \textcolor{comment}{ !! the geopotential anomaly at the top and bottom of the layer divided} 1091 \textcolor{comment}{ !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1092 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJB\_(HI))}, & 1093 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: inty\_dza\textcolor{comment}{ !< The integral in y of the difference between} 1094 \textcolor{comment}{ !! the geopotential anomaly at the top and bottom of the layer divided} 1095 \textcolor{comment}{ !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1096 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: halo\_size\textcolor{comment}{ !< The width of halo points on which to calculate dza.} 1097 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1098 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: bathyp\textcolor{comment}{ !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]} 1099 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: dp\_neglect\textcolor{comment}{ !< A minuscule pressure change with} 1100 \textcolor{comment}{ !! the same units as p\_t [R L2 T-2 ~> Pa] or [Pa]} 1101 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: usemasswghtinterp\textcolor{comment}{ !< If true, uses mass weighting} 1102 \textcolor{comment}{ !! to interpolate T/S for top and bottom integrals.} 1103 1104 \textcolor{comment}{! This subroutine calculates analytical and nearly-analytical integrals in} 1105 \textcolor{comment}{! pressure across layers of geopotential anomalies, which are required for} 1106 \textcolor{comment}{! calculating the finite-volume form pressure accelerations in a non-Boussinesq} 1107 \textcolor{comment}{! model. There are essentially no free assumptions, apart from the use of} 1108 \textcolor{comment}{! Boole's rule to do the horizontal integrals, and from a truncation in the} 1109 \textcolor{comment}{! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.} 1110 1111 \textcolor{comment}{! Local variables} 1112 \textcolor{keywordtype}{real} :: t5(5) \textcolor{comment}{! Temperatures at five quadrature points [degC]} 1113 \textcolor{keywordtype}{real} :: s5(5) \textcolor{comment}{! Salinities at five quadrature points [ppt]} 1114 \textcolor{keywordtype}{real} :: p5(5) \textcolor{comment}{! Pressures at five quadrature points, scaled back to Pa if necessary [Pa]} 1115 \textcolor{keywordtype}{real} :: a5(5) \textcolor{comment}{! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]} 1116 \textcolor{keywordtype}{real} :: alpha\_anom \textcolor{comment}{! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]} 1117 \textcolor{keywordtype}{real} :: dp \textcolor{comment}{! The pressure change through a layer [R L2 T-2 ~> Pa]} 1118 \textcolor{keywordtype}{real} :: hwght \textcolor{comment}{! A pressure-thickness below topography [R L2 T-2 ~> Pa]} 1119 \textcolor{keywordtype}{real} :: hl, hr \textcolor{comment}{! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]} 1120 \textcolor{keywordtype}{real} :: alpha\_ref\_mks \textcolor{comment}{! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1]} 1121 \textcolor{keywordtype}{real} :: idenom \textcolor{comment}{! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]} 1122 \textcolor{keywordtype}{real} :: hwt\_ll, hwt\_lr \textcolor{comment}{! hWt\_LA is the weighted influence of A on the left column [nondim]} 1123 \textcolor{keywordtype}{real} :: hwt\_rl, hwt\_rr \textcolor{comment}{! hWt\_RA is the weighted influence of A on the right column [nondim]} 1124 \textcolor{keywordtype}{real} :: wt\_l, wt\_r \textcolor{comment}{! The linear weights of the left and right columns [nondim]} 1125 \textcolor{keywordtype}{real} :: wtt\_l, wtt\_r \textcolor{comment}{! The weights for tracers from the left and right columns [nondim]} 1126 \textcolor{keywordtype}{real} :: intp(5) \textcolor{comment}{! The integrals of specific volume with pressure at the} 1127 \textcolor{comment}{! 5 sub-column locations [L2 T-2 ~> m2 s-2]} 1128 \textcolor{keywordtype}{real} :: rl2\_t2\_to\_pa \textcolor{comment}{! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1]} 1129 \textcolor{keywordtype}{real} :: sv\_scale \textcolor{comment}{! A multiplicative factor by which to scale specific} 1130 \textcolor{comment}{! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]} 1131 \textcolor{keywordtype}{logical} :: do\_massweight \textcolor{comment}{! Indicates whether to do mass weighting.} 1132 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_90 = 1.0/90.0 \textcolor{comment}{! A rational constant.} 1133 \textcolor{keywordtype}{integer} :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, n, halo 1134 1135 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB 1136 halo = 0 ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(halo\_size)) halo = max(halo\_size,0) 1137 ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo 1138 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intx\_dza)) \textcolor{keywordflow}{then} ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh);\textcolor{keywordflow}{ endif} 1139 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(inty\_dza)) \textcolor{keywordflow}{then} ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh);\textcolor{keywordflow}{ endif} 1140 1141 sv\_scale = us%R\_to\_kg\_m3 1142 rl2\_t2\_to\_pa = us%RL2\_T2\_to\_Pa 1143 alpha\_ref\_mks = alpha\_ref * us%kg\_m3\_to\_R 1144 1145 do\_massweight = .false. 1146 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(usemasswghtinterp)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (usemasswghtinterp) \textcolor{keywordflow}{then} 1147 do\_massweight = .true. 1148 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{present}(bathyp)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"int\_spec\_vol\_dp\_generic: "}//& 1149 \textcolor{stringliteral}{"bathyP must be present if useMassWghtInterp is present and true."}) 1150 \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{present}(dp\_neglect)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"int\_spec\_vol\_dp\_generic: "}//& 1151 \textcolor{stringliteral}{"dP\_neglect must be present if useMassWghtInterp is present and true."}) 1152 \textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ endif} 1153 1154 \textcolor{keywordflow}{do} j=jsh,jeh ; \textcolor{keywordflow}{do} i=ish,ieh 1155 dp = p\_b(i,j) - p\_t(i,j) 1156 \textcolor{keywordflow}{do} n=1,5 1157 t5(n) = t(i,j) ; s5(n) = s(i,j) 1158 p5(n) = rl2\_t2\_to\_pa * (p\_b(i,j) - 0.25*\textcolor{keywordtype}{real}(n-1)*dp) 1159 \textcolor{keywordflow}{ enddo} 1160 1161 \textcolor{keywordflow}{if} (sv\_scale /= 1.0) \textcolor{keywordflow}{then} 1162 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks, scale=sv\_scale) 1163 \textcolor{keywordflow}{else} 1164 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks) 1165 \textcolor{keywordflow}{ endif} 1166 1167 \textcolor{comment}{! Use Boole's rule to estimate the interface height anomaly change.} 1168 alpha\_anom = c1\_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3)) 1169 dza(i,j) = dp*alpha\_anom 1170 \textcolor{comment}{! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of} 1171 \textcolor{comment}{! the interface height anomaly.} 1172 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intp\_dza)) intp\_dza(i,j) = 0.5*dp**2 * & 1173 (alpha\_anom - c1\_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) ) 1174 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1175 1176 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intx\_dza)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=hi%jsc,hi%jec ; \textcolor{keywordflow}{do} i=isq,ieq 1177 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 1178 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation of} 1179 \textcolor{comment}{! T & S along the top and bottom integrals, akin to thickness weighting.} 1180 hwght = 0.0 1181 \textcolor{keywordflow}{if} (do\_massweight) & 1182 hwght = max(0., bathyp(i,j)-p\_t(i+1,j), bathyp(i+1,j)-p\_t(i,j)) 1183 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 1184 hl = (p\_b(i,j) - p\_t(i,j)) + dp\_neglect 1185 hr = (p\_b(i+1,j) - p\_t(i+1,j)) + dp\_neglect 1186 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 1187 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr ) 1188 hwt\_ll = (hwght*hl + hr*hl) * idenom ; hwt\_lr = (hwght*hr) * idenom 1189 hwt\_rr = (hwght*hr + hr*hl) * idenom ; hwt\_rl = (hwght*hl) * idenom 1190 \textcolor{keywordflow}{else} 1191 hwt\_ll = 1.0 ; hwt\_lr = 0.0 ; hwt\_rr = 1.0 ; hwt\_rl = 0.0 1192 \textcolor{keywordflow}{ endif} 1193 1194 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) 1195 \textcolor{keywordflow}{do} m=2,4 1196 wt\_l = 0.25*\textcolor{keywordtype}{real(5-m)} ; wt\_r = 1.0-wt\_l 1197 wtt\_l = wt\_l*hwt\_ll + wt\_r*hwt\_rl ; wtt\_r = wt\_l*hwt\_lr + wt\_r*hwt\_rr 1198 1199 \textcolor{comment}{! T, S, and p are interpolated in the horizontal. The p interpolation} 1200 \textcolor{comment}{! is linear, but for T and S it may be thickness weighted.} 1201 p5(1) = rl2\_t2\_to\_pa * (wt\_l*p\_b(i,j) + wt\_r*p\_b(i+1,j)) 1202 dp = wt\_l*(p\_b(i,j) - p\_t(i,j)) + wt\_r*(p\_b(i+1,j) - p\_t(i+1,j)) 1203 t5(1) = wtt\_l*t(i,j) + wtt\_r*t(i+1,j) 1204 s5(1) = wtt\_l*s(i,j) + wtt\_r*s(i+1,j) 1205 1206 \textcolor{keywordflow}{do} n=2,5 1207 t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - rl2\_t2\_to\_pa * 0.25*dp 1208 \textcolor{keywordflow}{ enddo} 1209 \textcolor{keywordflow}{if} (sv\_scale /= 1.0) \textcolor{keywordflow}{then} 1210 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks, scale=sv\_scale) 1211 \textcolor{keywordflow}{else} 1212 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks) 1213 \textcolor{keywordflow}{ endif} 1214 1215 \textcolor{comment}{! Use Boole's rule to estimate the interface height anomaly change.} 1216 intp(m) = dp*( c1\_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & 1217 12.0*a5(3))) 1218 \textcolor{keywordflow}{ enddo} 1219 \textcolor{comment}{! Use Boole's rule to integrate the interface height anomaly values in x.} 1220 intx\_dza(i,j) = c1\_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + & 1221 12.0*intp(3)) 1222 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 1223 1224 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(inty\_dza)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} i=hi%isc,hi%iec 1225 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 1226 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation of} 1227 \textcolor{comment}{! T & S along the top and bottom integrals, akin to thickness weighting.} 1228 hwght = 0.0 1229 \textcolor{keywordflow}{if} (do\_massweight) & 1230 hwght = max(0., bathyp(i,j)-p\_t(i,j+1), bathyp(i,j+1)-p\_t(i,j)) 1231 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 1232 hl = (p\_b(i,j) - p\_t(i,j)) + dp\_neglect 1233 hr = (p\_b(i,j+1) - p\_t(i,j+1)) + dp\_neglect 1234 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 1235 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr ) 1236 hwt\_ll = (hwght*hl + hr*hl) * idenom ; hwt\_lr = (hwght*hr) * idenom 1237 hwt\_rr = (hwght*hr + hr*hl) * idenom ; hwt\_rl = (hwght*hl) * idenom 1238 \textcolor{keywordflow}{else} 1239 hwt\_ll = 1.0 ; hwt\_lr = 0.0 ; hwt\_rr = 1.0 ; hwt\_rl = 0.0 1240 \textcolor{keywordflow}{ endif} 1241 1242 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) 1243 \textcolor{keywordflow}{do} m=2,4 1244 wt\_l = 0.25*\textcolor{keywordtype}{real(5-m)} ; wt\_r = 1.0-wt\_l 1245 wtt\_l = wt\_l*hwt\_ll + wt\_r*hwt\_rl ; wtt\_r = wt\_l*hwt\_lr + wt\_r*hwt\_rr 1246 1247 \textcolor{comment}{! T, S, and p are interpolated in the horizontal. The p interpolation} 1248 \textcolor{comment}{! is linear, but for T and S it may be thickness weighted.} 1249 p5(1) = rl2\_t2\_to\_pa * (wt\_l*p\_b(i,j) + wt\_r*p\_b(i,j+1)) 1250 dp = wt\_l*(p\_b(i,j) - p\_t(i,j)) + wt\_r*(p\_b(i,j+1) - p\_t(i,j+1)) 1251 t5(1) = wtt\_l*t(i,j) + wtt\_r*t(i,j+1) 1252 s5(1) = wtt\_l*s(i,j) + wtt\_r*s(i,j+1) 1253 \textcolor{keywordflow}{do} n=2,5 1254 t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = rl2\_t2\_to\_pa * (p5(n-1) - 0.25*dp) 1255 \textcolor{keywordflow}{ enddo} 1256 \textcolor{keywordflow}{if} (sv\_scale /= 1.0) \textcolor{keywordflow}{then} 1257 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks, scale=sv\_scale) 1258 \textcolor{keywordflow}{else} 1259 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks) 1260 \textcolor{keywordflow}{ endif} 1261 1262 \textcolor{comment}{! Use Boole's rule to estimate the interface height anomaly change.} 1263 intp(m) = dp*( c1\_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & 1264 12.0*a5(3))) 1265 \textcolor{keywordflow}{ enddo} 1266 \textcolor{comment}{! Use Boole's rule to integrate the interface height anomaly values in y.} 1267 inty\_dza(i,j) = c1\_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + & 1268 12.0*intp(3)) 1269 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 1270 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_a71a518f80b5f2ecf83fd9c638cad140d}\label{namespacemom__density__integrals_a71a518f80b5f2ecf83fd9c638cad140d}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+plm@{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+plm}} \index{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+plm@{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+plm}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+plm()}{int\_spec\_vol\_dp\_generic\_plm()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::int\+\_\+spec\+\_\+vol\+\_\+dp\+\_\+generic\+\_\+plm (\begin{DoxyParamCaption}\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{T\+\_\+t, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{T\+\_\+b, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{S\+\_\+t, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{S\+\_\+b, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{p\+\_\+t, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{p\+\_\+b, }\item[{real, intent(in)}]{alpha\+\_\+ref, }\item[{real, intent(in)}]{d\+P\+\_\+neglect, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(in)}]{bathyP, }\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout)}]{dza, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intp\+\_\+dza, }\item[{real, dimension( hi \%isdb\+: hi \%iedb, hi \%jsd\+: hi \%jed), intent(inout), optional}]{intx\+\_\+dza, }\item[{real, dimension( hi \%isd\+: hi \%ied, hi \%jsdb\+: hi \%jedb), intent(inout), optional}]{inty\+\_\+dza, }\item[{logical, intent(in), optional}]{use\+Mass\+Wght\+Interp }\end{DoxyParamCaption})} This subroutine calculates integrals of specific volume anomalies in pressure across layers, which are required for calculating the finite-\/volume form pressure accelerations in a non-\/\+Boussinesq model. There are essentially no free assumptions, apart from the use of Boole\textquotesingle{}s rule quadrature to do the integrals. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em hi} & A horizontal index type structure.\\ \hline \mbox{\tt in} & {\em t\+\_\+t} & Potential temperature at the top of the layer \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em t\+\_\+b} & Potential temperature at the bottom of the layer \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+t} & Salinity at the top the layer \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em s\+\_\+b} & Salinity at the bottom the layer \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+t} & Pressure atop the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+b} & Pressure below the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em alpha\+\_\+ref} & A mean specific volume that is subtracted out to reduce the magnitude of each of the integrals \mbox{[}R-\/1 $\sim$$>$ m3 kg-\/1\mbox{]} The calculation is mathematically identical with different values of alpha\+\_\+ref, but alpha\+\_\+ref alters the effects of roundoff, and answers do change.\\ \hline \mbox{\tt in} & {\em dp\+\_\+neglect} & !$<$ A miniscule pressure change with the same units as p\+\_\+t \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em bathyp} & The pressure at the bathymetry \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in,out} & {\em dza} & The change in the geopotential anomaly\\ \hline \mbox{\tt in,out} & {\em intp\+\_\+dza} & The integral in pressure through the layer of\\ \hline \mbox{\tt in,out} & {\em intx\+\_\+dza} & The integral in x of the difference between\\ \hline \mbox{\tt in,out} & {\em inty\+\_\+dza} & The integral in y of the difference between\\ \hline \mbox{\tt in} & {\em usemasswghtinterp} & If true, uses mass weighting to interpolate T/S for top and bottom integrals. \\ \hline \end{DoxyParams} Definition at line 1280 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 1280 \textcolor{keywordtype}{type}(hor\_index\_type), \textcolor{keywordtype}{intent(in)} :: hi\textcolor{comment}{ !< A horizontal index type structure.} 1281 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1282 \textcolor{keywordtype}{intent(in)} :: t\_t\textcolor{comment}{ !< Potential temperature at the top of the layer [degC]} 1283 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1284 \textcolor{keywordtype}{intent(in)} :: t\_b\textcolor{comment}{ !< Potential temperature at the bottom of the layer [degC]} 1285 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1286 \textcolor{keywordtype}{intent(in)} :: s\_t\textcolor{comment}{ !< Salinity at the top the layer [ppt]} 1287 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1288 \textcolor{keywordtype}{intent(in)} :: s\_b\textcolor{comment}{ !< Salinity at the bottom the layer [ppt]} 1289 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1290 \textcolor{keywordtype}{intent(in)} :: p\_t\textcolor{comment}{ !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]} 1291 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1292 \textcolor{keywordtype}{intent(in)} :: p\_b\textcolor{comment}{ !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]} 1293 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: alpha\_ref\textcolor{comment}{ !< A mean specific volume that is subtracted out} 1294 \textcolor{comment}{ !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]} 1295 \textcolor{comment}{ !! The calculation is mathematically identical with different values of} 1296 \textcolor{comment}{ !! alpha\_ref, but alpha\_ref alters the effects of roundoff, and} 1297 \textcolor{comment}{ !! answers do change.} 1298 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: dp\_neglect\textcolor{comment}{ ! Pa] or [Pa]} 1300 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1301 \textcolor{keywordtype}{intent(in)} :: bathyp\textcolor{comment}{ !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]} 1302 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 1303 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 1304 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1305 \textcolor{keywordtype}{intent(inout)} :: dza\textcolor{comment}{ !< The change in the geopotential anomaly} 1306 \textcolor{comment}{ !! across the layer [L2 T-2 ~> m2 s-2]} 1307 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1308 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intp\_dza\textcolor{comment}{ !< The integral in pressure through the layer of} 1309 \textcolor{comment}{ !! the geopotential anomaly relative to the anomaly at the bottom of the} 1310 \textcolor{comment}{ !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]} 1311 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(HI),SZJ\_(HI))}, & 1312 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intx\_dza\textcolor{comment}{ !< The integral in x of the difference between} 1313 \textcolor{comment}{ !! the geopotential anomaly at the top and bottom of the layer divided} 1314 \textcolor{comment}{ !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1315 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJB\_(HI))}, & 1316 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: inty\_dza\textcolor{comment}{ !< The integral in y of the difference between} 1317 \textcolor{comment}{ !! the geopotential anomaly at the top and bottom of the layer divided} 1318 \textcolor{comment}{ !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1319 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: usemasswghtinterp\textcolor{comment}{ !< If true, uses mass weighting} 1320 \textcolor{comment}{ !! to interpolate T/S for top and bottom integrals.} 1321 1322 \textcolor{comment}{! This subroutine calculates analytical and nearly-analytical integrals in} 1323 \textcolor{comment}{! pressure across layers of geopotential anomalies, which are required for} 1324 \textcolor{comment}{! calculating the finite-volume form pressure accelerations in a non-Boussinesq} 1325 \textcolor{comment}{! model. There are essentially no free assumptions, apart from the use of} 1326 \textcolor{comment}{! Boole's rule to do the horizontal integrals, and from a truncation in the} 1327 \textcolor{comment}{! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.} 1328 1329 \textcolor{keywordtype}{real} :: t5(5) \textcolor{comment}{! Temperatures at five quadrature points [degC]} 1330 \textcolor{keywordtype}{real} :: s5(5) \textcolor{comment}{! Salinities at five quadrature points [ppt]} 1331 \textcolor{keywordtype}{real} :: p5(5) \textcolor{comment}{! Pressures at five quadrature points, scaled back to Pa as necessary [Pa]} 1332 \textcolor{keywordtype}{real} :: a5(5) \textcolor{comment}{! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]} 1333 \textcolor{keywordtype}{real} :: t15(15) \textcolor{comment}{! Temperatures at fifteen interior quadrature points [degC]} 1334 \textcolor{keywordtype}{real} :: s15(15) \textcolor{comment}{! Salinities at fifteen interior quadrature points [ppt]} 1335 \textcolor{keywordtype}{real} :: p15(15) \textcolor{comment}{! Pressures at fifteen quadrature points, scaled back to Pa as necessary [Pa]} 1336 \textcolor{keywordtype}{real} :: a15(15) \textcolor{comment}{! Specific volumes at fifteen quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1]} 1337 \textcolor{keywordtype}{real} :: wt\_t(5), wt\_b(5) \textcolor{comment}{! Weights of top and bottom values at quadrature points [nondim]} 1338 \textcolor{keywordtype}{real} :: t\_top, t\_bot, s\_top, s\_bot, p\_top, p\_bot 1339 1340 \textcolor{keywordtype}{real} :: alpha\_anom \textcolor{comment}{! The depth averaged specific density anomaly [m3 kg-1]} 1341 \textcolor{keywordtype}{real} :: dp \textcolor{comment}{! The pressure change through a layer [R L2 T-2 ~> Pa]} 1342 \textcolor{keywordtype}{real} :: dp\_90(2:4) \textcolor{comment}{! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa]} 1343 \textcolor{keywordtype}{real} :: hwght \textcolor{comment}{! A pressure-thickness below topography [R L2 T-2 ~> Pa]} 1344 \textcolor{keywordtype}{real} :: hl, hr \textcolor{comment}{! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]} 1345 \textcolor{keywordtype}{real} :: alpha\_ref\_mks \textcolor{comment}{! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1]} 1346 \textcolor{keywordtype}{real} :: idenom \textcolor{comment}{! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]} 1347 \textcolor{keywordtype}{real} :: hwt\_ll, hwt\_lr \textcolor{comment}{! hWt\_LA is the weighted influence of A on the left column [nondim]} 1348 \textcolor{keywordtype}{real} :: hwt\_rl, hwt\_rr \textcolor{comment}{! hWt\_RA is the weighted influence of A on the right column [nondim]} 1349 \textcolor{keywordtype}{real} :: wt\_l, wt\_r \textcolor{comment}{! The linear weights of the left and right columns [nondim]} 1350 \textcolor{keywordtype}{real} :: wtt\_l, wtt\_r \textcolor{comment}{! The weights for tracers from the left and right columns [nondim]} 1351 \textcolor{keywordtype}{real} :: intp(5) \textcolor{comment}{! The integrals of specific volume with pressure at the} 1352 \textcolor{comment}{! 5 sub-column locations [L2 T-2 ~> m2 s-2]} 1353 \textcolor{keywordtype}{real} :: rl2\_t2\_to\_pa \textcolor{comment}{! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1]} 1354 \textcolor{keywordtype}{real} :: sv\_scale \textcolor{comment}{! A multiplicative factor by which to scale specific} 1355 \textcolor{comment}{! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]} 1356 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_90 = 1.0/90.0 \textcolor{comment}{! A rational constant.} 1357 \textcolor{keywordtype}{logical} :: do\_massweight \textcolor{comment}{! Indicates whether to do mass weighting.} 1358 \textcolor{keywordtype}{integer} :: isq, ieq, jsq, jeq, i, j, m, n, pos 1359 1360 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB 1361 1362 do\_massweight = .false. 1363 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(usemasswghtinterp)) do\_massweight = usemasswghtinterp 1364 1365 sv\_scale = us%R\_to\_kg\_m3 1366 rl2\_t2\_to\_pa = us%RL2\_T2\_to\_Pa 1367 alpha\_ref\_mks = alpha\_ref * us%kg\_m3\_to\_R 1368 1369 \textcolor{keywordflow}{do} n = 1, 5 \textcolor{comment}{! Note that these are reversed from int\_density\_dz.} 1370 wt\_t(n) = 0.25 * \textcolor{keywordtype}{real}(n-1) 1371 wt\_b(n) = 1.0 - wt\_t(n) 1372 \textcolor{keywordflow}{ enddo} 1373 1374 \textcolor{comment}{! 1. Compute vertical integrals} 1375 \textcolor{keywordflow}{do} j=jsq,jeq+1; \textcolor{keywordflow}{do} i=isq,ieq+1 1376 dp = p\_b(i,j) - p\_t(i,j) 1377 \textcolor{keywordflow}{do} n=1,5 \textcolor{comment}{! T, S and p are linearly interpolated in the vertical.} 1378 p5(n) = rl2\_t2\_to\_pa * (wt\_t(n) * p\_t(i,j) + wt\_b(n) * p\_b(i,j)) 1379 s5(n) = wt\_t(n) * s\_t(i,j) + wt\_b(n) * s\_b(i,j) 1380 t5(n) = wt\_t(n) * t\_t(i,j) + wt\_b(n) * t\_b(i,j) 1381 \textcolor{keywordflow}{ enddo} 1382 \textcolor{keywordflow}{if} (sv\_scale /= 1.0) \textcolor{keywordflow}{then} 1383 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks, scale=sv\_scale) 1384 \textcolor{keywordflow}{else} 1385 \textcolor{keyword}{call }calculate\_spec\_vol(t5, s5, p5, a5, 1, 5, eos, alpha\_ref\_mks) 1386 \textcolor{keywordflow}{ endif} 1387 1388 \textcolor{comment}{! Use Boole's rule to estimate the interface height anomaly change.} 1389 alpha\_anom = c1\_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3)) 1390 dza(i,j) = dp*alpha\_anom 1391 \textcolor{comment}{! Use a Boole's-rule-like fifth-order accurate estimate of the double integral of} 1392 \textcolor{comment}{! the interface height anomaly.} 1393 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intp\_dza)) intp\_dza(i,j) = 0.5*dp**2 * & 1394 (alpha\_anom - c1\_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) ) 1395 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} 1396 1397 \textcolor{comment}{! 2. Compute horizontal integrals in the x direction} 1398 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(intx\_dza)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=hi%jsc,hi%jec ; \textcolor{keywordflow}{do} i=isq,ieq 1399 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 1400 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation} 1401 \textcolor{comment}{! of T,S along the top and bottom integrals, almost like thickness} 1402 \textcolor{comment}{! weighting. Note: To work in terrain following coordinates we could} 1403 \textcolor{comment}{! offset this distance by the layer thickness to replicate other models.} 1404 hwght = 0.0 1405 \textcolor{keywordflow}{if} (do\_massweight) & 1406 hwght = max(0., bathyp(i,j)-p\_t(i+1,j), bathyp(i+1,j)-p\_t(i,j)) 1407 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 1408 hl = (p\_b(i,j) - p\_t(i,j)) + dp\_neglect 1409 hr = (p\_b(i+1,j) - p\_t(i+1,j)) + dp\_neglect 1410 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 1411 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr ) 1412 hwt\_ll = (hwght*hl + hr*hl) * idenom ; hwt\_lr = (hwght*hr) * idenom 1413 hwt\_rr = (hwght*hr + hr*hl) * idenom ; hwt\_rl = (hwght*hl) * idenom 1414 \textcolor{keywordflow}{else} 1415 hwt\_ll = 1.0 ; hwt\_lr = 0.0 ; hwt\_rr = 1.0 ; hwt\_rl = 0.0 1416 \textcolor{keywordflow}{ endif} 1417 1418 \textcolor{keywordflow}{do} m=2,4 1419 wt\_l = 0.25*\textcolor{keywordtype}{real(5-m)} ; wt\_r = 1.0-wt\_l 1420 wtt\_l = wt\_l*hwt\_ll + wt\_r*hwt\_rl ; wtt\_r = wt\_l*hwt\_lr + wt\_r*hwt\_rr 1421 1422 \textcolor{comment}{! T, S, and p are interpolated in the horizontal. The p interpolation} 1423 \textcolor{comment}{! is linear, but for T and S it may be thickness weighted.} 1424 p\_top = wt\_l*p\_t(i,j) + wt\_r*p\_t(i+1,j) 1425 p\_bot = wt\_l*p\_b(i,j) + wt\_r*p\_b(i+1,j) 1426 t\_top = wtt\_l*t\_t(i,j) + wtt\_r*t\_t(i+1,j) 1427 t\_bot = wtt\_l*t\_b(i,j) + wtt\_r*t\_b(i+1,j) 1428 s\_top = wtt\_l*s\_t(i,j) + wtt\_r*s\_t(i+1,j) 1429 s\_bot = wtt\_l*s\_b(i,j) + wtt\_r*s\_b(i+1,j) 1430 dp\_90(m) = c1\_90*(p\_bot - p\_top) 1431 1432 \textcolor{comment}{! Salinity, temperature and pressure with linear interpolation in the vertical.} 1433 pos = (m-2)*5 1434 \textcolor{keywordflow}{do} n=1,5 1435 p15(pos+n) = rl2\_t2\_to\_pa * (wt\_t(n) * p\_top + wt\_b(n) * p\_bot) 1436 s15(pos+n) = wt\_t(n) * s\_top + wt\_b(n) * s\_bot 1437 t15(pos+n) = wt\_t(n) * t\_top + wt\_b(n) * t\_bot 1438 \textcolor{keywordflow}{ enddo} 1439 \textcolor{keywordflow}{ enddo} 1440 1441 \textcolor{keywordflow}{if} (sv\_scale /= 1.0) \textcolor{keywordflow}{then} 1442 \textcolor{keyword}{call }calculate\_spec\_vol(t15, s15, p15, a15, 1, 15, eos, alpha\_ref\_mks, scale=sv\_scale) 1443 \textcolor{keywordflow}{else} 1444 \textcolor{keyword}{call }calculate\_spec\_vol(t15, s15, p15, a15, 1, 15, eos, alpha\_ref\_mks) 1445 \textcolor{keywordflow}{ endif} 1446 1447 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) 1448 \textcolor{keywordflow}{do} m=2,4 1449 \textcolor{comment}{! Use Boole's rule to estimate the interface height anomaly change.} 1450 \textcolor{comment}{! The integrals at the ends of the segment are already known.} 1451 pos = (m-2)*5 1452 intp(m) = dp\_90(m)*((7.0*(a15(pos+1)+a15(pos+5)) + & 1453 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3)) 1454 \textcolor{keywordflow}{ enddo} 1455 \textcolor{comment}{! Use Boole's rule to integrate the interface height anomaly values in x.} 1456 intx\_dza(i,j) = c1\_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + & 1457 12.0*intp(3)) 1458 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 1459 1460 \textcolor{comment}{! 3. Compute horizontal integrals in the y direction} 1461 \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(inty\_dza)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} i=hi%isc,hi%iec 1462 \textcolor{comment}{! hWght is the distance measure by which the cell is violation of} 1463 \textcolor{comment}{! hydrostatic consistency. For large hWght we bias the interpolation} 1464 \textcolor{comment}{! of T,S along the top and bottom integrals, like thickness weighting.} 1465 hwght = 0.0 1466 \textcolor{keywordflow}{if} (do\_massweight) & 1467 hwght = max(0., bathyp(i,j)-p\_t(i,j+1), bathyp(i,j+1)-p\_t(i,j)) 1468 \textcolor{keywordflow}{if} (hwght > 0.) \textcolor{keywordflow}{then} 1469 hl = (p\_b(i,j) - p\_t(i,j)) + dp\_neglect 1470 hr = (p\_b(i,j+1) - p\_t(i,j+1)) + dp\_neglect 1471 hwght = hwght * ( (hl-hr)/(hl+hr) )**2 1472 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr ) 1473 hwt\_ll = (hwght*hl + hr*hl) * idenom ; hwt\_lr = (hwght*hr) * idenom 1474 hwt\_rr = (hwght*hr + hr*hl) * idenom ; hwt\_rl = (hwght*hl) * idenom 1475 \textcolor{keywordflow}{else} 1476 hwt\_ll = 1.0 ; hwt\_lr = 0.0 ; hwt\_rr = 1.0 ; hwt\_rl = 0.0 1477 \textcolor{keywordflow}{ endif} 1478 1479 \textcolor{keywordflow}{do} m=2,4 1480 wt\_l = 0.25*\textcolor{keywordtype}{real(5-m)} ; wt\_r = 1.0-wt\_l 1481 wtt\_l = wt\_l*hwt\_ll + wt\_r*hwt\_rl ; wtt\_r = wt\_l*hwt\_lr + wt\_r*hwt\_rr 1482 1483 \textcolor{comment}{! T, S, and p are interpolated in the horizontal. The p interpolation} 1484 \textcolor{comment}{! is linear, but for T and S it may be thickness weighted.} 1485 p\_top = wt\_l*p\_t(i,j) + wt\_r*p\_t(i,j+1) 1486 p\_bot = wt\_l*p\_b(i,j) + wt\_r*p\_b(i,j+1) 1487 t\_top = wtt\_l*t\_t(i,j) + wtt\_r*t\_t(i,j+1) 1488 t\_bot = wtt\_l*t\_b(i,j) + wtt\_r*t\_b(i,j+1) 1489 s\_top = wtt\_l*s\_t(i,j) + wtt\_r*s\_t(i,j+1) 1490 s\_bot = wtt\_l*s\_b(i,j) + wtt\_r*s\_b(i,j+1) 1491 dp\_90(m) = c1\_90*(p\_bot - p\_top) 1492 1493 \textcolor{comment}{! Salinity, temperature and pressure with linear interpolation in the vertical.} 1494 pos = (m-2)*5 1495 \textcolor{keywordflow}{do} n=1,5 1496 p15(pos+n) = rl2\_t2\_to\_pa * (wt\_t(n) * p\_top + wt\_b(n) * p\_bot) 1497 s15(pos+n) = wt\_t(n) * s\_top + wt\_b(n) * s\_bot 1498 t15(pos+n) = wt\_t(n) * t\_top + wt\_b(n) * t\_bot 1499 \textcolor{keywordflow}{ enddo} 1500 \textcolor{keywordflow}{ enddo} 1501 1502 \textcolor{keywordflow}{if} (sv\_scale /= 1.0) \textcolor{keywordflow}{then} 1503 \textcolor{keyword}{call }calculate\_spec\_vol(t15, s15, p15, a15, 1, 15, eos, alpha\_ref\_mks, scale=sv\_scale) 1504 \textcolor{keywordflow}{else} 1505 \textcolor{keyword}{call }calculate\_spec\_vol(t15, s15, p15, a15, 1, 15, eos, alpha\_ref\_mks) 1506 \textcolor{keywordflow}{ endif} 1507 1508 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) 1509 \textcolor{keywordflow}{do} m=2,4 1510 \textcolor{comment}{! Use Boole's rule to estimate the interface height anomaly change.} 1511 \textcolor{comment}{! The integrals at the ends of the segment are already known.} 1512 pos = (m-2)*5 1513 intp(m) = dp\_90(m) * ((7.0*(a15(pos+1)+a15(pos+5)) + & 1514 32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3)) 1515 \textcolor{keywordflow}{ enddo} 1516 \textcolor{comment}{! Use Boole's rule to integrate the interface height anomaly values in x.} 1517 inty\_dza(i,j) = c1\_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + & 1518 12.0*intp(3)) 1519 \textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif} 1520 \end{DoxyCode} \mbox{\Hypertarget{namespacemom__density__integrals_a759c2ae7aec17c59d532050f68a1e518}\label{namespacemom__density__integrals_a759c2ae7aec17c59d532050f68a1e518}} \index{mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}!int\+\_\+specific\+\_\+vol\+\_\+dp@{int\+\_\+specific\+\_\+vol\+\_\+dp}} \index{int\+\_\+specific\+\_\+vol\+\_\+dp@{int\+\_\+specific\+\_\+vol\+\_\+dp}!mom\+\_\+density\+\_\+integrals@{mom\+\_\+density\+\_\+integrals}} \subsubsection{\texorpdfstring{int\+\_\+specific\+\_\+vol\+\_\+dp()}{int\_specific\_vol\_dp()}} {\footnotesize\ttfamily subroutine, public mom\+\_\+density\+\_\+integrals\+::int\+\_\+specific\+\_\+vol\+\_\+dp (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{T, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{S, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{p\+\_\+t, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in)}]{p\+\_\+b, }\item[{real, intent(in)}]{alpha\+\_\+ref, }\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI, }\item[{type(eos\+\_\+type), pointer}]{E\+OS, }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(inout)}]{dza, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(inout), optional}]{intp\+\_\+dza, }\item[{real, dimension(szib\+\_\+(hi),szj\+\_\+(hi)), intent(inout), optional}]{intx\+\_\+dza, }\item[{real, dimension(szi\+\_\+(hi),szjb\+\_\+(hi)), intent(inout), optional}]{inty\+\_\+dza, }\item[{integer, intent(in), optional}]{halo\+\_\+size, }\item[{real, dimension(szi\+\_\+(hi),szj\+\_\+(hi)), intent(in), optional}]{bathyP, }\item[{real, intent(in), optional}]{d\+P\+\_\+tiny, }\item[{logical, intent(in), optional}]{use\+Mass\+Wght\+Interp }\end{DoxyParamCaption})} Calls the appropriate subroutine to calculate analytical and nearly-\/analytical integrals in pressure across layers of geopotential anomalies, which are required for calculating the finite-\/volume form pressure accelerations in a non-\/\+Boussinesq model. There are essentially no free assumptions, apart from the use of Boole\textquotesingle{}s rule to do the horizontal integrals, and from a truncation in the series for log(1-\/eps/1+eps) that assumes that $\vert$eps$\vert$ $<$ 0.\+34. \begin{DoxyParams}[1]{Parameters} \mbox{\tt in} & {\em hi} & The horizontal index structure\\ \hline \mbox{\tt in} & {\em t} & Potential temperature referenced to the surface \mbox{[}degC\mbox{]}\\ \hline \mbox{\tt in} & {\em s} & Salinity \mbox{[}ppt\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+t} & Pressure at the top of the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em p\+\_\+b} & Pressure at the bottom of the layer \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em alpha\+\_\+ref} & A mean specific volume that is subtracted out to reduce the magnitude of each of the integrals \mbox{[}R-\/1 $\sim$$>$ m3 kg-\/1\mbox{]} The calculation is mathematically identical with different values of alpha\+\_\+ref, but this reduces the effects of roundoff.\\ \hline & {\em eos} & Equation of state structure\\ \hline \mbox{\tt in} & {\em us} & A dimensional unit scaling type\\ \hline \mbox{\tt in,out} & {\em dza} & The change in the geopotential anomaly across\\ \hline \mbox{\tt in,out} & {\em intp\+\_\+dza} & The integral in pressure through the layer of the\\ \hline \mbox{\tt in,out} & {\em intx\+\_\+dza} & The integral in x of the difference between the\\ \hline \mbox{\tt in,out} & {\em inty\+\_\+dza} & The integral in y of the difference between the\\ \hline \mbox{\tt in} & {\em halo\+\_\+size} & The width of halo points on which to calculate dza.\\ \hline \mbox{\tt in} & {\em bathyp} & The pressure at the bathymetry \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em dp\+\_\+tiny} & A minuscule pressure change with the same units as p\+\_\+t \mbox{[}R L2 T-\/2 $\sim$$>$ Pa\mbox{]} or \mbox{[}Pa\mbox{]}\\ \hline \mbox{\tt in} & {\em usemasswghtinterp} & If true, uses mass weighting to interpolate T/S for top and bottom integrals. \\ \hline \end{DoxyParams} Definition at line 1007 of file M\+O\+M\+\_\+density\+\_\+integrals.\+F90. \begin{DoxyCode} 1007 \textcolor{keywordtype}{type}(hor\_index\_type), \textcolor{keywordtype}{intent(in)} :: hi\textcolor{comment}{ !< The horizontal index structure} 1008 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1009 \textcolor{keywordtype}{intent(in)} :: t\textcolor{comment}{ !< Potential temperature referenced to the surface [degC]} 1010 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1011 \textcolor{keywordtype}{intent(in)} :: s\textcolor{comment}{ !< Salinity [ppt]} 1012 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1013 \textcolor{keywordtype}{intent(in)} :: p\_t\textcolor{comment}{ !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]} 1014 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1015 \textcolor{keywordtype}{intent(in)} :: p\_b\textcolor{comment}{ !< Pressure at the bottom of the layer [R L2 T-2 ~> Pa] or [Pa]} 1016 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{intent(in)} :: alpha\_ref\textcolor{comment}{ !< A mean specific volume that is subtracted out} 1017 \textcolor{comment}{ !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]} 1018 \textcolor{comment}{ !! The calculation is mathematically identical with different values of} 1019 \textcolor{comment}{ !! alpha\_ref, but this reduces the effects of roundoff.} 1020 \textcolor{keywordtype}{type}(eos\_type), \textcolor{keywordtype}{pointer} :: eos\textcolor{comment}{ !< Equation of state structure} 1021 \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{ !< A dimensional unit scaling type} 1022 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1023 \textcolor{keywordtype}{intent(inout)} :: dza\textcolor{comment}{ !< The change in the geopotential anomaly across} 1024 \textcolor{comment}{ !! the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1025 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1026 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intp\_dza\textcolor{comment}{ !< The integral in pressure through the layer of the} 1027 \textcolor{comment}{ !! geopotential anomaly relative to the anomaly at the bottom of the} 1028 \textcolor{comment}{ !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]} 1029 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(HI),SZJ\_(HI))}, & 1030 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: intx\_dza\textcolor{comment}{ !< The integral in x of the difference between the} 1031 \textcolor{comment}{ !! geopotential anomaly at the top and bottom of the layer divided by} 1032 \textcolor{comment}{ !! the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1033 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJB\_(HI))}, & 1034 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(inout)} :: inty\_dza\textcolor{comment}{ !< The integral in y of the difference between the} 1035 \textcolor{comment}{ !! geopotential anomaly at the top and bottom of the layer divided by} 1036 \textcolor{comment}{ !! the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]} 1037 \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: halo\_size\textcolor{comment}{ !< The width of halo points on which to calculate dza.} 1038 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(HI),SZJ\_(HI))}, & 1039 \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: bathyp\textcolor{comment}{ !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]} 1040 \textcolor{keywordtype}{real}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: dp\_tiny\textcolor{comment}{ !< A minuscule pressure change with} 1041 \textcolor{comment}{ !! the same units as p\_t [R L2 T-2 ~> Pa] or [Pa]} 1042 \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)} :: usemasswghtinterp\textcolor{comment}{ !< If true, uses mass weighting} 1043 \textcolor{comment}{ !! to interpolate T/S for top and bottom integrals.} 1044 1045 \textcolor{keywordflow}{if} (eos\_quadrature(eos)) \textcolor{keywordflow}{then} 1046 \textcolor{keyword}{call }int\_spec\_vol\_dp\_generic\_pcm(t, s, p\_t, p\_b, alpha\_ref, hi, eos, us, & 1047 dza, intp\_dza, intx\_dza, inty\_dza, halo\_size, & 1048 bathyp, dp\_tiny, usemasswghtinterp) 1049 \textcolor{keywordflow}{else} 1050 \textcolor{keyword}{call }analytic\_int\_specific\_vol\_dp(t, s, p\_t, p\_b, alpha\_ref, hi, eos, & 1051 dza, intp\_dza, intx\_dza, inty\_dza, halo\_size, & 1052 bathyp, dp\_tiny, usemasswghtinterp) 1053 \textcolor{keywordflow}{ endif} 1054 \end{DoxyCode}