\hypertarget{namespacemom__thickness__diffuse}{}\section{mom\+\_\+thickness\+\_\+diffuse Module Reference}
\label{namespacemom__thickness__diffuse}\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}


\subsection{Detailed Description}
Thickness diffusion (or Gent Mc\+Williams) 

\hypertarget{namespacemom__thickness__diffuse_section_gm}{}\subsection{Thickness diffusion (aka Gent-\/\+Mc\+Williams)}\label{namespacemom__thickness__diffuse_section_gm}
Thickness diffusion is implemented via along-\/layer mass fluxes \[ h^\dagger \leftarrow h^n - \Delta t \nabla \cdot ( \vec{uh}^* ) \] where the mass fluxes are cast as the difference in vector streamfunction

\[ \vec{uh}^* = \delta_k \vec{\psi} . \]

The GM implementation of thickness diffusion made the streamfunction proportional to the potential density slope \[ \vec{\psi} = - \kappa_h \frac{\nabla_z \rho}{\partial_z \rho} = \frac{g\kappa_h}{\rho_o} \frac{\nabla \rho}{N^2} = \kappa_h \frac{M^2}{N^2} \] but for robustness the scheme is implemented as \[ \vec{\psi} = \kappa_h \frac{M^2}{\sqrt{N^4 + M^4}} \] since the quantity $\frac{M^2}{\sqrt{N^2 + M^2}}$ is bounded between \$-\/1\$ and \$1\$ and does not change sign if $N^2<0$.

Optionally, the method of Ferrari et al, 2010, can be used to obtain the streamfunction which solves the vertically elliptic equation\+: \[ \gamma_F \partial_z c^2 \partial_z \psi - N_*^2 \psi = ( 1 + \gamma_F ) \kappa_h N_*^2 \frac{M^2}{\sqrt{N^4+M^4}} \] which recovers the previous streamfunction relation in the limit that $ c \rightarrow 0 $. Here, $c=\max(c_{min},c_g)$ is the maximum of either $c_{min}$ and either the first baroclinic mode wave-\/speed or the equivalent barotropic mode wave-\/speed. $N_*^2 = \max(N^2,0)$ is a non-\/negative form of the square of the Brunt-\/\+Vaisala frequency. The parameter $\gamma_F$ is used to reduce the vertical smoothing length scale. \[ \kappa_h = \left( \kappa_o + \alpha_{s} L_{s}^2 < S N > + \alpha_{M} \kappa_{M} \right) r(\Delta x,L_d) \] where $ S $ is the isoneutral slope magnitude, $ N $ is the square root of Brunt-\/\+Vaisala frequency, $\kappa_{M}$ is the diffusivity calculated by the M\+E\+KE parameterization (\hyperlink{namespacemom__meke}{mom\+\_\+meke} module) and $ r(\Delta x,L_d) $ is a function of the local resolution (ratio of grid-\/spacing, $\Delta x$, to deformation radius, $L_d$). The length $L_s$ is provided by the \hyperlink{namespacemom__lateral__mixing__coeffs}{mom\+\_\+lateral\+\_\+mixing\+\_\+coeffs} module (enabled with {\ttfamily U\+S\+E\+\_\+\+V\+A\+R\+I\+A\+B\+L\+E\+\_\+\+M\+I\+X\+I\+NG=True} and the term $<SN>$ is the vertical average slope times the Brunt-\/\+Vaisala frequency prescribed by Visbeck et al., 1996.

The result of the above expression is subsequently bounded by minimum and maximum values, including an upper diffusivity consistent with numerical stability ( $ \kappa_{cfl} $ is calculated internally). \[ \kappa_h \leftarrow \min{\left( \kappa_{max}, \kappa_{cfl}, \max{\left( \kappa_{min}, \kappa_h \right)} \right)} f(c_g,z) \]

where $f(c_g,z)$ is a vertical structure function. $f(c_g,z)$ is calculated in module \hyperlink{namespacemom__lateral__mixing__coeffs}{mom\+\_\+lateral\+\_\+mixing\+\_\+coeffs}. If {\ttfamily K\+H\+T\+H\+\_\+\+U\+S\+E\+\_\+\+E\+B\+T\+\_\+\+S\+T\+R\+U\+CT=True} then $f(c_g,z)$ is set to look like the equivalent barotropic modal velocity structure. Otherwise $f(c_g,z)=1$ and the diffusivity is independent of depth.

In order to calculate meaningful slopes in vanished layers, temporary copies of the thermodynamic variables are passed through a vertical smoother, function vert\+\_\+fill\+\_\+ts()\+: \begin{eqnarray*} \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] \theta & \leftarrow & \theta \\ \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] s & \leftarrow & s \end{eqnarray*}\hypertarget{namespacemom__thickness__diffuse_section_khth_module_parameters}{}\subsubsection{Module mom\+\_\+thickness\+\_\+diffuse parameters}\label{namespacemom__thickness__diffuse_section_khth_module_parameters}
\tabulinesep=1mm
\begin{longtabu} spread 0pt [c]{*{2}{|X[-1]}|}
\hline
\rowcolor{\tableheadbgcolor}\textbf{ Symbol }&\textbf{ Module parameter  }\\\cline{1-2}
\endfirsthead
\hline
\endfoot
\hline
\rowcolor{\tableheadbgcolor}\textbf{ Symbol }&\textbf{ Module parameter  }\\\cline{1-2}
\endhead
-\/ &{\ttfamily T\+H\+I\+C\+K\+N\+E\+S\+S\+D\+I\+F\+F\+U\+SE} \\\cline{1-2}
$ \kappa_o $ &{\ttfamily K\+H\+TH} \\\cline{1-2}
$ \alpha_{s} $ &{\ttfamily K\+H\+T\+H\+\_\+\+S\+L\+O\+P\+E\+\_\+\+C\+FF} \\\cline{1-2}
$ \kappa_{min} $ &{\ttfamily K\+H\+T\+H\+\_\+\+M\+IN} \\\cline{1-2}
$ \kappa_{max} $ &{\ttfamily K\+H\+T\+H\+\_\+\+M\+AX} \\\cline{1-2}
-\/ &{\ttfamily K\+H\+T\+H\+\_\+\+M\+A\+X\+\_\+\+C\+FL} \\\cline{1-2}
$ \kappa_{smth} $ &{\ttfamily K\+D\+\_\+\+S\+M\+O\+O\+TH} \\\cline{1-2}
$ \alpha_{M} $ &{\ttfamily M\+E\+K\+E\+\_\+\+K\+H\+T\+H\+\_\+\+F\+AC} (from \hyperlink{namespacemom__meke}{mom\+\_\+meke} module) \\\cline{1-2}
-\/ &{\ttfamily K\+H\+T\+H\+\_\+\+U\+S\+E\+\_\+\+E\+B\+T\+\_\+\+S\+T\+R\+U\+CT} (from \hyperlink{namespacemom__lateral__mixing__coeffs}{mom\+\_\+lateral\+\_\+mixing\+\_\+coeffs} module) \\\cline{1-2}
-\/ &{\ttfamily K\+H\+T\+H\+\_\+\+U\+S\+E\+\_\+\+F\+G\+N\+V\+\_\+\+S\+T\+R\+E\+A\+M\+F\+U\+N\+C\+T\+I\+ON} \\\cline{1-2}
$ \gamma_F $ &{\ttfamily F\+G\+N\+V\+\_\+\+F\+I\+L\+T\+E\+R\+\_\+\+S\+C\+A\+LE} \\\cline{1-2}
$ c_{min} $ &{\ttfamily F\+G\+N\+V\+\_\+\+C\+\_\+\+M\+IN} \\\cline{1-2}
\end{longtabu}
\hypertarget{namespacemom__thickness__diffuse_section_khth_module_reference}{}\subsubsection{References}\label{namespacemom__thickness__diffuse_section_khth_module_reference}
Ferrari, R., S.\+M. Griffies, A.\+J.\+G. Nurser and G.\+K. Vallis, 2010\+: A boundary-\/value problem for the parameterized mesoscale eddy transport. Ocean Modelling, 32, 143-\/156. \href{http://doi.org/10.1016/j.ocemod.2010.01.004}{\tt http\+://doi.\+org/10.\+1016/j.\+ocemod.\+2010.\+01.\+004}

Viscbeck, M., J.\+C. Marshall, H. Jones, 1996\+: Dynamics of isolated convective regions in the ocean. J. Phys. Oceangr., 26, 1721-\/1734. \href{http://dx.doi.org/10.1175/1520-0485(1996)026%3C1721:DOICRI%3E2.0.CO;2}{\tt http\+://dx.\+doi.\+org/10.\+1175/1520-\/0485(1996)026\%3\+C1721\+:\+D\+O\+I\+C\+R\+I\%3\+E2.\+0.\+C\+O;2} \subsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \hyperlink{structmom__thickness__diffuse_1_1thickness__diffuse__cs}{thickness\+\_\+diffuse\+\_\+cs}
\begin{DoxyCompactList}\small\item\em Control structure for thickness diffusion. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \hyperlink{namespacemom__thickness__diffuse_a8a538b778a567f489bfd9c5eadeeebef}{thickness\+\_\+diffuse} (h, uhtr, vhtr, tv, dt, G, GV, US, M\+E\+KE, Var\+Mix, C\+Dp, CS)
\begin{DoxyCompactList}\small\item\em Calculates thickness diffusion coefficients and applies thickness diffusion to layer thicknesses, h. Diffusivities are limited to ensure stability. Also returns along-\/layer mass fluxes used in the continuity equation. \end{DoxyCompactList}\item 
subroutine \hyperlink{namespacemom__thickness__diffuse_ae9909642254fcf0160afe46997e10c30}{thickness\+\_\+diffuse\+\_\+full} (h, e, Kh\+\_\+u, Kh\+\_\+v, tv, uhD, vhD, cg1, dt, G, GV, US, M\+E\+KE, CS, int\+\_\+slope\+\_\+u, int\+\_\+slope\+\_\+v, slope\+\_\+x, slope\+\_\+y)
\begin{DoxyCompactList}\small\item\em Calculates parameterized layer transports for use in the continuity equation. Fluxes are limited to give positive definite thicknesses. Called by \hyperlink{namespacemom__thickness__diffuse_a8a538b778a567f489bfd9c5eadeeebef}{thickness\+\_\+diffuse()}. \end{DoxyCompactList}\item 
subroutine \hyperlink{namespacemom__thickness__diffuse_a52d5fe57d53414fdc05f669723c9774e}{streamfn\+\_\+solver} (nk, c2\+\_\+h, h\+N2, sfn)
\begin{DoxyCompactList}\small\item\em Tridiagonal solver for streamfunction at interfaces. \end{DoxyCompactList}\item 
subroutine \hyperlink{namespacemom__thickness__diffuse_ab6206370a3f642ad57c63b6e268ee0fb}{add\+\_\+detangling\+\_\+kh} (h, e, Kh\+\_\+u, Kh\+\_\+v, K\+H\+\_\+u\+\_\+\+C\+FL, K\+H\+\_\+v\+\_\+\+C\+FL, tv, dt, G, GV, US, CS, int\+\_\+slope\+\_\+u, int\+\_\+slope\+\_\+v)
\begin{DoxyCompactList}\small\item\em Modifies thickness diffusivities to untangle layer structures. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__thickness__diffuse_affd209ab72d54799680c78d8ba1ad779}{thickness\+\_\+diffuse\+\_\+init} (Time, G, GV, US, param\+\_\+file, diag, C\+Dp, CS)
\begin{DoxyCompactList}\small\item\em Initialize the thickness diffusion module/structure. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__thickness__diffuse_a96e66e8491141614027c1583aeaecd7a}{thickness\+\_\+diffuse\+\_\+get\+\_\+kh} (CS, K\+H\+\_\+u\+\_\+\+G\+ME, K\+H\+\_\+v\+\_\+\+G\+ME, G)
\begin{DoxyCompactList}\small\item\em Copies ubtav and vbtav from private type into arrays. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__thickness__diffuse_a1a3a5c750e8a2323afa82b118329378c}{thickness\+\_\+diffuse\+\_\+end} (CS)
\begin{DoxyCompactList}\small\item\em Deallocate the thickness diffusion control structure. \end{DoxyCompactList}\end{DoxyCompactItemize}


\subsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespacemom__thickness__diffuse_ab6206370a3f642ad57c63b6e268ee0fb}\label{namespacemom__thickness__diffuse_ab6206370a3f642ad57c63b6e268ee0fb}} 
\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}!add\+\_\+detangling\+\_\+kh@{add\+\_\+detangling\+\_\+kh}}
\index{add\+\_\+detangling\+\_\+kh@{add\+\_\+detangling\+\_\+kh}!mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}
\subsubsection{\texorpdfstring{add\+\_\+detangling\+\_\+kh()}{add\_detangling\_kh()}}
{\footnotesize\ttfamily subroutine mom\+\_\+thickness\+\_\+diffuse\+::add\+\_\+detangling\+\_\+kh (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(in)}]{e,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(inout)}]{Kh\+\_\+u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)+1), intent(inout)}]{Kh\+\_\+v,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g)), intent(in)}]{K\+H\+\_\+u\+\_\+\+C\+FL,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g)), intent(in)}]{K\+H\+\_\+v\+\_\+\+C\+FL,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\hyperlink{structmom__thickness__diffuse_1_1thickness__diffuse__cs}{thickness\+\_\+diffuse\+\_\+cs}), pointer}]{CS,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(inout)}]{int\+\_\+slope\+\_\+u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)+1), intent(inout)}]{int\+\_\+slope\+\_\+v }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Modifies thickness diffusivities to untangle layer structures. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em e} & Interface positions \mbox{[}Z $\sim$$>$ m\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em kh\+\_\+u} & Thickness diffusivity on interfaces at u points \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em kh\+\_\+v} & Thickness diffusivity on interfaces at v points \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em kh\+\_\+u\+\_\+cfl} & Maximum stable thickness diffusivity at u points \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em kh\+\_\+v\+\_\+cfl} & Maximum stable thickness diffusivity at v points \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & Thermodynamics structure\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
 & {\em cs} & Control structure for thickness diffusion\\
\hline
\mbox{\tt in,out}  & {\em int\+\_\+slope\+\_\+u} & Ratio that determine how much of the isopycnal slopes are taken directly from the interface slopes without consideration of density gradients.\\
\hline
\mbox{\tt in,out}  & {\em int\+\_\+slope\+\_\+v} & Ratio that determine how much of the isopycnal slopes are taken directly from the interface slopes without consideration of density gradients. \\
\hline
\end{DoxyParams}


Definition at line 1468 of file M\+O\+M\+\_\+thickness\+\_\+diffuse.\+F90.


\begin{DoxyCode}
1468   \textcolor{keywordtype}{type}(ocean\_grid\_type),                       \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{    !< Ocean grid structure}
1469   \textcolor{keywordtype}{type}(verticalgrid\_type),                     \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< Vertical grid structure}
1470   \textcolor{keywordtype}{type}(unit\_scale\_type),                       \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{   !< A dimensional unit scaling type}
1471   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))},    \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{    !< Layer thickness [H ~> m or kg m-2]}
1472   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G)+1)},  \textcolor{keywordtype}{intent(in)}    :: e\textcolor{comment}{    !< Interface positions [Z ~> m]}
1473   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(inout)} :: kh\_u\textcolor{comment}{ !< Thickness diffusivity on interfaces}
1474 \textcolor{comment}{                                                                     !! at u points [L2 T-1 ~> m2 s-1]}
1475   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(inout)} :: kh\_v\textcolor{comment}{ !< Thickness diffusivity on interfaces}
1476 \textcolor{comment}{                                                                     !! at v points [L2 T-1 ~> m2 s-1]}
1477   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G))},           \textcolor{keywordtype}{intent(in)}    :: kh\_u\_cfl\textcolor{comment}{ !< Maximum stable thickness
       diffusivity}
1478 \textcolor{comment}{                                                                     !! at u points [L2 T-1 ~> m2 s-1]}
1479   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G))},           \textcolor{keywordtype}{intent(in)}    :: kh\_v\_cfl\textcolor{comment}{ !< Maximum stable thickness
       diffusivity}
1480 \textcolor{comment}{                                                                     !! at v points [L2 T-1 ~> m2 s-1]}
1481   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),                       \textcolor{keywordtype}{intent(in)}    :: tv\textcolor{comment}{   !< Thermodynamics structure}
1482   \textcolor{keywordtype}{real},                                        \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{   !< Time increment [T ~> s]}
1483   \textcolor{keywordtype}{type}(thickness\_diffuse\_cs),                  \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{   !< Control structure for thickness
       diffusion}
1484   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(inout)} :: int\_slope\_u\textcolor{comment}{ !< Ratio that determine how
       much of}
1485 \textcolor{comment}{                                                                     !! the isopycnal slopes are taken
       directly from}
1486 \textcolor{comment}{                                                                     !! the interface slopes without
       consideration}
1487 \textcolor{comment}{                                                                     !! of density gradients.}
1488   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(inout)} :: int\_slope\_v\textcolor{comment}{ !< Ratio that determine how
       much of}
1489 \textcolor{comment}{                                                                     !! the isopycnal slopes are taken
       directly from}
1490 \textcolor{comment}{                                                                     !! the interface slopes without
       consideration}
1491 \textcolor{comment}{                                                                     !! of density gradients.}
1492   \textcolor{comment}{! Local variables}
1493   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))} :: &
1494     de\_top     \textcolor{comment}{! The distances between the top of a layer and the top of the}
1495                \textcolor{comment}{! region where the detangling is applied [H ~> m or kg m-2].}
1496   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))} :: &
1497     kh\_lay\_u   \textcolor{comment}{! The tentative interface height diffusivity for each layer at}
1498                \textcolor{comment}{! u points [L2 T-1 ~> m2 s-1].}
1499   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))} :: &
1500     kh\_lay\_v   \textcolor{comment}{! The tentative interface height diffusivity for each layer at}
1501                \textcolor{comment}{! v points [L2 T-1 ~> m2 s-1].}
1502   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G))} :: &
1503     de\_bot     \textcolor{comment}{! The distances from the bottom of the region where the}
1504                \textcolor{comment}{! detangling is applied [H ~> m or kg m-2].}
1505   \textcolor{keywordtype}{real} :: h1, h2    \textcolor{comment}{! The thinner and thicker surrounding thicknesses [H ~> m or kg m-2],}
1506                     \textcolor{comment}{! with the thinner modified near the boundaries to mask out}
1507                     \textcolor{comment}{! thickness variations due to topography, etc.}
1508   \textcolor{keywordtype}{real} :: jag\_rat   \textcolor{comment}{! The nondimensional jaggedness ratio for a layer, going}
1509                     \textcolor{comment}{! from 0 (smooth) to 1 (jagged).  This is the difference}
1510                     \textcolor{comment}{! between the arithmetic and harmonic mean thicknesses}
1511                     \textcolor{comment}{! normalized by the arithmetic mean thickness.}
1512   \textcolor{keywordtype}{real} :: kh\_scale  \textcolor{comment}{! A ratio by which Kh\_u\_CFL is scaled for maximally jagged}
1513                     \textcolor{comment}{! layers [nondim].}
1514   \textcolor{keywordtype}{real} :: h\_neglect \textcolor{comment}{! A thickness that is so small it is usually lost}
1515                     \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
1516 
1517   \textcolor{keywordtype}{real} :: i\_sl      \textcolor{comment}{! The absolute value of the larger in magnitude of the slopes}
1518                     \textcolor{comment}{! above and below [L Z-1 ~> nondim].}
1519   \textcolor{keywordtype}{real} :: rsl       \textcolor{comment}{! The ratio of the smaller magnitude slope to the larger}
1520                     \textcolor{comment}{! magnitude one [nondim]. 0 <= Rsl <1.}
1521   \textcolor{keywordtype}{real} :: irsl      \textcolor{comment}{! The (limited) inverse of Rsl [nondim]. 1 < IRsl <= 1e9.}
1522   \textcolor{keywordtype}{real} :: dh        \textcolor{comment}{! The thickness gradient divided by the damping timescale}
1523                     \textcolor{comment}{! and the ratio of the face length to the adjacent cell}
1524                     \textcolor{comment}{! areas for comparability with the diffusivities [L Z T-1 ~> m2 s-1].}
1525   \textcolor{keywordtype}{real} :: adh       \textcolor{comment}{! The absolute value of dH [L Z T-1 ~> m2 s-1].}
1526   \textcolor{keywordtype}{real} :: sign      \textcolor{comment}{! 1 or -1, with the same sign as the layer thickness gradient.}
1527   \textcolor{keywordtype}{real} :: sl\_k      \textcolor{comment}{! The sign-corrected slope of the interface above [Z L-1 ~> nondim].}
1528   \textcolor{keywordtype}{real} :: sl\_kp1    \textcolor{comment}{! The sign-corrected slope of the interface below [Z L-1 ~> nondim].}
1529   \textcolor{keywordtype}{real} :: i\_sl\_k    \textcolor{comment}{! The (limited) inverse of sl\_K [L Z-1 ~> nondim].}
1530   \textcolor{keywordtype}{real} :: i\_sl\_kp1  \textcolor{comment}{! The (limited) inverse of sl\_Kp1 [L Z-1 ~> nondim].}
1531   \textcolor{keywordtype}{real} :: i\_4t      \textcolor{comment}{! A quarter of a flux scaling factor divided by}
1532                     \textcolor{comment}{! the damping timescale [T-1 ~> s-1].}
1533   \textcolor{keywordtype}{real} :: fn\_r      \textcolor{comment}{! A function of Rsl, such that Rsl < Fn\_R < 1.}
1534   \textcolor{keywordtype}{real} :: denom, i\_denom \textcolor{comment}{! A denominator and its inverse, various units.}
1535   \textcolor{keywordtype}{real} :: kh\_max    \textcolor{comment}{! A local ceiling on the diffusivity [L2 T-1 ~> m2 s-1].}
1536   \textcolor{keywordtype}{real} :: wt1, wt2  \textcolor{comment}{! Nondimensional weights.}
1537   \textcolor{comment}{!   Variables used only in testing code.}
1538   \textcolor{comment}{! real, dimension(SZK\_(G)) :: uh\_here}
1539   \textcolor{comment}{! real, dimension(SZK\_(G)+1) :: Sfn}
1540   \textcolor{keywordtype}{real} :: dkh       \textcolor{comment}{! An increment in the diffusivity [L2 T-1 ~> m2 s-1].}
1541 
1542   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(G)+1)} :: &
1543     kh\_bg, &        \textcolor{comment}{! The background (floor) value of Kh [L2 T-1 ~> m2 s-1].}
1544     kh, &           \textcolor{comment}{! The tentative value of Kh [L2 T-1 ~> m2 s-1].}
1545     kh\_detangle, &  \textcolor{comment}{! The detangling diffusivity that could be used [L2 T-1 ~> m2 s-1].}
1546     kh\_min\_max\_p, & \textcolor{comment}{! The smallest ceiling that can be placed on Kh(I,K)}
1547                     \textcolor{comment}{! based on the value of Kh(I,K+1) [L2 T-1 ~> m2 s-1].}
1548     kh\_min\_max\_m, & \textcolor{comment}{! The smallest ceiling that can be placed on Kh(I,K)}
1549                     \textcolor{comment}{! based on the value of Kh(I,K-1) [L2 T-1 ~> m2 s-1].}
1550     \textcolor{comment}{! The following are variables that define the relationships between}
1551     \textcolor{comment}{! successive values of Kh.}
1552     \textcolor{comment}{! Search for Kh that satisfy...}
1553     \textcolor{comment}{!    Kh(I,K) >= Kh\_min\_m(I,K)*Kh(I,K-1) + Kh0\_min\_m(I,K)}
1554     \textcolor{comment}{!    Kh(I,K) >= Kh\_min\_p(I,K)*Kh(I,K+1) + Kh0\_min\_p(I,K)}
1555     \textcolor{comment}{!    Kh(I,K) <= Kh\_max\_m(I,K)*Kh(I,K-1) + Kh0\_max\_m(I,K)}
1556     \textcolor{comment}{!    Kh(I,K) <= Kh\_max\_p(I,K)*Kh(I,K+1) + Kh0\_max\_p(I,K)}
1557     kh\_min\_m , &   \textcolor{comment}{! See above [nondim].}
1558     kh0\_min\_m , &  \textcolor{comment}{! See above [L2 T-1 ~> m2 s-1].}
1559     kh\_max\_m , &   \textcolor{comment}{! See above [nondim].}
1560     kh0\_max\_m, &   \textcolor{comment}{! See above [L2 T-1 ~> m2 s-1].}
1561     kh\_min\_p , &   \textcolor{comment}{! See above [nondim].}
1562     kh0\_min\_p , &  \textcolor{comment}{! See above [L2 T-1 ~> m2 s-1].}
1563     kh\_max\_p , &   \textcolor{comment}{! See above [nondim].}
1564     kh0\_max\_p      \textcolor{comment}{! See above [L2 T-1 ~> m2 s-1].}
1565   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
1566     kh\_max\_max  \textcolor{comment}{! The maximum diffusivity permitted in a column [L2 T-1 ~> m2 s-1]..}
1567   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
1568     do\_i        \textcolor{comment}{! If true, work on a column.}
1569   \textcolor{keywordtype}{integer} :: i, j, k, n, ish, jsh, is, ie, js, je, nz, k\_top
1570   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1571 
1572   k\_top = gv%nk\_rho\_varies + 1
1573   h\_neglect = gv%H\_subroundoff
1574   \textcolor{comment}{!   The 0.5 is because we are not using uniform weightings, but are}
1575   \textcolor{comment}{! distributing the diffusivities more effectively (with wt1 & wt2), but this}
1576   \textcolor{comment}{! means that the additions to a single interface can be up to twice as large.}
1577   kh\_scale = 0.5
1578   \textcolor{keywordflow}{if} (cs%detangle\_time > dt) kh\_scale = 0.5 * dt / cs%detangle\_time
1579 
1580   \textcolor{keywordflow}{do} j=js-1,je+1 ; \textcolor{keywordflow}{do} i=is-1,ie+1
1581     de\_top(i,j,k\_top) = 0.0 ; de\_bot(i,j) = 0.0
1582 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
1583   \textcolor{keywordflow}{do} k=k\_top+1,nz ; \textcolor{keywordflow}{do} j=js-1,je+1 ; \textcolor{keywordflow}{do} i=is-1,ie+1
1584     de\_top(i,j,k) = de\_top(i,j,k-1) + h(i,j,k-1)
1585 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1586 
1587   \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
1588     kh\_lay\_u(i,j,nz) = 0.0 ; kh\_lay\_u(i,j,k\_top) = 0.0
1589 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
1590   \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
1591     kh\_lay\_v(i,j,nz) = 0.0 ; kh\_lay\_v(i,j,k\_top) = 0.0
1592 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
1593 
1594   \textcolor{keywordflow}{do} k=nz-1,k\_top+1,-1
1595     \textcolor{comment}{! Find the diffusivities associated with each layer.}
1596     \textcolor{keywordflow}{do} j=js-1,je+1 ; \textcolor{keywordflow}{do} i=is-1,ie+1
1597       de\_bot(i,j) = de\_bot(i,j) + h(i,j,k+1)
1598 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
1599 
1600     \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie ; \textcolor{keywordflow}{if} (g%mask2dCu(i,j) > 0.0) \textcolor{keywordflow}{then}
1601       \textcolor{keywordflow}{if} (h(i,j,k) > h(i+1,j,k)) \textcolor{keywordflow}{then}
1602         h2 = h(i,j,k)
1603         h1 = max( h(i+1,j,k), h2 - min(de\_bot(i+1,j), de\_top(i+1,j,k)) )
1604       \textcolor{keywordflow}{else}
1605         h2 = h(i+1,j,k)
1606         h1 = max( h(i,j,k), h2 - min(de\_bot(i,j), de\_top(i,j,k)) )
1607 \textcolor{keywordflow}{      endif}
1608       jag\_rat = (h2 - h1)**2 / (h2 + h1 + h\_neglect)**2
1609       kh\_lay\_u(i,j,k) = (kh\_scale * kh\_u\_cfl(i,j)) * jag\_rat**2
1610 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1611 
1612     \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dCv(i,j) > 0.0) \textcolor{keywordflow}{then}
1613       \textcolor{keywordflow}{if} (h(i,j,k) > h(i,j+1,k)) \textcolor{keywordflow}{then}
1614         h2 = h(i,j,k)
1615         h1 = max( h(i,j+1,k), h2 - min(de\_bot(i,j+1), de\_top(i,j+1,k)) )
1616       \textcolor{keywordflow}{else}
1617         h2 = h(i,j+1,k)
1618         h1 = max( h(i,j,k), h2 - min(de\_bot(i,j), de\_top(i,j,k)) )
1619 \textcolor{keywordflow}{      endif}
1620       jag\_rat = (h2 - h1)**2 / (h2 + h1 + h\_neglect)**2
1621       kh\_lay\_v(i,j,k) = (kh\_scale * kh\_v\_cfl(i,j)) * jag\_rat**2
1622 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1623 \textcolor{keywordflow}{  enddo}
1624 
1625   \textcolor{comment}{! Limit the diffusivities}
1626 
1627   i\_4t = kh\_scale / (4.0 * dt)
1628 
1629   \textcolor{keywordflow}{do} n=1,2
1630     \textcolor{keywordflow}{if} (n==1) \textcolor{keywordflow}{then} ; jsh = js ; ish = is-1
1631     \textcolor{keywordflow}{else} ; jsh = js-1 ; ish = is ;\textcolor{keywordflow}{ endif}
1632 
1633     \textcolor{keywordflow}{do} j=jsh,je
1634 
1635       \textcolor{comment}{! First, populate the diffusivities}
1636       \textcolor{keywordflow}{if} (n==1) \textcolor{keywordflow}{then} \textcolor{comment}{! This is a u-column.}
1637         \textcolor{keywordflow}{do} i=ish,ie
1638           do\_i(i) = (g%mask2dCu(i,j) > 0.0)
1639           kh\_max\_max(i) = kh\_u\_cfl(i,j)
1640 \textcolor{keywordflow}{        enddo}
1641         \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} i=ish,ie
1642           kh\_bg(i,k) = kh\_u(i,j,k) ; kh(i,k) = kh\_bg(i,k)
1643           kh\_min\_max\_p(i,k) = kh\_bg(i,k) ; kh\_min\_max\_m(i,k) = kh\_bg(i,k)
1644           kh\_detangle(i,k) = 0.0
1645 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1646       \textcolor{keywordflow}{else} \textcolor{comment}{! This is a v-column.}
1647         \textcolor{keywordflow}{do} i=ish,ie
1648           do\_i(i) = (g%mask2dCv(i,j) > 0.0) ; kh\_max\_max(i) = kh\_v\_cfl(i,j)
1649 \textcolor{keywordflow}{        enddo}
1650         \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} i=ish,ie
1651           kh\_bg(i,k) = kh\_v(i,j,k) ; kh(i,k) = kh\_bg(i,k)
1652           kh\_min\_max\_p(i,k) = kh\_bg(i,k) ; kh\_min\_max\_m(i,k) = kh\_bg(i,k)
1653           kh\_detangle(i,k) = 0.0
1654 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1655 \textcolor{keywordflow}{      endif}
1656 
1657       \textcolor{comment}{! Determine the limits on the diffusivities.}
1658       \textcolor{keywordflow}{do} k=k\_top,nz ; \textcolor{keywordflow}{do} i=ish,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1659         \textcolor{keywordflow}{if} (n==1) \textcolor{keywordflow}{then} \textcolor{comment}{! This is a u-column.}
1660           dh = 0.0
1661           denom = ((g%IareaT(i+1,j) + g%IareaT(i,j)) * g%dy\_Cu(i,j))
1662           \textcolor{comment}{!   This expression uses differences in e in place of h for better}
1663           \textcolor{comment}{! consistency with the slopes.}
1664           \textcolor{keywordflow}{if} (denom > 0.0) &
1665             dh = i\_4t * ((e(i+1,j,k) - e(i+1,j,k+1)) - &
1666                          (e(i,j,k) - e(i,j,k+1))) / denom
1667            \textcolor{comment}{! dH = I\_4t * (h(i+1,j,k) - h(i,j,k)) / denom}
1668 
1669           adh = abs(dh)
1670           sign = 1.0 ; \textcolor{keywordflow}{if} (dh < 0) sign = -1.0
1671           sl\_k = sign * (e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j)
1672           sl\_kp1 = sign * (e(i+1,j,k+1)-e(i,j,k+1)) * g%IdxCu(i,j)
1673 
1674           \textcolor{comment}{! Add the incremental diffusivites to the surrounding interfaces.}
1675           \textcolor{comment}{! Adding more to the more steeply sloping layers (as below) makes}
1676           \textcolor{comment}{! the diffusivities more than twice as effective.}
1677           denom = (sl\_k**2 + sl\_kp1**2)
1678           wt1 = 0.5 ; wt2 = 0.5
1679           \textcolor{keywordflow}{if} (denom > 0.0) \textcolor{keywordflow}{then}
1680             wt1 = sl\_k**2 / denom ; wt2 = sl\_kp1**2 / denom
1681 \textcolor{keywordflow}{          endif}
1682           kh\_detangle(i,k) = kh\_detangle(i,k) + wt1*kh\_lay\_u(i,j,k)
1683           kh\_detangle(i,k+1) = kh\_detangle(i,k+1) + wt2*kh\_lay\_u(i,j,k)
1684         \textcolor{keywordflow}{else} \textcolor{comment}{! This is a v-column.}
1685           dh = 0.0
1686           denom = ((g%IareaT(i,j+1) + g%IareaT(i,j)) * g%dx\_Cv(i,j))
1687           \textcolor{keywordflow}{if} (denom > 0.0) &
1688             dh = i\_4t * ((e(i,j+1,k) - e(i,j+1,k+1)) - &
1689                          (e(i,j,k) - e(i,j,k+1))) / denom
1690            \textcolor{comment}{! dH = I\_4t * (h(i,j+1,k) - h(i,j,k)) / denom}
1691 
1692           adh = abs(dh)
1693           sign = 1.0 ; \textcolor{keywordflow}{if} (dh < 0) sign = -1.0
1694           sl\_k = sign * (e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j)
1695           sl\_kp1 = sign * (e(i,j+1,k+1)-e(i,j,k+1)) * g%IdyCv(i,j)
1696 
1697           \textcolor{comment}{! Add the incremental diffusviites to the surrounding interfaces.}
1698           \textcolor{comment}{! Adding more to the more steeply sloping layers (as below) makes}
1699           \textcolor{comment}{! the diffusivities more than twice as effective.}
1700           denom = (sl\_k**2 + sl\_kp1**2)
1701           wt1 = 0.5 ; wt2 = 0.5
1702           \textcolor{keywordflow}{if} (denom > 0.0) \textcolor{keywordflow}{then}
1703             wt1 = sl\_k**2 / denom ; wt2 = sl\_kp1**2 / denom
1704 \textcolor{keywordflow}{          endif}
1705           kh\_detangle(i,k) = kh\_detangle(i,k) + wt1*kh\_lay\_v(i,j,k)
1706           kh\_detangle(i,k+1) = kh\_detangle(i,k+1) + wt2*kh\_lay\_v(i,j,k)
1707 \textcolor{keywordflow}{        endif}
1708 
1709         \textcolor{keywordflow}{if} (adh == 0.0) \textcolor{keywordflow}{then}
1710           kh\_min\_m(i,k+1) = 1.0 ; kh0\_min\_m(i,k+1) = 0.0
1711           kh\_max\_m(i,k+1) = 1.0 ; kh0\_max\_m(i,k+1) = 0.0
1712           kh\_min\_p(i,k) = 1.0 ; kh0\_min\_p(i,k) = 0.0
1713           kh\_max\_p(i,k) = 1.0 ; kh0\_max\_p(i,k) = 0.0
1714         \textcolor{keywordflow}{elseif} (adh > 0.0) \textcolor{keywordflow}{then}
1715           \textcolor{keywordflow}{if} (sl\_k <= sl\_kp1) \textcolor{keywordflow}{then}
1716             \textcolor{comment}{! This case should only arise from nonlinearities in the equation of state.}
1717             \textcolor{comment}{! Treat it as though dedx(K) = dedx(K+1) & dH = 0.}
1718             kh\_min\_m(i,k+1) = 1.0 ; kh0\_min\_m(i,k+1) = 0.0
1719             kh\_max\_m(i,k+1) = 1.0 ; kh0\_max\_m(i,k+1) = 0.0
1720             kh\_min\_p(i,k) = 1.0 ; kh0\_min\_p(i,k) = 0.0
1721             kh\_max\_p(i,k) = 1.0 ; kh0\_max\_p(i,k) = 0.0
1722           \textcolor{keywordflow}{elseif} (sl\_k <= 0.0) \textcolor{keywordflow}{then}   \textcolor{comment}{! Both slopes are opposite to dH}
1723             i\_sl = -1.0 / sl\_kp1
1724             rsl = -sl\_k * i\_sl                            \textcolor{comment}{! 0 <= Rsl < 1}
1725             irsl = 1e9 ; \textcolor{keywordflow}{if} (rsl > 1e-9) irsl = 1.0/rsl   \textcolor{comment}{! 1 < IRsl <= 1e9}
1726 
1727             fn\_r = rsl
1728             \textcolor{keywordflow}{if} (kh\_max\_max(i) > 0) &
1729               fn\_r = min(sqrt(rsl), rsl + (adh * i\_sl) / (kh\_max\_max(i)))
1730 
1731             kh\_min\_m(i,k+1) = fn\_r ; kh0\_min\_m(i,k+1) = 0.0
1732             kh\_max\_m(i,k+1) = rsl ; kh0\_max\_m(i,k+1) = adh * i\_sl
1733             kh\_min\_p(i,k) = irsl ; kh0\_min\_p(i,k) = -adh * (i\_sl*irsl)
1734             kh\_max\_p(i,k) = 1.0/(fn\_r + 1.0e-30) ; kh0\_max\_p(i,k) = 0.0
1735           \textcolor{keywordflow}{elseif} (sl\_kp1 < 0.0) \textcolor{keywordflow}{then}  \textcolor{comment}{! Opposite (nonzero) signs of slopes.}
1736             i\_sl\_k = 1e18*us%Z\_to\_L ; \textcolor{keywordflow}{if} (sl\_k > 1e-18*us%L\_to\_Z) i\_sl\_k = 1.0 / sl\_k
1737             i\_sl\_kp1 = 1e18*us%Z\_to\_L ; \textcolor{keywordflow}{if} (-sl\_kp1 > 1e-18*us%L\_to\_Z) i\_sl\_kp1 = -1.0 / sl\_kp1
1738 
1739             kh\_min\_m(i,k+1) = 0.0 ; kh0\_min\_m(i,k+1) = 0.0
1740             kh\_max\_m(i,k+1) = - sl\_k*i\_sl\_kp1 ; kh0\_max\_m(i,k+1) = adh*i\_sl\_kp1
1741             kh\_min\_p(i,k) = 0.0 ; kh0\_min\_p(i,k) = 0.0
1742             kh\_max\_p(i,k) = sl\_kp1*i\_sl\_k ; kh0\_max\_p(i,k) = adh*i\_sl\_k
1743 
1744             \textcolor{comment}{! This limit does not use the slope weighting so that potentially}
1745             \textcolor{comment}{! sharp gradients in diffusivities are not forced to occur.}
1746             kh\_max = adh / (sl\_k - sl\_kp1)
1747             kh\_min\_max\_p(i,k) = max(kh\_min\_max\_p(i,k), kh\_max)
1748             kh\_min\_max\_m(i,k+1) = max(kh\_min\_max\_m(i,k+1), kh\_max)
1749           \textcolor{keywordflow}{else} \textcolor{comment}{! Both slopes are of the same sign as dH.}
1750             i\_sl = 1.0 / sl\_k
1751             rsl = sl\_kp1 * i\_sl                           \textcolor{comment}{! 0 <= Rsl < 1}
1752             irsl = 1e9 ; \textcolor{keywordflow}{if} (rsl > 1e-9) irsl = 1.0/rsl   \textcolor{comment}{! 1 < IRsl <= 1e9}
1753 
1754             \textcolor{comment}{! Rsl <= Fn\_R <= 1}
1755             fn\_r = rsl
1756             \textcolor{keywordflow}{if} (kh\_max\_max(i) > 0) &
1757               fn\_r = min(sqrt(rsl), rsl + (adh * i\_sl) / kh\_max\_max(i))
1758 
1759             kh\_min\_m(i,k+1) = irsl ; kh0\_min\_m(i,k+1) = -adh * (i\_sl*irsl)
1760             kh\_max\_m(i,k+1) = 1.0/(fn\_r + 1.0e-30) ; kh0\_max\_m(i,k+1) = 0.0
1761             kh\_min\_p(i,k) = fn\_r ; kh0\_min\_p(i,k) = 0.0
1762             kh\_max\_p(i,k) = rsl ; kh0\_max\_p(i,k) = adh * i\_sl
1763 \textcolor{keywordflow}{          endif}
1764 \textcolor{keywordflow}{        endif}
1765 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! I-loop & k-loop}
1766 
1767       \textcolor{keywordflow}{do} k=k\_top,nz+1,nz+1-k\_top ; \textcolor{keywordflow}{do} i=ish,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1768         \textcolor{comment}{! The diffusivities at k\_top and nz+1 are both fixed.}
1769         kh\_min\_m(i,k) = 0.0 ; kh0\_min\_m(i,k) = 0.0
1770         kh\_max\_m(i,k) = 0.0 ; kh0\_max\_m(i,k) = 0.0
1771         kh\_min\_p(i,k) = 0.0 ; kh0\_min\_p(i,k) = 0.0
1772         kh\_max\_p(i,k) = 0.0 ; kh0\_max\_p(i,k) = 0.0
1773         kh\_min\_max\_p(i,k) = kh\_bg(i,k)
1774         kh\_min\_max\_m(i,k) = kh\_bg(i,k)
1775 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! I-loop and k\_top/nz+1 loop}
1776 
1777       \textcolor{comment}{! Search for Kh that satisfy...}
1778       \textcolor{comment}{!    Kh(I,K) >= Kh\_min\_m(I,K)*Kh(I,K-1) + Kh0\_min\_m(I,K)}
1779       \textcolor{comment}{!    Kh(I,K) >= Kh\_min\_p(I,K)*Kh(I,K+1) + Kh0\_min\_p(I,K)}
1780       \textcolor{comment}{!    Kh(I,K) <= Kh\_max\_m(I,K)*Kh(I,K-1) + Kh0\_max\_m(I,K)}
1781       \textcolor{comment}{!    Kh(I,K) <= Kh\_max\_p(I,K)*Kh(I,K+1) + Kh0\_max\_p(I,K)}
1782 
1783       \textcolor{comment}{! Increase the diffusivies to satisfy the min constraints.}
1784       \textcolor{comment}{! All non-zero min constraints on one diffusivity are max constraints on another.}
1785       \textcolor{keywordflow}{do} k=k\_top+1,nz ; \textcolor{keywordflow}{do} i=ish,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1786         kh(i,k) = max(kh\_bg(i,k), kh\_detangle(i,k), &
1787                       min(kh\_min\_m(i,k)*kh(i,k-1) + kh0\_min\_m(i,k), kh(i,k-1)))
1788 
1789         \textcolor{keywordflow}{if} (kh0\_max\_m(i,k) > kh\_bg(i,k)) kh(i,k) = min(kh(i,k), kh0\_max\_m(i,k))
1790         \textcolor{keywordflow}{if} (kh0\_max\_p(i,k) > kh\_bg(i,k)) kh(i,k) = min(kh(i,k), kh0\_max\_p(i,k))
1791 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! I-loop & k-loop}
1792       \textcolor{comment}{! This is still true... do i=ish,ie ; Kh(I,nz+1) = Kh\_bg(I,nz+1) ; enddo}
1793       \textcolor{keywordflow}{do} k=nz,k\_top+1,-1 ; \textcolor{keywordflow}{do} i=ish,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1794         kh(i,k) = max(kh(i,k), min(kh\_min\_p(i,k)*kh(i,k+1) + kh0\_min\_p(i,k), kh(i,k+1)))
1795 
1796         kh\_max = max(kh\_min\_max\_p(i,k), kh\_max\_p(i,k)*kh(i,k+1) + kh0\_max\_p(i,k))
1797         kh(i,k) = min(kh(i,k), kh\_max)
1798 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! I-loop & k-loop}
1799       \textcolor{comment}{!  All non-zero min constraints on one diffusivity are max constraints on}
1800       \textcolor{comment}{! another layer, so the min constraints can now be discounted.}
1801 
1802       \textcolor{comment}{! Decrease the diffusivities to satisfy the max constraints.}
1803         \textcolor{keywordflow}{do} k=k\_top+1,nz ; \textcolor{keywordflow}{do} i=ish,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1804           kh\_max = max(kh\_min\_max\_m(i,k), kh\_max\_m(i,k)*kh(i,k-1) + kh0\_max\_m(i,k))
1805           \textcolor{keywordflow}{if} (kh(i,k) > kh\_max) kh(i,k) = kh\_max
1806 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}  \textcolor{comment}{! i- and K-loops}
1807 
1808       \textcolor{comment}{! This code tests the solutions...}
1809 \textcolor{comment}{!     do i=ish,ie}
1810 \textcolor{comment}{!       Sfn(:) = 0.0 ; uh\_here(:) = 0.0}
1811 \textcolor{comment}{!       do K=k\_top,nz}
1812 \textcolor{comment}{!         if ((Kh(i,K) > Kh\_bg(i,K)) .or. (Kh(i,K+1) > Kh\_bg(i,K+1))) then}
1813 \textcolor{comment}{!           if (n==1) then ! u-point.}
1814 \textcolor{comment}{!             if ((h(i+1,j,k) - h(i,j,k)) * &}
1815 \textcolor{comment}{!                 ((e(i+1,j,K)-e(i+1,j,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then}
1816 \textcolor{comment}{!               Sfn(K) = -Kh(i,K) * (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)}
1817 \textcolor{comment}{!               Sfn(K+1) = -Kh(i,K+1) * (e(i+1,j,K+1)-e(i,j,K+1)) * G%IdxCu(I,j)}
1818 \textcolor{comment}{!               uh\_here(k) = (Sfn(K) - Sfn(K+1))*G%dy\_Cu(I,j)}
1819 \textcolor{comment}{!               if (abs(uh\_here(k)) * min(G%IareaT(i,j), G%IareaT(i+1,j)) > &}
1820 \textcolor{comment}{!                   (1e-10*GV%m\_to\_H)) then}
1821 \textcolor{comment}{!                 if (uh\_here(k) * (h(i+1,j,k) - h(i,j,k)) > 0.0) then}
1822 \textcolor{comment}{!                   call MOM\_error(WARNING, "Corrective u-transport is up the thickness gradient.", .true.)}
1823 \textcolor{comment}{!                 endif}
1824 \textcolor{comment}{!                 if (((h(i,j,k) - 4.0*dt*G%IareaT(i,j)*uh\_here(k)) - &}
1825 \textcolor{comment}{!                      (h(i+1,j,k) + 4.0*dt*G%IareaT(i+1,j)*uh\_here(k))) * &}
1826 \textcolor{comment}{!                     (h(i,j,k) - h(i+1,j,k)) < 0.0) then}
1827 \textcolor{comment}{!                   call MOM\_error(WARNING, "Corrective u-transport is too large.", .true.)}
1828 \textcolor{comment}{!                 endif}
1829 \textcolor{comment}{!               endif}
1830 \textcolor{comment}{!             endif}
1831 \textcolor{comment}{!           else ! v-point}
1832 \textcolor{comment}{!             if ((h(i,j+1,k) - h(i,j,k)) * &}
1833 \textcolor{comment}{!                 ((e(i,j+1,K)-e(i,j+1,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then}
1834 \textcolor{comment}{!               Sfn(K) = -Kh(i,K) * (e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)}
1835 \textcolor{comment}{!               Sfn(K+1) = -Kh(i,K+1) * (e(i,j+1,K+1)-e(i,j,K+1)) * G%IdyCv(i,J)}
1836 \textcolor{comment}{!               uh\_here(k) = (Sfn(K) - Sfn(K+1))*G%dx\_Cv(i,J)}
1837 \textcolor{comment}{!               if (abs(uh\_here(K)) * min(G%IareaT(i,j), G%IareaT(i,j+1)) > &}
1838 \textcolor{comment}{!                   (1e-10*GV%m\_to\_H)) then}
1839 \textcolor{comment}{!                 if (uh\_here(K) * (h(i,j+1,k) - h(i,j,k)) > 0.0) then}
1840 \textcolor{comment}{!                   call MOM\_error(WARNING, &}
1841 \textcolor{comment}{!                          "Corrective v-transport is up the thickness gradient.", .true.)}
1842 \textcolor{comment}{!                 endif}
1843 \textcolor{comment}{!                 if (((h(i,j,k) - 4.0*dt*G%IareaT(i,j)*uh\_here(K)) - &}
1844 \textcolor{comment}{!                      (h(i,j+1,k) + 4.0*dt*G%IareaT(i,j+1)*uh\_here(K))) * &}
1845 \textcolor{comment}{!                     (h(i,j,k) - h(i,j+1,k)) < 0.0) then}
1846 \textcolor{comment}{!                   call MOM\_error(WARNING, &}
1847 \textcolor{comment}{!                          "Corrective v-transport is too large.", .true.)}
1848 \textcolor{comment}{!                 endif}
1849 \textcolor{comment}{!               endif}
1850 \textcolor{comment}{!             endif}
1851 \textcolor{comment}{!           endif ! u- or v- selection.}
1852 \textcolor{comment}{!          !  de\_dx(I,K) = (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)}
1853 \textcolor{comment}{!         endif}
1854 \textcolor{comment}{!       enddo}
1855 \textcolor{comment}{!     enddo}
1856 
1857       \textcolor{keywordflow}{if} (n==1) \textcolor{keywordflow}{then} \textcolor{comment}{! This is a u-column.}
1858         \textcolor{keywordflow}{do} k=k\_top+1,nz ; \textcolor{keywordflow}{do} i=ish,ie
1859           \textcolor{keywordflow}{if} (kh(i,k) > kh\_u(i,j,k)) \textcolor{keywordflow}{then}
1860             dkh = (kh(i,k) - kh\_u(i,j,k))
1861             int\_slope\_u(i,j,k) = dkh / kh(i,k)
1862             kh\_u(i,j,k) = kh(i,k)
1863 \textcolor{keywordflow}{          endif}
1864 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1865       \textcolor{keywordflow}{else} \textcolor{comment}{! This is a v-column.}
1866         \textcolor{keywordflow}{do} k=k\_top+1,nz ; \textcolor{keywordflow}{do} i=ish,ie
1867           \textcolor{keywordflow}{if} (kh(i,k) > kh\_v(i,j,k)) \textcolor{keywordflow}{then}
1868             dkh = kh(i,k) - kh\_v(i,j,k)
1869             int\_slope\_v(i,j,k) = dkh / kh(i,k)
1870             kh\_v(i,j,k) = kh(i,k)
1871 \textcolor{keywordflow}{          endif}
1872 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1873 \textcolor{keywordflow}{      endif}
1874 
1875 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! j-loop}
1876 \textcolor{keywordflow}{  enddo}  \textcolor{comment}{! n-loop over u- and v- directions.}
1877 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__thickness__diffuse_a52d5fe57d53414fdc05f669723c9774e}\label{namespacemom__thickness__diffuse_a52d5fe57d53414fdc05f669723c9774e}} 
\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}!streamfn\+\_\+solver@{streamfn\+\_\+solver}}
\index{streamfn\+\_\+solver@{streamfn\+\_\+solver}!mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}
\subsubsection{\texorpdfstring{streamfn\+\_\+solver()}{streamfn\_solver()}}
{\footnotesize\ttfamily subroutine mom\+\_\+thickness\+\_\+diffuse\+::streamfn\+\_\+solver (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{nk,  }\item[{real, dimension(nk), intent(in)}]{c2\+\_\+h,  }\item[{real, dimension(nk+1), intent(in)}]{h\+N2,  }\item[{real, dimension(nk+1), intent(inout)}]{sfn }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Tridiagonal solver for streamfunction at interfaces. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em nk} & Number of layers\\
\hline
\mbox{\tt in}  & {\em c2\+\_\+h} & Wave speed squared over thickness in layers \mbox{[}L2 Z-\/1 T-\/2 $\sim$$>$ m s-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em hn2} & Thickness times N2 at interfaces \mbox{[}L2 Z-\/1 T-\/2 $\sim$$>$ m s-\/2\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em sfn} & Streamfunction \mbox{[}Z L2 T-\/1 $\sim$$>$ m3 s-\/1\mbox{]} or arbitrary units On entry, equals diffusivity times slope. On exit, equals the streamfunction. \\
\hline
\end{DoxyParams}


Definition at line 1434 of file M\+O\+M\+\_\+thickness\+\_\+diffuse.\+F90.


\begin{DoxyCode}
1434   \textcolor{keywordtype}{integer},               \textcolor{keywordtype}{intent(in)}    :: nk\textcolor{comment}{   !< Number of layers}
1435   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(nk)},   \textcolor{keywordtype}{intent(in)}    :: c2\_h\textcolor{comment}{ !< Wave speed squared over thickness in layers [L2 Z-1 T-2
       ~> m s-2]}
1436   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(nk+1)}, \textcolor{keywordtype}{intent(in)}    :: hn2\textcolor{comment}{  !< Thickness times N2 at interfaces [L2 Z-1 T-2 ~> m s-2]}
1437   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(nk+1)}, \textcolor{keywordtype}{intent(inout)} :: sfn\textcolor{comment}{  !< Streamfunction [Z L2 T-1 ~> m3 s-1] or arbitrary units}
1438 \textcolor{comment}{                                               !! On entry, equals diffusivity times slope.}
1439 \textcolor{comment}{                                               !! On exit, equals the streamfunction.}
1440   \textcolor{comment}{! Local variables}
1441   \textcolor{keywordtype}{integer} :: k
1442 
1443   \textcolor{keywordtype}{real} :: b\_denom, beta, d1, c1(nk)
1444 
1445   sfn(1) = 0.
1446   b\_denom = hn2(2) + c2\_h(1)
1447   beta = 1.0 / ( b\_denom + c2\_h(2) )
1448   d1 = beta * b\_denom
1449   sfn(2) = ( beta * hn2(2) )*sfn(2)
1450   \textcolor{keywordflow}{do} k=3,nk
1451     c1(k-1) = beta * c2\_h(k-1)
1452     b\_denom = hn2(k) + d1*c2\_h(k-1)
1453     beta = 1.0 / (b\_denom + c2\_h(k))
1454     d1 = beta * b\_denom
1455     sfn(k) = beta * (hn2(k)*sfn(k) + c2\_h(k-1)*sfn(k-1))
1456 \textcolor{keywordflow}{  enddo}
1457   c1(nk) = beta * c2\_h(nk)
1458   sfn(nk+1) = 0.
1459   \textcolor{keywordflow}{do} k=nk,2,-1
1460     sfn(k) = sfn(k) + c1(k)*sfn(k+1)
1461 \textcolor{keywordflow}{  enddo}
1462 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__thickness__diffuse_a8a538b778a567f489bfd9c5eadeeebef}\label{namespacemom__thickness__diffuse_a8a538b778a567f489bfd9c5eadeeebef}} 
\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}!thickness\+\_\+diffuse@{thickness\+\_\+diffuse}}
\index{thickness\+\_\+diffuse@{thickness\+\_\+diffuse}!mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}
\subsubsection{\texorpdfstring{thickness\+\_\+diffuse()}{thickness\_diffuse()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+thickness\+\_\+diffuse\+::thickness\+\_\+diffuse (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{h,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{uhtr,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)), intent(inout)}]{vhtr,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(meke\+\_\+type), pointer}]{M\+E\+KE,  }\item[{type(varmix\+\_\+cs), pointer}]{Var\+Mix,  }\item[{type(cont\+\_\+diag\+\_\+ptrs), intent(inout)}]{C\+Dp,  }\item[{type(\hyperlink{structmom__thickness__diffuse_1_1thickness__diffuse__cs}{thickness\+\_\+diffuse\+\_\+cs}), pointer}]{CS }\end{DoxyParamCaption})}



Calculates thickness diffusion coefficients and applies thickness diffusion to layer thicknesses, h. Diffusivities are limited to ensure stability. Also returns along-\/layer mass fluxes used in the continuity equation. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in,out}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em uhtr} & Accumulated zonal mass flux \mbox{[}L2 H $\sim$$>$ m3 or kg\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em vhtr} & Accumulated meridional mass flux \mbox{[}L2 H $\sim$$>$ m3 or kg\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & Thermodynamics structure\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
 & {\em meke} & M\+E\+KE control structure\\
\hline
 & {\em varmix} & Variable mixing coefficients\\
\hline
\mbox{\tt in,out}  & {\em cdp} & Diagnostics for the continuity equation\\
\hline
 & {\em cs} & Control structure for thickness diffusion \\
\hline
\end{DoxyParams}


Definition at line 109 of file M\+O\+M\+\_\+thickness\+\_\+diffuse.\+F90.


\begin{DoxyCode}
109   \textcolor{keywordtype}{type}(ocean\_grid\_type),                     \textcolor{keywordtype}{intent(in)}    :: g\textcolor{comment}{      !< Ocean grid structure}
110   \textcolor{keywordtype}{type}(verticalgrid\_type),                   \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{     !< Vertical grid structure}
111   \textcolor{keywordtype}{type}(unit\_scale\_type),                     \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{     !< A dimensional unit scaling type}
112   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))},  \textcolor{keywordtype}{intent(inout)} :: h\textcolor{comment}{      !< Layer thickness [H ~> m or kg m-2]}
113   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: uhtr\textcolor{comment}{   !< Accumulated zonal mass flux}
114 \textcolor{comment}{                                                                     !! [L2 H ~> m3 or kg]}
115   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))}, \textcolor{keywordtype}{intent(inout)} :: vhtr\textcolor{comment}{   !< Accumulated meridional mass flux}
116 \textcolor{comment}{                                                                     !! [L2 H ~> m3 or kg]}
117   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),                     \textcolor{keywordtype}{intent(in)}    :: tv\textcolor{comment}{     !< Thermodynamics structure}
118   \textcolor{keywordtype}{real},                                      \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{     !< Time increment [T ~> s]}
119   \textcolor{keywordtype}{type}(meke\_type),                           \textcolor{keywordtype}{pointer}       :: meke\textcolor{comment}{   !< MEKE control structure}
120   \textcolor{keywordtype}{type}(varmix\_cs),                           \textcolor{keywordtype}{pointer}       :: varmix\textcolor{comment}{ !< Variable mixing coefficients}
121   \textcolor{keywordtype}{type}(cont\_diag\_ptrs),                      \textcolor{keywordtype}{intent(inout)} :: cdp\textcolor{comment}{    !< Diagnostics for the continuity
       equation}
122   \textcolor{keywordtype}{type}(thickness\_diffuse\_cs),                \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{     !< Control structure for thickness
       diffusion}
123   \textcolor{comment}{! Local variables}
124   \textcolor{keywordtype}{real} :: e(szi\_(g), szj\_(g), szk\_(g)+1) \textcolor{comment}{! heights of interfaces, relative to mean}
125                                          \textcolor{comment}{! sea level [Z ~> m], positive up.}
126   \textcolor{keywordtype}{real} :: uhd(szib\_(g), szj\_(g), szk\_(g)) \textcolor{comment}{! Diffusive u*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]}
127   \textcolor{keywordtype}{real} :: vhd(szi\_(g), szjb\_(g), szk\_(g)) \textcolor{comment}{! Diffusive v*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]}
128 
129   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G), SZJ\_(G), SZK\_(G)+1)} :: &
130     kh\_u, &       \textcolor{comment}{! interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1]}
131     int\_slope\_u   \textcolor{comment}{! A nondimensional ratio from 0 to 1 that gives the relative}
132                   \textcolor{comment}{! weighting of the interface slopes to that calculated also}
133                   \textcolor{comment}{! using density gradients at u points.  The physically correct}
134                   \textcolor{comment}{! slopes occur at 0, while 1 is used for numerical closures [nondim].}
135   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJB\_(G), SZK\_(G)+1)} :: &
136     kh\_v, &       \textcolor{comment}{! interface height diffusivities in v-columns [L2 T-1 ~> m2 s-1]}
137     int\_slope\_v   \textcolor{comment}{! A nondimensional ratio from 0 to 1 that gives the relative}
138                   \textcolor{comment}{! weighting of the interface slopes to that calculated also}
139                   \textcolor{comment}{! using density gradients at v points.  The physically correct}
140                   \textcolor{comment}{! slopes occur at 0, while 1 is used for numerical closures [nondim].}
141   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G), SZK\_(G))} :: &
142     kh\_t          \textcolor{comment}{! diagnosed diffusivity at tracer points [L2 T-1 ~> m2 s-1]}
143 
144   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G), SZJ\_(G))} :: &
145     kh\_u\_cfl      \textcolor{comment}{! The maximum stable interface height diffusivity at u grid points [L2 T-1 ~> m2 s-1]}
146   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJB\_(G))} :: &
147     kh\_v\_cfl      \textcolor{comment}{! The maximum stable interface height diffusivity at v grid points [L2 T-1 ~> m2 s-1]}
148   \textcolor{keywordtype}{real} :: khth\_loc\_u(szib\_(g), szj\_(g))
149   \textcolor{keywordtype}{real} :: khth\_loc\_v(szi\_(g), szjb\_(g))
150   \textcolor{keywordtype}{real} :: khth\_loc(szib\_(g), szjb\_(g))  \textcolor{comment}{! locally calculated thickness diffusivity [L2 T-1 ~> m2 s-1]}
151   \textcolor{keywordtype}{real} :: h\_neglect \textcolor{comment}{! A thickness that is so small it is usually lost}
152                     \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
153   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{pointer} :: cg1 => null() \textcolor{comment}{!< Wave speed [L T-1 ~> m s-1]}
154   \textcolor{keywordtype}{logical} :: use\_varmix, resoln\_scaled, depth\_scaled, use\_stored\_slopes, khth\_use\_ebt\_struct, use\_visbeck
155   \textcolor{keywordtype}{logical} :: use\_qg\_leith
156   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, nz
157   \textcolor{keywordtype}{real} :: hu(szi\_(g), szj\_(g))       \textcolor{comment}{! u-thickness [H ~> m or kg m-2]}
158   \textcolor{keywordtype}{real} :: hv(szi\_(g), szj\_(g))       \textcolor{comment}{! v-thickness [H ~> m or kg m-2]}
159   \textcolor{keywordtype}{real} :: kh\_u\_lay(szi\_(g), szj\_(g)) \textcolor{comment}{! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1]}
160   \textcolor{keywordtype}{real} :: kh\_v\_lay(szi\_(g), szj\_(g)) \textcolor{comment}{! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1]}
161 
162   \textcolor{keywordflow}{if} (.not. \textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"MOM\_thickness\_diffuse: "}//&
163          \textcolor{stringliteral}{"Module must be initialized before it is used."})
164 
165   \textcolor{keywordflow}{if} ((.not.cs%thickness\_diffuse) .or. &
166        .not.( cs%Khth > 0.0 .or. \textcolor{keyword}{associated}(varmix) .or. \textcolor{keyword}{associated}(meke) ) ) \textcolor{keywordflow}{return}
167 
168   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
169   h\_neglect = gv%H\_subroundoff
170 
171   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke)) \textcolor{keywordflow}{then}
172     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke%GM\_src)) \textcolor{keywordflow}{then}
173       \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie ; meke%GM\_src(i,j) = 0. ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
174 \textcolor{keywordflow}{    endif}
175 \textcolor{keywordflow}{  endif}
176 
177   use\_varmix = .false. ; resoln\_scaled = .false. ; use\_stored\_slopes = .false.
178   khth\_use\_ebt\_struct = .false. ; use\_visbeck = .false. ; use\_qg\_leith = .false.
179   depth\_scaled = .false.
180 
181   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(varmix)) \textcolor{keywordflow}{then}
182     use\_varmix = varmix%use\_variable\_mixing .and. (cs%KHTH\_Slope\_Cff > 0.)
183     resoln\_scaled = varmix%Resoln\_scaled\_KhTh
184     depth\_scaled = varmix%Depth\_scaled\_KhTh
185     use\_stored\_slopes = varmix%use\_stored\_slopes
186     khth\_use\_ebt\_struct = varmix%khth\_use\_ebt\_struct
187     use\_visbeck = varmix%use\_Visbeck
188     use\_qg\_leith = varmix%use\_QG\_Leith\_GM
189     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(varmix%cg1)) cg1 => varmix%cg1
190   \textcolor{keywordflow}{else}
191     cg1 => null()
192 \textcolor{keywordflow}{  endif}
193 
194 
195 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,KH\_u\_CFL,dt,G,CS)}
196   \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
197     kh\_u\_cfl(i,j) = (0.25*cs%max\_Khth\_CFL) /  &
198       (dt * (g%IdxCu(i,j)*g%IdxCu(i,j) + g%IdyCu(i,j)*g%IdyCu(i,j)))
199 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
200 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,KH\_v\_CFL,dt,G,CS)}
201   \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
202     kh\_v\_cfl(i,j) = (0.25*cs%max\_Khth\_CFL) / &
203       (dt * (g%IdxCv(i,j)*g%IdxCv(i,j) + g%IdyCv(i,j)*g%IdyCv(i,j)))
204 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
205 
206   \textcolor{comment}{! Calculates interface heights, e, in [Z ~> m].}
207   \textcolor{keyword}{call }find\_eta(h, tv, g, gv, us, e, halo\_size=1)
208 
209   \textcolor{comment}{! Set the diffusivities.}
210 \textcolor{comment}{!$OMP parallel default(none) shared(is,ie,js,je,Khth\_Loc\_u,CS,use\_VarMix,VarMix,    &}
211 \textcolor{comment}{!$OMP                               MEKE,Resoln\_scaled,KH\_u,G,use\_QG\_Leith,use\_Visbeck,&}
212 \textcolor{comment}{!$OMP                               KH\_u\_CFL,nz,Khth\_Loc,KH\_v,KH\_v\_CFL,int\_slope\_u, &}
213 \textcolor{comment}{!$OMP                               int\_slope\_v,khth\_use\_ebt\_struct, Depth\_scaled, &}
214 \textcolor{comment}{!$OMP                               Khth\_loc\_v)}
215 \textcolor{comment}{!$OMP do}
216   \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is-1,ie
217     khth\_loc\_u(i,j) = cs%Khth
218 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
219 
220   \textcolor{keywordflow}{if} (use\_varmix) \textcolor{keywordflow}{then}
221     \textcolor{keywordflow}{if} (use\_visbeck) \textcolor{keywordflow}{then}
222 \textcolor{comment}{!$OMP do}
223       \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
224         khth\_loc\_u(i,j) = khth\_loc\_u(i,j) + &
225           cs%KHTH\_Slope\_Cff*varmix%L2u(i,j) * varmix%SN\_u(i,j)
226 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
227 \textcolor{keywordflow}{    endif}
228 \textcolor{keywordflow}{  endif}
229 
230   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke%Kh)) \textcolor{keywordflow}{then}
231     \textcolor{keywordflow}{if} (cs%MEKE\_GEOMETRIC) \textcolor{keywordflow}{then}
232 \textcolor{comment}{!$OMP do}
233       \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
234         khth\_loc\_u(i,j) = khth\_loc\_u(i,j) + g%mask2dCu(i,j) * cs%MEKE\_GEOMETRIC\_alpha * &
235                           0.5*(meke%MEKE(i,j)+meke%MEKE(i+1,j)) / &
236                           (varmix%SN\_u(i,j) + cs%MEKE\_GEOMETRIC\_epsilon)
237 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
238     \textcolor{keywordflow}{else}
239       \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
240         khth\_loc\_u(i,j) = khth\_loc\_u(i,j) + meke%KhTh\_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
241 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
242 \textcolor{keywordflow}{    endif}
243 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
244 
245   \textcolor{keywordflow}{if} (resoln\_scaled) \textcolor{keywordflow}{then}
246 \textcolor{comment}{!$OMP do}
247     \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is-1,ie
248       khth\_loc\_u(i,j) = khth\_loc\_u(i,j) * varmix%Res\_fn\_u(i,j)
249 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
250 \textcolor{keywordflow}{  endif}
251 
252   \textcolor{keywordflow}{if} (depth\_scaled) \textcolor{keywordflow}{then}
253 \textcolor{comment}{!$OMP do}
254     \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is-1,ie
255       khth\_loc\_u(i,j) = khth\_loc\_u(i,j) * varmix%Depth\_fn\_u(i,j)
256 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
257 \textcolor{keywordflow}{  endif}
258 
259   \textcolor{keywordflow}{if} (cs%Khth\_Max > 0) \textcolor{keywordflow}{then}
260 \textcolor{comment}{!$OMP do}
261     \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is-1,ie
262       khth\_loc\_u(i,j) = max(cs%Khth\_Min, min(khth\_loc\_u(i,j), cs%Khth\_Max))
263 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
264   \textcolor{keywordflow}{else}
265 \textcolor{comment}{!$OMP do}
266     \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is-1,ie
267       khth\_loc\_u(i,j) = max(cs%Khth\_Min, khth\_loc\_u(i,j))
268 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
269 \textcolor{keywordflow}{  endif}
270 \textcolor{comment}{!$OMP do}
271   \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=is-1,ie
272     kh\_u(i,j,1) = min(kh\_u\_cfl(i,j), khth\_loc\_u(i,j))
273 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
274 
275   \textcolor{keywordflow}{if} (khth\_use\_ebt\_struct) \textcolor{keywordflow}{then}
276 \textcolor{comment}{!$OMP do}
277     \textcolor{keywordflow}{do} k=2,nz+1 ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
278       kh\_u(i,j,k) = kh\_u(i,j,1) * 0.5 * ( varmix%ebt\_struct(i,j,k-1) + varmix%ebt\_struct(i+1,j,k-1) )
279 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
280   \textcolor{keywordflow}{else}
281 \textcolor{comment}{!$OMP do}
282     \textcolor{keywordflow}{do} k=2,nz+1 ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
283       kh\_u(i,j,k) = kh\_u(i,j,1)
284 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
285 \textcolor{keywordflow}{  endif}
286 
287   \textcolor{keywordflow}{if} (use\_varmix) \textcolor{keywordflow}{then}
288     \textcolor{keywordflow}{if} (use\_qg\_leith) \textcolor{keywordflow}{then}
289 \textcolor{comment}{!$OMP do}
290       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
291         kh\_u(i,j,k) = varmix%KH\_u\_QG(i,j,k)
292 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
293 \textcolor{keywordflow}{    endif}
294 \textcolor{keywordflow}{  endif}
295 
296   \textcolor{keywordflow}{if} (cs%use\_GME\_thickness\_diffuse) \textcolor{keywordflow}{then}
297 \textcolor{comment}{!$OMP do}
298     \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
299       cs%KH\_u\_GME(i,j,k) = kh\_u(i,j,k)
300 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
301 \textcolor{keywordflow}{  endif}
302 
303 \textcolor{comment}{!$OMP do}
304   \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
305     khth\_loc\_v(i,j) = cs%Khth
306 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
307 
308   \textcolor{keywordflow}{if} (use\_varmix) \textcolor{keywordflow}{then}
309     \textcolor{keywordflow}{if} (use\_visbeck) \textcolor{keywordflow}{then}
310 \textcolor{comment}{!$OMP do}
311       \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
312         khth\_loc\_v(i,j) = khth\_loc\_v(i,j) + cs%KHTH\_Slope\_Cff*varmix%L2v(i,j)*varmix%SN\_v(i,j)
313 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
314 \textcolor{keywordflow}{    endif}
315 \textcolor{keywordflow}{  endif}
316   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke%Kh)) \textcolor{keywordflow}{then}
317     \textcolor{keywordflow}{if} (cs%MEKE\_GEOMETRIC) \textcolor{keywordflow}{then}
318 \textcolor{comment}{!$OMP do}
319       \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
320         khth\_loc\_v(i,j) = khth\_loc\_v(i,j) + g%mask2dCv(i,j) * cs%MEKE\_GEOMETRIC\_alpha * &
321                         0.5*(meke%MEKE(i,j)+meke%MEKE(i,j+1)) / &
322                         (varmix%SN\_v(i,j) + cs%MEKE\_GEOMETRIC\_epsilon)
323 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
324     \textcolor{keywordflow}{else}
325       \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
326         khth\_loc\_v(i,j) = khth\_loc\_v(i,j) + meke%KhTh\_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
327 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
328 \textcolor{keywordflow}{    endif}
329 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
330 
331   \textcolor{keywordflow}{if} (resoln\_scaled) \textcolor{keywordflow}{then}
332 \textcolor{comment}{!$OMP do}
333     \textcolor{keywordflow}{do} j=js-1,je; \textcolor{keywordflow}{do} i=is,ie
334       khth\_loc\_v(i,j) = khth\_loc\_v(i,j) * varmix%Res\_fn\_v(i,j)
335 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
336 \textcolor{keywordflow}{  endif}
337 
338   \textcolor{keywordflow}{if} (depth\_scaled) \textcolor{keywordflow}{then}
339 \textcolor{comment}{!$OMP do}
340     \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
341       khth\_loc\_v(i,j) = khth\_loc\_v(i,j) * varmix%Depth\_fn\_v(i,j)
342 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
343 \textcolor{keywordflow}{  endif}
344 
345   \textcolor{keywordflow}{if} (cs%Khth\_Max > 0) \textcolor{keywordflow}{then}
346 \textcolor{comment}{!$OMP do}
347     \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
348       khth\_loc\_v(i,j) = max(cs%Khth\_Min, min(khth\_loc\_v(i,j), cs%Khth\_Max))
349 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
350   \textcolor{keywordflow}{else}
351 \textcolor{comment}{!$OMP do}
352     \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
353       khth\_loc\_v(i,j) = max(cs%Khth\_Min, khth\_loc\_v(i,j))
354 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
355 \textcolor{keywordflow}{  endif}
356 
357   \textcolor{keywordflow}{if} (cs%max\_Khth\_CFL > 0.0) \textcolor{keywordflow}{then}
358 \textcolor{comment}{!$OMP do}
359     \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
360       kh\_v(i,j,1) = min(kh\_v\_cfl(i,j), khth\_loc\_v(i,j))
361 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
362 \textcolor{keywordflow}{  endif}
363 
364   \textcolor{keywordflow}{if} (khth\_use\_ebt\_struct) \textcolor{keywordflow}{then}
365 \textcolor{comment}{!$OMP do}
366     \textcolor{keywordflow}{do} k=2,nz+1 ; \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
367       kh\_v(i,j,k) = kh\_v(i,j,1) * 0.5 * ( varmix%ebt\_struct(i,j,k-1) + varmix%ebt\_struct(i,j+1,k-1) )
368 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
369   \textcolor{keywordflow}{else}
370 \textcolor{comment}{!$OMP do}
371     \textcolor{keywordflow}{do} k=2,nz+1 ; \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
372       kh\_v(i,j,k) = kh\_v(i,j,1)
373 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
374 \textcolor{keywordflow}{  endif}
375 
376   \textcolor{keywordflow}{if} (use\_varmix) \textcolor{keywordflow}{then}
377     \textcolor{keywordflow}{if} (use\_qg\_leith) \textcolor{keywordflow}{then}
378 \textcolor{comment}{!$OMP do}
379       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
380         kh\_v(i,j,k) = varmix%KH\_v\_QG(i,j,k)
381 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
382 \textcolor{keywordflow}{    endif}
383 \textcolor{keywordflow}{  endif}
384 
385   \textcolor{keywordflow}{if} (cs%use\_GME\_thickness\_diffuse) \textcolor{keywordflow}{then}
386 \textcolor{comment}{!$OMP do}
387     \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
388       cs%KH\_v\_GME(i,j,k) = kh\_v(i,j,k)
389 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
390 \textcolor{keywordflow}{  endif}
391 
392   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke%Kh)) \textcolor{keywordflow}{then}
393     \textcolor{keywordflow}{if} (cs%MEKE\_GEOMETRIC) \textcolor{keywordflow}{then}
394       \textcolor{keywordflow}{if} (cs%MEKE\_GEOM\_answers\_2018) \textcolor{keywordflow}{then}
395         \textcolor{comment}{!$OMP do}
396         \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
397           \textcolor{comment}{! This does not give bitwise rotational symmetry.}
398           meke%Kh(i,j) = cs%MEKE\_GEOMETRIC\_alpha * meke%MEKE(i,j) / &
399                          (0.25*(varmix%SN\_u(i,j)+varmix%SN\_u(i-1,j) + &
400                                 varmix%SN\_v(i,j)+varmix%SN\_v(i,j-1)) + &
401                           cs%MEKE\_GEOMETRIC\_epsilon)
402 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
403       \textcolor{keywordflow}{else}
404         \textcolor{comment}{!$OMP do}
405         \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
406           \textcolor{comment}{! With the additional parentheses this gives bitwise rotational symmetry.}
407           meke%Kh(i,j) = cs%MEKE\_GEOMETRIC\_alpha * meke%MEKE(i,j) / &
408                          (0.25*((varmix%SN\_u(i,j)+varmix%SN\_u(i-1,j)) + &
409                                 (varmix%SN\_v(i,j)+varmix%SN\_v(i,j-1))) + &
410                           cs%MEKE\_GEOMETRIC\_epsilon)
411 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
412 \textcolor{keywordflow}{      endif}
413 \textcolor{keywordflow}{    endif}
414 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
415 
416 
417 \textcolor{comment}{!$OMP do}
418   \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie ; int\_slope\_u(i,j,k) = 0.0 ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
419 \textcolor{comment}{!$OMP do}
420   \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie ; int\_slope\_v(i,j,k) = 0.0 ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
421 \textcolor{comment}{!$OMP end parallel}
422 
423   \textcolor{keywordflow}{if} (cs%detangle\_interfaces) \textcolor{keywordflow}{then}
424     \textcolor{keyword}{call }add\_detangling\_kh(h, e, kh\_u, kh\_v, kh\_u\_cfl, kh\_v\_cfl, tv, dt, g, gv, us, &
425                            cs, int\_slope\_u, int\_slope\_v)
426 \textcolor{keywordflow}{  endif}
427 
428   \textcolor{keywordflow}{if} (cs%debug) \textcolor{keywordflow}{then}
429     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"Kh\_[uv]"}, kh\_u, kh\_v, g%HI, haloshift=0, &
430                   scale=(us%L\_to\_m**2)*us%s\_to\_T, scalar\_pair=.true.)
431     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"int\_slope\_[uv]"}, int\_slope\_u, int\_slope\_v, g%HI, haloshift=0)
432     \textcolor{keyword}{call }hchksum(h, \textcolor{stringliteral}{"thickness\_diffuse\_1 h"}, g%HI, haloshift=1, scale=gv%H\_to\_m)
433     \textcolor{keyword}{call }hchksum(e, \textcolor{stringliteral}{"thickness\_diffuse\_1 e"}, g%HI, haloshift=1, scale=us%Z\_to\_m)
434     \textcolor{keywordflow}{if} (use\_stored\_slopes) \textcolor{keywordflow}{then}
435       \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"VarMix%slope\_[xy]"}, varmix%slope\_x, varmix%slope\_y, &
436                     g%HI, haloshift=0)
437 \textcolor{keywordflow}{    endif}
438     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%eqn\_of\_state)) \textcolor{keywordflow}{then}
439       \textcolor{keyword}{call }hchksum(tv%T, \textcolor{stringliteral}{"thickness\_diffuse T"}, g%HI, haloshift=1)
440       \textcolor{keyword}{call }hchksum(tv%S, \textcolor{stringliteral}{"thickness\_diffuse S"}, g%HI, haloshift=1)
441 \textcolor{keywordflow}{    endif}
442 \textcolor{keywordflow}{  endif}
443 
444   \textcolor{comment}{! Calculate uhD, vhD from h, e, KH\_u, KH\_v, tv%T/S}
445   \textcolor{keywordflow}{if} (use\_stored\_slopes) \textcolor{keywordflow}{then}
446     \textcolor{keyword}{call }thickness\_diffuse\_full(h, e, kh\_u, kh\_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
447                                 int\_slope\_u, int\_slope\_v, varmix%slope\_x, varmix%slope\_y)
448   \textcolor{keywordflow}{else}
449     \textcolor{keyword}{call }thickness\_diffuse\_full(h, e, kh\_u, kh\_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
450                                 int\_slope\_u, int\_slope\_v)
451 \textcolor{keywordflow}{  endif}
452 
453   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke) .AND. \textcolor{keyword}{associated}(varmix)) \textcolor{keywordflow}{then}
454     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke%Rd\_dx\_h) .and. \textcolor{keyword}{associated}(varmix%Rd\_dx\_h)) \textcolor{keywordflow}{then}
455 \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,MEKE,VarMix)}
456       \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
457         meke%Rd\_dx\_h(i,j) = varmix%Rd\_dx\_h(i,j)
458 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
459 \textcolor{keywordflow}{    endif}
460 \textcolor{keywordflow}{  endif}
461 
462   \textcolor{comment}{! offer diagnostic fields for averaging}
463   \textcolor{keywordflow}{if} (query\_averaging\_enabled(cs%diag)) \textcolor{keywordflow}{then}
464     \textcolor{keywordflow}{if} (cs%id\_uhGM > 0)   \textcolor{keyword}{call }post\_data(cs%id\_uhGM, uhd, cs%diag)
465     \textcolor{keywordflow}{if} (cs%id\_vhGM > 0)   \textcolor{keyword}{call }post\_data(cs%id\_vhGM, vhd, cs%diag)
466     \textcolor{keywordflow}{if} (cs%id\_GMwork > 0) \textcolor{keyword}{call }post\_data(cs%id\_GMwork, cs%GMwork, cs%diag)
467     \textcolor{keywordflow}{if} (cs%id\_KH\_u > 0)   \textcolor{keyword}{call }post\_data(cs%id\_KH\_u, kh\_u, cs%diag)
468     \textcolor{keywordflow}{if} (cs%id\_KH\_v > 0)   \textcolor{keyword}{call }post\_data(cs%id\_KH\_v, kh\_v, cs%diag)
469     \textcolor{keywordflow}{if} (cs%id\_KH\_u1 > 0)  \textcolor{keyword}{call }post\_data(cs%id\_KH\_u1, kh\_u(:,:,1), cs%diag)
470     \textcolor{keywordflow}{if} (cs%id\_KH\_v1 > 0)  \textcolor{keyword}{call }post\_data(cs%id\_KH\_v1, kh\_v(:,:,1), cs%diag)
471 
472     \textcolor{comment}{! Diagnose diffusivity at T-cell point.  Do simple average, rather than}
473     \textcolor{comment}{! thickness-weighted average, in order that KH\_t is depth-independent}
474     \textcolor{comment}{! in the case where KH\_u and KH\_v are depth independent.  Otherwise,}
475     \textcolor{comment}{! if use thickness weighted average, the variations of thickness with}
476     \textcolor{comment}{! depth will place a spurious depth dependence to the diagnosed KH\_t.}
477     \textcolor{keywordflow}{if} (cs%id\_KH\_t > 0 .or. cs%id\_KH\_t1 > 0 .or. cs%Use\_KH\_in\_MEKE) \textcolor{keywordflow}{then}
478       \textcolor{keywordflow}{do} k=1,nz
479         \textcolor{comment}{! thicknesses across u and v faces, converted to 0/1 mask}
480         \textcolor{comment}{! layer average of the interface diffusivities KH\_u and KH\_v}
481         \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
482           hu(i,j)       = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h\_neglect)
483           \textcolor{keywordflow}{if} (hu(i,j) /= 0.0) hu(i,j) = 1.0
484           kh\_u\_lay(i,j) = 0.5*(kh\_u(i,j,k)+kh\_u(i,j,k+1))
485 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
486         \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
487           hv(i,j)       = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h\_neglect)
488           \textcolor{keywordflow}{if} (hv(i,j) /= 0.0) hv(i,j) = 1.0
489           kh\_v\_lay(i,j) = 0.5*(kh\_v(i,j,k)+kh\_v(i,j,k+1))
490 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
491         \textcolor{comment}{! diagnose diffusivity at T-point}
492         \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
493           kh\_t(i,j,k) = ((hu(i-1,j)*kh\_u\_lay(i-1,j)+hu(i,j)*kh\_u\_lay(i,j))  &
494                         +(hv(i,j-1)*kh\_v\_lay(i,j-1)+hv(i,j)*kh\_v\_lay(i,j))) &
495                        / (hu(i-1,j)+hu(i,j)+hv(i,j-1)+hv(i,j)+h\_neglect)
496 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
497 \textcolor{keywordflow}{      enddo}
498 
499       \textcolor{keywordflow}{if} (cs%Use\_KH\_in\_MEKE) \textcolor{keywordflow}{then}
500         meke%Kh\_diff(:,:) = 0.0
501         \textcolor{keywordflow}{do} k=1,nz
502           \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
503             meke%Kh\_diff(i,j) = meke%Kh\_diff(i,j) + kh\_t(i,j,k) * h(i,j,k)
504 \textcolor{keyword}{          end}do; enddo
505 \textcolor{keywordflow}{        enddo}
506 
507         \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
508           meke%Kh\_diff(i,j) = meke%Kh\_diff(i,j) / max(1.0,g%bathyT(i,j))
509 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
510 \textcolor{keywordflow}{      endif}
511 
512       \textcolor{keywordflow}{if} (cs%id\_KH\_t  > 0) \textcolor{keyword}{call }post\_data(cs%id\_KH\_t,  kh\_t,        cs%diag)
513       \textcolor{keywordflow}{if} (cs%id\_KH\_t1 > 0) \textcolor{keyword}{call }post\_data(cs%id\_KH\_t1, kh\_t(:,:,1), cs%diag)
514 \textcolor{keywordflow}{    endif}
515 
516 \textcolor{keywordflow}{  endif}
517 
518   \textcolor{comment}{!$OMP parallel do default(none) shared(is,ie,js,je,nz,uhtr,uhD,dt,vhtr,CDp,vhD,h,G,GV)}
519   \textcolor{keywordflow}{do} k=1,nz
520     \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
521       uhtr(i,j,k) = uhtr(i,j,k) + uhd(i,j,k) * dt
522       \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cdp%uhGM)) cdp%uhGM(i,j,k) = uhd(i,j,k)
523 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
524     \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
525       vhtr(i,j,k) = vhtr(i,j,k) + vhd(i,j,k) * dt
526       \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cdp%vhGM)) cdp%vhGM(i,j,k) = vhd(i,j,k)
527 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
528     \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
529       h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * &
530           ((uhd(i,j,k) - uhd(i-1,j,k)) + (vhd(i,j,k) - vhd(i,j-1,k)))
531       \textcolor{keywordflow}{if} (h(i,j,k) < gv%Angstrom\_H) h(i,j,k) = gv%Angstrom\_H
532 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
533 \textcolor{keywordflow}{  enddo}
534 
535   \textcolor{comment}{! Whenever thickness changes let the diag manager know, target grids}
536   \textcolor{comment}{! for vertical remapping may need to be regenerated.}
537   \textcolor{comment}{! This needs to happen after the H update and before the next post\_data.}
538   \textcolor{keyword}{call }diag\_update\_remap\_grids(cs%diag)
539 
540   \textcolor{keywordflow}{if} (cs%debug) \textcolor{keywordflow}{then}
541     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"thickness\_diffuse [uv]hD"}, uhd, vhd, &
542                   g%HI, haloshift=0, scale=gv%H\_to\_m*us%L\_to\_m**2*us%s\_to\_T)
543     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"thickness\_diffuse [uv]htr"}, uhtr, vhtr, &
544                   g%HI, haloshift=0, scale=us%L\_to\_m**2*gv%H\_to\_m)
545     \textcolor{keyword}{call }hchksum(h, \textcolor{stringliteral}{"thickness\_diffuse h"}, g%HI, haloshift=0, scale=gv%H\_to\_m)
546 \textcolor{keywordflow}{  endif}
547 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__thickness__diffuse_a1a3a5c750e8a2323afa82b118329378c}\label{namespacemom__thickness__diffuse_a1a3a5c750e8a2323afa82b118329378c}} 
\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}!thickness\+\_\+diffuse\+\_\+end@{thickness\+\_\+diffuse\+\_\+end}}
\index{thickness\+\_\+diffuse\+\_\+end@{thickness\+\_\+diffuse\+\_\+end}!mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}
\subsubsection{\texorpdfstring{thickness\+\_\+diffuse\+\_\+end()}{thickness\_diffuse\_end()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+thickness\+\_\+diffuse\+::thickness\+\_\+diffuse\+\_\+end (\begin{DoxyParamCaption}\item[{type(\hyperlink{structmom__thickness__diffuse_1_1thickness__diffuse__cs}{thickness\+\_\+diffuse\+\_\+cs}), pointer}]{CS }\end{DoxyParamCaption})}



Deallocate the thickness diffusion control structure. 


\begin{DoxyParams}{Parameters}
{\em cs} & Control structure for thickness diffusion \\
\hline
\end{DoxyParams}


Definition at line 2113 of file M\+O\+M\+\_\+thickness\+\_\+diffuse.\+F90.


\begin{DoxyCode}
2113   \textcolor{keywordtype}{type}(thickness\_diffuse\_cs), \textcolor{keywordtype}{pointer} :: cs\textcolor{comment}{ !< Control structure for thickness diffusion}
2114 
2115   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{deallocate}(cs)
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__thickness__diffuse_ae9909642254fcf0160afe46997e10c30}\label{namespacemom__thickness__diffuse_ae9909642254fcf0160afe46997e10c30}} 
\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}!thickness\+\_\+diffuse\+\_\+full@{thickness\+\_\+diffuse\+\_\+full}}
\index{thickness\+\_\+diffuse\+\_\+full@{thickness\+\_\+diffuse\+\_\+full}!mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}
\subsubsection{\texorpdfstring{thickness\+\_\+diffuse\+\_\+full()}{thickness\_diffuse\_full()}}
{\footnotesize\ttfamily subroutine mom\+\_\+thickness\+\_\+diffuse\+::thickness\+\_\+diffuse\+\_\+full (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(in)}]{e,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(in)}]{Kh\+\_\+u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)+1), intent(in)}]{Kh\+\_\+v,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(out)}]{uhD,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)), intent(out)}]{vhD,  }\item[{real, dimension(\+:,\+:), pointer}]{cg1,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(meke\+\_\+type), pointer}]{M\+E\+KE,  }\item[{type(\hyperlink{structmom__thickness__diffuse_1_1thickness__diffuse__cs}{thickness\+\_\+diffuse\+\_\+cs}), pointer}]{CS,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(in), optional}]{int\+\_\+slope\+\_\+u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)+1), intent(in), optional}]{int\+\_\+slope\+\_\+v,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(in), optional}]{slope\+\_\+x,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)+1), intent(in), optional}]{slope\+\_\+y }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Calculates parameterized layer transports for use in the continuity equation. Fluxes are limited to give positive definite thicknesses. Called by \hyperlink{namespacemom__thickness__diffuse_a8a538b778a567f489bfd9c5eadeeebef}{thickness\+\_\+diffuse()}. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em e} & Interface positions \mbox{[}Z $\sim$$>$ m\mbox{]}\\
\hline
\mbox{\tt in}  & {\em kh\+\_\+u} & Thickness diffusivity on interfaces at u points \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em kh\+\_\+v} & Thickness diffusivity on interfaces at v points \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em tv} & Thermodynamics structure\\
\hline
\mbox{\tt out}  & {\em uhd} & Zonal mass fluxes \mbox{[}H L2 T-\/1 $\sim$$>$ m3 s-\/1 or kg s-\/1\mbox{]}\\
\hline
\mbox{\tt out}  & {\em vhd} & Meridional mass fluxes \mbox{[}H L2 T-\/1 $\sim$$>$ m3 s-\/1 or kg s-\/1\mbox{]}\\
\hline
 & {\em cg1} & Wave speed \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
 & {\em meke} & M\+E\+KE control structure\\
\hline
 & {\em cs} & Control structure for thickness diffusion\\
\hline
\mbox{\tt in}  & {\em int\+\_\+slope\+\_\+u} & Ratio that determine how much of the isopycnal slopes are taken directly from the interface slopes without consideration of density gradients \mbox{[}nondim\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em int\+\_\+slope\+\_\+v} & Ratio that determine how much of the isopycnal slopes are taken directly from the interface slopes without consideration of density gradients \mbox{[}nondim\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em slope\+\_\+x} & Isopycnal slope at u-\/points\\
\hline
\mbox{\tt in}  & {\em slope\+\_\+y} & Isopycnal slope at v-\/points \\
\hline
\end{DoxyParams}


Definition at line 555 of file M\+O\+M\+\_\+thickness\+\_\+diffuse.\+F90.


\begin{DoxyCode}
555   \textcolor{keywordtype}{type}(ocean\_grid\_type),                       \textcolor{keywordtype}{intent(in)}  :: g\textcolor{comment}{      !< Ocean grid structure}
556   \textcolor{keywordtype}{type}(verticalgrid\_type),                     \textcolor{keywordtype}{intent(in)}  :: gv\textcolor{comment}{     !< Vertical grid structure}
557   \textcolor{keywordtype}{type}(unit\_scale\_type),                       \textcolor{keywordtype}{intent(in)}  :: us\textcolor{comment}{     !< A dimensional unit scaling type}
558   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))},    \textcolor{keywordtype}{intent(in)}  :: h\textcolor{comment}{      !< Layer thickness [H ~> m or kg m-2]}
559   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G)+1)},  \textcolor{keywordtype}{intent(in)}  :: e\textcolor{comment}{      !< Interface positions [Z ~> m]}
560   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(in)}  :: kh\_u\textcolor{comment}{   !< Thickness diffusivity on interfaces}
561 \textcolor{comment}{                                                                     !! at u points [L2 T-1 ~> m2 s-1]}
562   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(in)}  :: kh\_v\textcolor{comment}{   !< Thickness diffusivity on interfaces}
563 \textcolor{comment}{                                                                     !! at v points [L2 T-1 ~> m2 s-1]}
564   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),                       \textcolor{keywordtype}{intent(in)}  :: tv\textcolor{comment}{     !< Thermodynamics structure}
565   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))},   \textcolor{keywordtype}{intent(out)} :: uhd\textcolor{comment}{    !< Zonal mass fluxes}
566 \textcolor{comment}{                                                                     !! [H L2 T-1 ~> m3 s-1 or kg s-1]}
567   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))},   \textcolor{keywordtype}{intent(out)} :: vhd\textcolor{comment}{    !< Meridional mass fluxes}
568 \textcolor{comment}{                                                                     !! [H L2 T-1 ~> m3 s-1 or kg s-1]}
569   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)},                        \textcolor{keywordtype}{pointer}     :: cg1\textcolor{comment}{    !< Wave speed [L T-1 ~> m s-1]}
570   \textcolor{keywordtype}{real},                                        \textcolor{keywordtype}{intent(in)}  :: dt\textcolor{comment}{     !< Time increment [T ~> s]}
571   \textcolor{keywordtype}{type}(meke\_type),                             \textcolor{keywordtype}{pointer}     :: meke\textcolor{comment}{   !< MEKE control structure}
572   \textcolor{keywordtype}{type}(thickness\_diffuse\_cs),                  \textcolor{keywordtype}{pointer}     :: cs\textcolor{comment}{     !< Control structure for thickness
       diffusion}
573   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: int\_slope\_u\textcolor{comment}{ !< Ratio that determine
       how much of}
574 \textcolor{comment}{                                                                     !! the isopycnal slopes are taken
       directly from}
575 \textcolor{comment}{                                                                     !! the interface slopes without
       consideration of}
576 \textcolor{comment}{                                                                     !! density gradients [nondim].}
577   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: int\_slope\_v\textcolor{comment}{ !< Ratio that determine
       how much of}
578 \textcolor{comment}{                                                                     !! the isopycnal slopes are taken
       directly from}
579 \textcolor{comment}{                                                                     !! the interface slopes without
       consideration of}
580 \textcolor{comment}{                                                                     !! density gradients [nondim].}
581   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: slope\_x\textcolor{comment}{ !< Isopycnal slope at
       u-points}
582   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: slope\_y\textcolor{comment}{ !< Isopycnal slope at
       v-points}
583   \textcolor{comment}{! Local variables}
584   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G), SZK\_(G))} :: &
585     t, &          \textcolor{comment}{! The temperature (or density) [degC], with the values in}
586                   \textcolor{comment}{! in massless layers filled vertically by diffusion.}
587     s, &          \textcolor{comment}{! The filled salinity [ppt], with the values in}
588                   \textcolor{comment}{! in massless layers filled vertically by diffusion.}
589     h\_avail, &    \textcolor{comment}{! The mass available for diffusion out of each face, divided}
590                   \textcolor{comment}{! by dt [H L2 T-1 ~> m3 s-1 or kg s-1].}
591     h\_frac        \textcolor{comment}{! The fraction of the mass in the column above the bottom}
592                   \textcolor{comment}{! interface of a layer that is within a layer [nondim]. 0<h\_frac<=1}
593   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJB\_(G), SZK\_(G)+1)} :: &
594     slope\_y\_pe, &  \textcolor{comment}{! 3D array of neutral slopes at v-points, set equal to Slope (below, nondim)}
595     hn2\_y\_pe       \textcolor{comment}{! thickness in m times Brunt-Vaisala freqeuncy at v-points [L2 Z-1 T-2 ~> m s-2],}
596                    \textcolor{comment}{! used for calculating PE release}
597   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G), SZJ\_(G), SZK\_(G)+1)} :: &
598     slope\_x\_pe, &  \textcolor{comment}{! 3D array of neutral slopes at u-points, set equal to Slope (below, nondim)}
599     hn2\_x\_pe       \textcolor{comment}{! thickness in m times Brunt-Vaisala freqeuncy at u-points [L2 Z-1 T-2 ~> m s-2],}
600                    \textcolor{comment}{! used for calculating PE release}
601   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G), SZK\_(G)+1)} :: &
602     pres, &       \textcolor{comment}{! The pressure at an interface [R L2 T-2 ~> Pa].}
603     h\_avail\_rsum  \textcolor{comment}{! The running sum of h\_avail above an interface [H L2 T-1 ~> m3 s-1 or kg s-1].}
604   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
605     drho\_dt\_u, &  \textcolor{comment}{! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1]}
606     drho\_ds\_u, &  \textcolor{comment}{! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1].}
607     drho\_dt\_dt\_u  \textcolor{comment}{! The second derivative of density with temperature at u points [R degC-2 ~> kg m-3
       degC-2]}
608   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: scrap \textcolor{comment}{! An array to pass to calculate\_density\_second\_derivs() that will be
       ingored.}
609   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
610     drho\_dt\_v, &  \textcolor{comment}{! The derivative of density with temperature at v points [R degC-1 ~> kg m-3 degC-1]}
611     drho\_ds\_v, &  \textcolor{comment}{! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1].}
612     drho\_dt\_dt\_v  \textcolor{comment}{! The second derivative of density with temperature at v points [R degC-2 ~> kg m-3
       degC-2]}
613   \textcolor{keywordtype}{real} :: uhtot(szib\_(g), szj\_(g))  \textcolor{comment}{! The vertical sum of uhD [H L2 T-1 ~> m3 s-1 or kg s-1].}
614   \textcolor{keywordtype}{real} :: vhtot(szi\_(g), szjb\_(g))  \textcolor{comment}{! The vertical sum of vhD [H L2 T-1 ~> m3 s-1 or kg s-1].}
615   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
616     t\_u, &        \textcolor{comment}{! Temperature on the interface at the u-point [degC].}
617     s\_u, &        \textcolor{comment}{! Salinity on the interface at the u-point [ppt].}
618     pres\_u        \textcolor{comment}{! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].}
619   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G))} :: &
620     t\_v, &        \textcolor{comment}{! Temperature on the interface at the v-point [degC].}
621     s\_v, &        \textcolor{comment}{! Salinity on the interface at the v-point [ppt].}
622     pres\_v        \textcolor{comment}{! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].}
623   \textcolor{keywordtype}{real} :: work\_u(szib\_(g), szj\_(g)) \textcolor{comment}{! The work being done by the thickness}
624   \textcolor{keywordtype}{real} :: work\_v(szi\_(g), szjb\_(g)) \textcolor{comment}{! diffusion integrated over a cell [R Z L4 T-3  ~> W ]}
625   \textcolor{keywordtype}{real} :: work\_h        \textcolor{comment}{! The work averaged over an h-cell [R Z L2 T-3 ~> W m-2].}
626   \textcolor{keywordtype}{real} :: pe\_release\_h  \textcolor{comment}{! The amount of potential energy released by GM averaged over an h-cell [L4 Z-1 T-3
       ~> m3 s-3]}
627                         \textcolor{comment}{! The calculation is equal to h * S^2 * N^2 * kappa\_GM.}
628   \textcolor{keywordtype}{real} :: i4dt          \textcolor{comment}{! 1 / 4 dt [T-1 ~> s-1].}
629   \textcolor{keywordtype}{real} :: drdia, drdib  \textcolor{comment}{! Along layer zonal- and meridional- potential density}
630   \textcolor{keywordtype}{real} :: drdja, drdjb  \textcolor{comment}{! gradients in the layers above (A) and below(B) the}
631                         \textcolor{comment}{! interface times the grid spacing [R ~> kg m-3].}
632   \textcolor{keywordtype}{real} :: drdkl, drdkr  \textcolor{comment}{! Vertical density differences across an interface [R ~> kg m-3].}
633   \textcolor{keywordtype}{real} :: drdi\_u(szib\_(g), szk\_(g)) \textcolor{comment}{! Copy of drdi at u-points [R ~> kg m-3].}
634   \textcolor{keywordtype}{real} :: drdj\_v(szi\_(g), szk\_(g))  \textcolor{comment}{! Copy of drdj at v-points [R ~> kg m-3].}
635   \textcolor{keywordtype}{real} :: drdkde\_u(szib\_(g),szk\_(g)+1) \textcolor{comment}{! Lateral difference of product of drdk and e at u-points}
636                                        \textcolor{comment}{! [Z R ~> kg m-2].}
637   \textcolor{keywordtype}{real} :: drdkde\_v(szi\_(g),szk\_(g)+1)  \textcolor{comment}{! Lateral difference of product of drdk and e at v-points}
638                                        \textcolor{comment}{! [Z R ~> kg m-2].}
639   \textcolor{keywordtype}{real} :: hg2a, hg2b, hg2l, hg2r \textcolor{comment}{! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].}
640   \textcolor{keywordtype}{real} :: haa, hab, hal, har     \textcolor{comment}{! Arithmetic mean thicknesses [H ~> m or kg m-2].}
641   \textcolor{keywordtype}{real} :: dzal, dzar    \textcolor{comment}{! Temporary thicknesses [Z ~> m].}
642   \textcolor{keywordtype}{real} :: wta, wtb, wtl, wtr  \textcolor{comment}{! Unscaled weights, with various units.}
643   \textcolor{keywordtype}{real} :: drdx, drdy    \textcolor{comment}{! Zonal and meridional density gradients [R L-1 ~> kg m-4].}
644   \textcolor{keywordtype}{real} :: drdz          \textcolor{comment}{! Vertical density gradient [R Z-1 ~> kg m-4].}
645   \textcolor{keywordtype}{real} :: h\_harm        \textcolor{comment}{! Harmonic mean layer thickness [H ~> m or kg m-2].}
646   \textcolor{keywordtype}{real} :: c2\_h\_u(szib\_(g), szk\_(g)+1) \textcolor{comment}{! Wave speed squared divided by h at u-points [L2 Z-1 T-2 ~> m s-2].}
647   \textcolor{keywordtype}{real} :: c2\_h\_v(szi\_(g), szk\_(g)+1)  \textcolor{comment}{! Wave speed squared divided by h at v-points [L2 Z-1 T-2 ~> m s-2].}
648   \textcolor{keywordtype}{real} :: hn2\_u(szib\_(g), szk\_(g)+1)  \textcolor{comment}{! Thickness in m times N2 at interfaces above u-points [L2 Z-1 T-2 ~>
       m s-2].}
649   \textcolor{keywordtype}{real} :: hn2\_v(szi\_(g), szk\_(g)+1)   \textcolor{comment}{! Thickness in m times N2 at interfaces above v-points [L2 Z-1 T-2 ~>
       m s-2].}
650   \textcolor{keywordtype}{real} :: sfn\_est       \textcolor{comment}{! A preliminary estimate (before limiting) of the overturning}
651                         \textcolor{comment}{! streamfunction [Z L2 T-1 ~> m3 s-1].}
652   \textcolor{keywordtype}{real} :: sfn\_unlim\_u(szib\_(g), szk\_(g)+1) \textcolor{comment}{! Streamfunction for u-points [Z L2 T-1 ~> m3 s-1].}
653   \textcolor{keywordtype}{real} :: sfn\_unlim\_v(szi\_(g), szk\_(g)+1)  \textcolor{comment}{! Streamfunction for v-points [Z L2 T-1 ~> m3 s-1].}
654   \textcolor{keywordtype}{real} :: slope2\_ratio\_u(szib\_(g), szk\_(g)+1) \textcolor{comment}{! The ratio of the slope squared to slope\_max squared.}
655   \textcolor{keywordtype}{real} :: slope2\_ratio\_v(szi\_(g), szk\_(g)+1)  \textcolor{comment}{! The ratio of the slope squared to slope\_max squared.}
656   \textcolor{keywordtype}{real} :: sfn\_in\_h      \textcolor{comment}{! The overturning streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1] (note that}
657                         \textcolor{comment}{! the units are different from other Sfn vars).}
658   \textcolor{keywordtype}{real} :: sfn\_safe      \textcolor{comment}{! The streamfunction that goes linearly back to 0 at the surface.  This is a}
659                         \textcolor{comment}{! good thing to use when the slope is so large as to be meaningless [Z L2 T-1 ~> m3
       s-1].}
660   \textcolor{keywordtype}{real} :: slope         \textcolor{comment}{! The slope of density surfaces, calculated in a way}
661                         \textcolor{comment}{! that is always between -1 and 1, nondimensional.}
662   \textcolor{keywordtype}{real} :: mag\_grad2     \textcolor{comment}{! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8].}
663   \textcolor{keywordtype}{real} :: i\_slope\_max2  \textcolor{comment}{! The inverse of slope\_max squared, nondimensional.}
664   \textcolor{keywordtype}{real} :: h\_neglect     \textcolor{comment}{! A thickness that is so small it is usually lost}
665                         \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
666   \textcolor{keywordtype}{real} :: h\_neglect2    \textcolor{comment}{! h\_neglect^2 [H2 ~> m2 or kg2 m-4].}
667   \textcolor{keywordtype}{real} :: dz\_neglect    \textcolor{comment}{! A thickness [Z ~> m], that is so small it is usually lost}
668                         \textcolor{comment}{! in roundoff and can be neglected [Z ~> m].}
669   \textcolor{keywordtype}{real} :: g\_scale       \textcolor{comment}{! The gravitational acceleration times a unit conversion}
670                         \textcolor{comment}{! factor [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].}
671   \textcolor{keywordtype}{logical} :: use\_eos    \textcolor{comment}{! If true, density is calculated from T & S using an}
672                         \textcolor{comment}{! equation of state.}
673   \textcolor{keywordtype}{logical} :: find\_work  \textcolor{comment}{! If true, find the change in energy due to the fluxes.}
674   \textcolor{keywordtype}{integer} :: nk\_linear  \textcolor{comment}{! The number of layers over which the streamfunction goes to 0.}
675   \textcolor{keywordtype}{real} :: g\_rho0        \textcolor{comment}{! g/Rho0 [L2 R-1 Z-1 T-2 ~> m4 kg-1 s-2].}
676   \textcolor{keywordtype}{real} :: n2\_floor      \textcolor{comment}{! A floor for N2 to avoid degeneracy in the elliptic solver}
677                         \textcolor{comment}{! times unit conversion factors [T-2 L2 Z-2 ~> s-2]}
678   \textcolor{keywordtype}{real} :: tl(5)         \textcolor{comment}{! copy and T in local stencil [degC]}
679   \textcolor{keywordtype}{real} :: mn\_t          \textcolor{comment}{! mean of T in local stencil [degC]}
680   \textcolor{keywordtype}{real} :: mn\_t2         \textcolor{comment}{! mean of T**2 in local stencil [degC]}
681   \textcolor{keywordtype}{real} :: hl(5)         \textcolor{comment}{! Copy of local stencil of H [H ~> m]}
682   \textcolor{keywordtype}{real} :: r\_sm\_h        \textcolor{comment}{! Reciprocal of sum of H in local stencil [H-1 ~> m-1]}
683   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJ\_(G), SZK\_(G))} :: tsgs2 \textcolor{comment}{! Sub-grid temperature variance [degC2]}
684 
685   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G), SZJ\_(G), SZK\_(G)+1)} :: diag\_sfn\_x, diag\_sfn\_unlim\_x \textcolor{comment}{! Diagnostics}
686   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G), SZJB\_(G), SZK\_(G)+1)} :: diag\_sfn\_y, diag\_sfn\_unlim\_y \textcolor{comment}{! Diagnostics}
687   \textcolor{keywordtype}{logical} :: present\_int\_slope\_u, present\_int\_slope\_v
688   \textcolor{keywordtype}{logical} :: present\_slope\_x, present\_slope\_y, calc\_derivatives
689   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(2)} ::  eosdom\_u \textcolor{comment}{! The shifted i-computational domain to use for equation of}
690                                      \textcolor{comment}{! state calculations at u-points.}
691   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(2)} ::  eosdom\_v \textcolor{comment}{! The shifted I-computational domain to use for equation of}
692                                      \textcolor{comment}{! state calculations at v-points.}
693   \textcolor{keywordtype}{logical} :: use\_stanley
694   \textcolor{keywordtype}{integer} :: is, ie, js, je, nz, isdb, halo
695   \textcolor{keywordtype}{integer} :: i, j, k
696   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke ; isdb = g%IsdB
697 
698   i4dt = 0.25 / dt
699   i\_slope\_max2 = 1.0 / (cs%slope\_max**2)
700   g\_scale = gv%g\_Earth * gv%H\_to\_Z
701 
702   h\_neglect = gv%H\_subroundoff ; h\_neglect2 = h\_neglect**2
703   dz\_neglect = gv%H\_subroundoff*gv%H\_to\_Z
704   g\_rho0 = gv%g\_Earth / gv%Rho0
705   n2\_floor = cs%N2\_floor*us%Z\_to\_L**2
706 
707   use\_eos = \textcolor{keyword}{associated}(tv%eqn\_of\_state)
708   present\_int\_slope\_u = \textcolor{keyword}{PRESENT}(int\_slope\_u)
709   present\_int\_slope\_v = \textcolor{keyword}{PRESENT}(int\_slope\_v)
710   present\_slope\_x = \textcolor{keyword}{PRESENT}(slope\_x)
711   present\_slope\_y = \textcolor{keyword}{PRESENT}(slope\_y)
712   use\_stanley = cs%Stanley\_det\_coeff >= 0.
713 
714   nk\_linear = max(gv%nkml, 1)
715 
716   slope\_x\_pe(:,:,:) = 0.0
717   slope\_y\_pe(:,:,:) = 0.0
718   hn2\_x\_pe(:,:,:) = 0.0
719   hn2\_y\_pe(:,:,:) = 0.0
720 
721   find\_work = .false.
722   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke)) find\_work = \textcolor{keyword}{associated}(meke%GM\_src)
723   find\_work = (\textcolor{keyword}{associated}(cs%GMwork) .or. find\_work)
724 
725   \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
726     halo = 1 \textcolor{comment}{! Default halo to fill is 1}
727     \textcolor{keywordflow}{if} (use\_stanley) halo = 2 \textcolor{comment}{! Need wider valid halo for gradients of T}
728     \textcolor{keyword}{call }vert\_fill\_ts(h, tv%T, tv%S, cs%kappa\_smooth*dt, t, s, g, gv, halo, larger\_h\_denom=.true.)
729 \textcolor{keywordflow}{  endif}
730 
731   \textcolor{keywordflow}{if} (cs%use\_FGNV\_streamfn .and. .not. \textcolor{keyword}{associated}(cg1)) \textcolor{keyword}{call }mom\_error(fatal, &
732        \textcolor{stringliteral}{"cg1 must be associated when using FGNV streamfunction."})
733 
734 \textcolor{comment}{!$OMP parallel default(none) shared(is,ie,js,je,h\_avail\_rsum,pres,h\_avail,I4dt, use\_Stanley, &}
735 \textcolor{comment}{!$OMP                               CS,G,GV,tv,h,h\_frac,nz,uhtot,Work\_u,vhtot,Work\_v,Tsgs2,T, &}
736 \textcolor{comment}{!$OMP                               diag\_sfn\_x, diag\_sfn\_y, diag\_sfn\_unlim\_x, diag\_sfn\_unlim\_y ) &}
737 \textcolor{comment}{!$OMP          private(hl,r\_sm\_H,Tl,mn\_T,mn\_T2)}
738   \textcolor{comment}{! Find the maximum and minimum permitted streamfunction.}
739 \textcolor{comment}{!$OMP do}
740   \textcolor{keywordflow}{do} j=js-1,je+1 ; \textcolor{keywordflow}{do} i=is-1,ie+1
741     h\_avail\_rsum(i,j,1) = 0.0
742     pres(i,j,1) = 0.0
743     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%p\_surf)) \textcolor{keywordflow}{then} ; pres(i,j,1) = tv%p\_surf(i,j) ;\textcolor{keywordflow}{ endif}
744 
745     h\_avail(i,j,1) = max(i4dt*g%areaT(i,j)*(h(i,j,1)-gv%Angstrom\_H),0.0)
746     h\_avail\_rsum(i,j,2) = h\_avail(i,j,1)
747     h\_frac(i,j,1) = 1.0
748     pres(i,j,2) = pres(i,j,1) + (gv%g\_Earth*gv%H\_to\_RZ) * h(i,j,1)
749 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
750   \textcolor{keywordflow}{if} (use\_stanley) \textcolor{keywordflow}{then}
751 \textcolor{comment}{!$OMP do}
752     \textcolor{keywordflow}{do} k=1, nz ; \textcolor{keywordflow}{do} j=js-1,je+1 ; \textcolor{keywordflow}{do} i=is-1,ie+1
753       \textcolor{comment}{!! SGS variance in i-direction [degC2]}
754       \textcolor{comment}{!dTdi2 = ( ( G%mask2dCu(I  ,j) * G%IdxCu(I  ,j) * ( T(i+1,j,k) - T(i,j,k) ) &}
755       \textcolor{comment}{!          + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) ) &}
756       \textcolor{comment}{!          ) * G%dxT(i,j) * 0.5 )**2}
757       \textcolor{comment}{!! SGS variance in j-direction [degC2]}
758       \textcolor{comment}{!dTdj2 = ( ( G%mask2dCv(i,J  ) * G%IdyCv(i,J  ) * ( T(i,j+1,k) - T(i,j,k) ) &}
759       \textcolor{comment}{!          + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) ) &}
760       \textcolor{comment}{!          ) * G%dyT(i,j) * 0.5 )**2}
761       \textcolor{comment}{!Tsgs2(i,j,k) = CS%Stanley\_det\_coeff * 0.5 * ( dTdi2 + dTdj2 )}
762       \textcolor{comment}{! This block does a thickness weighted variance calculation and helps control for}
763       \textcolor{comment}{! extreme gradients along layers which are vanished against topography. It is}
764       \textcolor{comment}{! still a poor approximation in the interior when coordinates are strongly tilted.}
765       hl(1) = h(i,j,k) * g%mask2dT(i,j)
766       hl(2) = h(i-1,j,k) * g%mask2dCu(i-1,j)
767       hl(3) = h(i+1,j,k) * g%mask2dCu(i,j)
768       hl(4) = h(i,j-1,k) * g%mask2dCv(i,j-1)
769       hl(5) = h(i,j+1,k) * g%mask2dCv(i,j)
770       r\_sm\_h = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + gv%H\_subroundoff )
771       \textcolor{comment}{! Mean of T}
772       tl(1) = t(i,j,k) ; tl(2) = t(i-1,j,k) ; tl(3) = t(i+1,j,k)
773       tl(4) = t(i,j-1,k) ; tl(5) = t(i,j+1,k)
774       mn\_t = ( hl(1)*tl(1) + ( ( hl(2)*tl(2) + hl(3)*tl(3) ) + ( hl(4)*tl(4) + hl(5)*tl(5) ) ) ) * r\_sm\_h
775       \textcolor{comment}{! Adjust T vectors to have zero mean}
776       tl(:) = tl(:) - mn\_t ; mn\_t = 0.
777       \textcolor{comment}{! Variance of T}
778       mn\_t2 = ( hl(1)*tl(1)*tl(1) + ( ( hl(2)*tl(2)*tl(2) + hl(3)*tl(3)*tl(3) ) &
779                                     + ( hl(4)*tl(4)*tl(4) + hl(5)*tl(5)*tl(5) ) ) ) * r\_sm\_h
780       \textcolor{comment}{! Variance should be positive but round-off can violate this. Calculating}
781       \textcolor{comment}{! variance directly would fix this but requires more operations.}
782       tsgs2(i,j,k) = cs%Stanley\_det\_coeff * max(0., mn\_t2)
783 \textcolor{keywordflow}{     enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
784 \textcolor{keywordflow}{  endif}
785 \textcolor{comment}{!$OMP do}
786   \textcolor{keywordflow}{do} j=js-1,je+1
787     \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is-1,ie+1
788       h\_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom\_H),0.0)
789       h\_avail\_rsum(i,j,k+1) = h\_avail\_rsum(i,j,k) + h\_avail(i,j,k)
790       h\_frac(i,j,k) = 0.0 ; \textcolor{keywordflow}{if} (h\_avail(i,j,k) > 0.0) &
791         h\_frac(i,j,k) = h\_avail(i,j,k) / h\_avail\_rsum(i,j,k+1)
792       pres(i,j,k+1) = pres(i,j,k) + (gv%g\_Earth*gv%H\_to\_RZ) * h(i,j,k)
793 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
794 \textcolor{keywordflow}{  enddo}
795 \textcolor{comment}{!$OMP do}
796   \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie
797     uhtot(i,j) = 0.0 ; work\_u(i,j) = 0.0
798     diag\_sfn\_x(i,j,1) = 0.0 ; diag\_sfn\_unlim\_x(i,j,1) = 0.0
799     diag\_sfn\_x(i,j,nz+1) = 0.0 ; diag\_sfn\_unlim\_x(i,j,nz+1) = 0.0
800 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
801 \textcolor{comment}{!$OMP do}
802   \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie
803     vhtot(i,j) = 0.0 ; work\_v(i,j) = 0.0
804     diag\_sfn\_y(i,j,1) = 0.0 ; diag\_sfn\_unlim\_y(i,j,1) = 0.0
805     diag\_sfn\_y(i,j,nz+1) = 0.0 ; diag\_sfn\_unlim\_y(i,j,nz+1) = 0.0
806 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
807 \textcolor{comment}{!$OMP end parallel}
808 
809     eosdom\_u(1) = (is-1) - (g%IsdB-1) ; eosdom\_u(2) = ie - (g%IsdB-1)
810 \textcolor{comment}{!$OMP parallel do default(none) shared(nz,is,ie,js,je,find\_work,use\_EOS,G,GV,US,pres,T,S, &}
811 \textcolor{comment}{!$OMP                                  nk\_linear,IsdB,tv,h,h\_neglect,e,dz\_neglect,  &}
812 \textcolor{comment}{!$OMP                                  I\_slope\_max2,h\_neglect2,present\_int\_slope\_u, &}
813 \textcolor{comment}{!$OMP                                  int\_slope\_u,KH\_u,uhtot,h\_frac,h\_avail\_rsum,  &}
814 \textcolor{comment}{!$OMP                                  uhD,h\_avail,G\_scale,Work\_u,CS,slope\_x,cg1,   &}
815 \textcolor{comment}{!$OMP                                  diag\_sfn\_x, diag\_sfn\_unlim\_x,N2\_floor,EOSdom\_u, &}
816 \textcolor{comment}{!$OMP                                  use\_stanley, Tsgs2,                          &}
817 \textcolor{comment}{!$OMP                                  present\_slope\_x,G\_rho0,Slope\_x\_PE,hN2\_x\_PE)  &}
818 \textcolor{comment}{!$OMP                          private(drdiA,drdiB,drdkL,drdkR,pres\_u,T\_u,S\_u,      &}
819 \textcolor{comment}{!$OMP                                  drho\_dT\_u,drho\_dS\_u,hg2A,hg2B,hg2L,hg2R,haA, &}
820 \textcolor{comment}{!$OMP                                  drho\_dT\_dT\_u,scrap,                          &}
821 \textcolor{comment}{!$OMP                                  haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,  &}
822 \textcolor{comment}{!$OMP                                  drdx,mag\_grad2,Slope,slope2\_Ratio\_u,hN2\_u,   &}
823 \textcolor{comment}{!$OMP                                  Sfn\_unlim\_u,drdi\_u,drdkDe\_u,h\_harm,c2\_h\_u,   &}
824 \textcolor{comment}{!$OMP                                  Sfn\_safe,Sfn\_est,Sfn\_in\_h,calc\_derivatives)}
825   \textcolor{keywordflow}{do} j=js,je
826     \textcolor{keywordflow}{do} i=is-1,ie ; hn2\_u(i,1) = 0. ; hn2\_u(i,nz+1) = 0. ;\textcolor{keywordflow}{ enddo}
827     \textcolor{keywordflow}{do} k=nz,2,-1
828       \textcolor{keywordflow}{if} (find\_work .and. .not.(use\_eos)) \textcolor{keywordflow}{then}
829         drdia = 0.0 ; drdib = 0.0
830         drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
831 \textcolor{keywordflow}{      endif}
832 
833       calc\_derivatives = use\_eos .and. (k >= nk\_linear) .and. &
834          (find\_work .or. .not. present\_slope\_x .or. cs%use\_FGNV\_streamfn .or. use\_stanley)
835 
836       \textcolor{comment}{! Calculate the zonal fluxes and gradients.}
837       \textcolor{keywordflow}{if} (calc\_derivatives) \textcolor{keywordflow}{then}
838         \textcolor{keywordflow}{do} i=is-1,ie
839           pres\_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
840           t\_u(i) = 0.25*((t(i,j,k) + t(i+1,j,k)) + (t(i,j,k-1) + t(i+1,j,k-1)))
841           s\_u(i) = 0.25*((s(i,j,k) + s(i+1,j,k)) + (s(i,j,k-1) + s(i+1,j,k-1)))
842 \textcolor{keywordflow}{        enddo}
843         \textcolor{keyword}{call }calculate\_density\_derivs(t\_u, s\_u, pres\_u, drho\_dt\_u, drho\_ds\_u, &
844                                       tv%eqn\_of\_state, eosdom\_u)
845 \textcolor{keywordflow}{      endif}
846       \textcolor{keywordflow}{if} (use\_stanley) \textcolor{keywordflow}{then}
847         \textcolor{comment}{! The second line below would correspond to arguments}
848         \textcolor{comment}{!            drho\_dS\_dS, drho\_dS\_dT, drho\_dT\_dT, drho\_dS\_dP, drho\_dT\_dP, &}
849         \textcolor{keyword}{call }calculate\_density\_second\_derivs(t\_u, s\_u, pres\_u, &
850                      scrap, scrap, drho\_dt\_dt\_u, scrap, scrap, &
851                      (is-isdb+1)-1, ie-is+2, tv%eqn\_of\_state)
852 \textcolor{keywordflow}{      endif}
853 
854       \textcolor{keywordflow}{do} i=is-1,ie
855         \textcolor{keywordflow}{if} (calc\_derivatives) \textcolor{keywordflow}{then}
856           \textcolor{comment}{! Estimate the horizontal density gradients along layers.}
857           drdia = drho\_dt\_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
858                   drho\_ds\_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
859           drdib = drho\_dt\_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
860                   drho\_ds\_u(i) * (s(i+1,j,k)-s(i,j,k))
861 
862           \textcolor{comment}{! Estimate the vertical density gradients times the grid spacing.}
863           drdkl = (drho\_dt\_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
864                    drho\_ds\_u(i) * (s(i,j,k)-s(i,j,k-1)))
865           drdkr = (drho\_dt\_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
866                    drho\_ds\_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
867           drdkde\_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
868         \textcolor{keywordflow}{elseif} (find\_work) \textcolor{keywordflow}{then} \textcolor{comment}{! This is used in pure stacked SW mode}
869           drdkde\_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
870 \textcolor{keywordflow}{        endif}
871         \textcolor{keywordflow}{if} (use\_stanley) \textcolor{keywordflow}{then}
872           \textcolor{comment}{! Correction to the horizontal density gradient due to nonlinearity in}
873           \textcolor{comment}{! the EOS rectifying SGS temperature anomalies}
874           drdia = drdia + drho\_dt\_dt\_u(i) * 0.5 * ( tsgs2(i+1,j,k-1)-tsgs2(i,j,k-1) )
875           drdib = drdib + drho\_dt\_dt\_u(i) * 0.5 * ( tsgs2(i+1,j,k)-tsgs2(i,j,k) )
876 \textcolor{keywordflow}{        endif}
877         \textcolor{keywordflow}{if} (find\_work) drdi\_u(i,k) = drdib
878 
879         \textcolor{keywordflow}{if} (k > nk\_linear) \textcolor{keywordflow}{then}
880           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
881             \textcolor{keywordflow}{if} (cs%use\_FGNV\_streamfn .or. find\_work .or. .not.present\_slope\_x) \textcolor{keywordflow}{then}
882               hg2l = h(i,j,k-1)*h(i,j,k) + h\_neglect2
883               hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h\_neglect2
884               hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h\_neglect
885               har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h\_neglect
886               \textcolor{keywordflow}{if} (gv%Boussinesq) \textcolor{keywordflow}{then}
887                 dzal = hal * gv%H\_to\_Z ; dzar = har * gv%H\_to\_Z
888               \textcolor{keywordflow}{else}
889                 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz\_neglect
890                 dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz\_neglect
891 \textcolor{keywordflow}{              endif}
892               \textcolor{comment}{! Use the harmonic mean thicknesses to weight the horizontal gradients.}
893               \textcolor{comment}{! These unnormalized weights have been rearranged to minimize divisions.}
894               wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
895 
896               drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
897               \textcolor{comment}{! The expression for drdz above is mathematically equivalent to:}
898               \textcolor{comment}{!   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &}
899               \textcolor{comment}{!          ((hg2L/haL) + (hg2R/haR))}
900               hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h\_neglect2
901               hg2b = h(i,j,k)*h(i+1,j,k) + h\_neglect2
902               haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h\_neglect
903               hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h\_neglect
904 
905               \textcolor{comment}{! hN2\_u is used with the FGNV streamfunction formulation}
906               hn2\_u(i,k) = (0.5 * gv%H\_to\_Z * ( hg2a / haa + hg2b / hab )) * &
907                            max(drdz*g\_rho0, n2\_floor)
908 \textcolor{keywordflow}{            endif}
909             \textcolor{keywordflow}{if} (present\_slope\_x) \textcolor{keywordflow}{then}
910               slope = slope\_x(i,j,k)
911               slope2\_ratio\_u(i,k) = slope**2 * i\_slope\_max2
912             \textcolor{keywordflow}{else}
913               \textcolor{comment}{! Use the harmonic mean thicknesses to weight the horizontal gradients.}
914               \textcolor{comment}{! These unnormalized weights have been rearranged to minimize divisions.}
915               wta = hg2a*hab ; wtb = hg2b*haa
916               \textcolor{comment}{! This is the gradient of density along geopotentials.}
917               drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
918                       drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
919 
920               \textcolor{comment}{! This estimate of slope is accurate for small slopes, but bounded}
921               \textcolor{comment}{! to be between -1 and 1.}
922               mag\_grad2 = drdx**2 + (us%L\_to\_Z*drdz)**2
923               \textcolor{keywordflow}{if} (mag\_grad2 > 0.0) \textcolor{keywordflow}{then}
924                 slope = drdx / sqrt(mag\_grad2)
925                 slope2\_ratio\_u(i,k) = slope**2 * i\_slope\_max2
926               \textcolor{keywordflow}{else} \textcolor{comment}{! Just in case mag\_grad2 = 0 ever.}
927                 slope = 0.0
928                 slope2\_ratio\_u(i,k) = 1.0e20  \textcolor{comment}{! Force the use of the safe streamfunction.}
929 \textcolor{keywordflow}{              endif}
930 \textcolor{keywordflow}{            endif}
931 
932             \textcolor{comment}{! Adjust real slope by weights that bias towards slope of interfaces}
933             \textcolor{comment}{! that ignore density gradients along layers.}
934             \textcolor{keywordflow}{if} (present\_int\_slope\_u) \textcolor{keywordflow}{then}
935               slope = (1.0 - int\_slope\_u(i,j,k)) * slope + &
936                       int\_slope\_u(i,j,k) * us%Z\_to\_L*((e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j))
937               slope2\_ratio\_u(i,k) = (1.0 - int\_slope\_u(i,j,k)) * slope2\_ratio\_u(i,k)
938 \textcolor{keywordflow}{            endif}
939 
940             slope\_x\_pe(i,j,k) = min(slope,cs%slope\_max)
941             hn2\_x\_pe(i,j,k) = hn2\_u(i,k)
942             \textcolor{keywordflow}{if} (cs%id\_slope\_x > 0) cs%diagSlopeX(i,j,k) = slope
943 
944             \textcolor{comment}{! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1].}
945             sfn\_unlim\_u(i,k) = -((kh\_u(i,j,k)*g%dy\_Cu(i,j))*us%L\_to\_Z*slope)
946 
947             \textcolor{comment}{! Avoid moving dense water upslope from below the level of}
948             \textcolor{comment}{! the bottom on the receiving side.}
949             \textcolor{keywordflow}{if} (sfn\_unlim\_u(i,k) > 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! The flow below this interface is positive.}
950               \textcolor{keywordflow}{if} (e(i,j,k) < e(i+1,j,nz+1)) \textcolor{keywordflow}{then}
951                 sfn\_unlim\_u(i,k) = 0.0 \textcolor{comment}{! This is not uhtot, because it may compensate for}
952                                 \textcolor{comment}{! deeper flow in very unusual cases.}
953               \textcolor{keywordflow}{elseif} (e(i+1,j,nz+1) > e(i,j,k+1)) \textcolor{keywordflow}{then}
954                 \textcolor{comment}{! Scale the transport with the fraction of the donor layer above}
955                 \textcolor{comment}{! the bottom on the receiving side.}
956                 sfn\_unlim\_u(i,k) = sfn\_unlim\_u(i,k) * ((e(i,j,k) - e(i+1,j,nz+1)) / &
957                                          ((e(i,j,k) - e(i,j,k+1)) + dz\_neglect))
958 \textcolor{keywordflow}{              endif}
959             \textcolor{keywordflow}{else}
960               \textcolor{keywordflow}{if} (e(i+1,j,k) < e(i,j,nz+1)) \textcolor{keywordflow}{then} ; sfn\_unlim\_u(i,k) = 0.0
961               \textcolor{keywordflow}{elseif} (e(i,j,nz+1) > e(i+1,j,k+1)) \textcolor{keywordflow}{then}
962                 sfn\_unlim\_u(i,k) = sfn\_unlim\_u(i,k) * ((e(i+1,j,k) - e(i,j,nz+1)) / &
963                                        ((e(i+1,j,k) - e(i+1,j,k+1)) + dz\_neglect))
964 \textcolor{keywordflow}{              endif}
965 \textcolor{keywordflow}{            endif}
966 
967           \textcolor{keywordflow}{else} \textcolor{comment}{! .not. use\_EOS}
968             \textcolor{keywordflow}{if} (present\_slope\_x) \textcolor{keywordflow}{then}
969               slope = slope\_x(i,j,k)
970             \textcolor{keywordflow}{else}
971               slope = us%Z\_to\_L*((e(i,j,k)-e(i+1,j,k))*g%IdxCu(i,j)) * g%mask2dCu(i,j)
972 \textcolor{keywordflow}{            endif}
973             \textcolor{keywordflow}{if} (cs%id\_slope\_x > 0) cs%diagSlopeX(i,j,k) = slope
974             sfn\_unlim\_u(i,k) = ((kh\_u(i,j,k)*g%dy\_Cu(i,j))*us%L\_to\_Z*slope)
975             hn2\_u(i,k) = gv%g\_prime(k)
976 \textcolor{keywordflow}{          endif} \textcolor{comment}{! if (use\_EOS)}
977         \textcolor{keywordflow}{else} \textcolor{comment}{! if (k > nk\_linear)}
978           hn2\_u(i,k) = n2\_floor * dz\_neglect
979           sfn\_unlim\_u(i,k) = 0.
980 \textcolor{keywordflow}{        endif} \textcolor{comment}{! if (k > nk\_linear)}
981         \textcolor{keywordflow}{if} (cs%id\_sfn\_unlim\_x>0) diag\_sfn\_unlim\_x(i,j,k) = sfn\_unlim\_u(i,k)
982 \textcolor{keywordflow}{      enddo} \textcolor{comment}{! i-loop}
983 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! k-loop}
984 
985     \textcolor{keywordflow}{if} (cs%use\_FGNV\_streamfn) \textcolor{keywordflow}{then}
986       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is-1,ie ; \textcolor{keywordflow}{if} (g%mask2dCu(i,j)>0.) \textcolor{keywordflow}{then}
987         h\_harm = max( h\_neglect, &
988               2. * h(i,j,k) * h(i+1,j,k) / ( ( h(i,j,k) + h(i+1,j,k) ) + h\_neglect ) )
989         c2\_h\_u(i,k) = cs%FGNV\_scale * &
990             ( 0.5*( cg1(i,j) + cg1(i+1,j) ) )**2 / (gv%H\_to\_Z*h\_harm)
991 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
992 
993       \textcolor{comment}{! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.}
994       \textcolor{keywordflow}{do} i=is-1,ie
995         \textcolor{keywordflow}{if} (g%mask2dCu(i,j)>0.) \textcolor{keywordflow}{then}
996           sfn\_unlim\_u(i,:) = ( 1. + cs%FGNV\_scale ) * sfn\_unlim\_u(i,:)
997           \textcolor{keyword}{call }streamfn\_solver(nz, c2\_h\_u(i,:), hn2\_u(i,:), sfn\_unlim\_u(i,:))
998         \textcolor{keywordflow}{else}
999           sfn\_unlim\_u(i,:) = 0.
1000 \textcolor{keywordflow}{        endif}
1001 \textcolor{keywordflow}{      enddo}
1002 \textcolor{keywordflow}{    endif}
1003 
1004     \textcolor{keywordflow}{do} k=nz,2,-1
1005       \textcolor{keywordflow}{do} i=is-1,ie
1006         \textcolor{keywordflow}{if} (k > nk\_linear) \textcolor{keywordflow}{then}
1007           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1008 
1009             \textcolor{keywordflow}{if} (uhtot(i,j) <= 0.0) \textcolor{keywordflow}{then}
1010               \textcolor{comment}{! The transport that must balance the transport below is positive.}
1011               sfn\_safe = uhtot(i,j) * (1.0 - h\_frac(i,j,k)) * gv%H\_to\_Z
1012             \textcolor{keywordflow}{else} \textcolor{comment}{!  (uhtot(I,j) > 0.0)}
1013               sfn\_safe = uhtot(i,j) * (1.0 - h\_frac(i+1,j,k)) * gv%H\_to\_Z
1014 \textcolor{keywordflow}{            endif}
1015 
1016             \textcolor{comment}{! The actual streamfunction at each interface.}
1017             sfn\_est = (sfn\_unlim\_u(i,k) + slope2\_ratio\_u(i,k)*sfn\_safe) / (1.0 + slope2\_ratio\_u(i,k))
1018           \textcolor{keywordflow}{else}  \textcolor{comment}{! With .not.use\_EOS, the layers are constant density.}
1019             sfn\_est = sfn\_unlim\_u(i,k)
1020 \textcolor{keywordflow}{          endif}
1021 
1022           \textcolor{comment}{! Make sure that there is enough mass above to allow the streamfunction}
1023           \textcolor{comment}{! to satisfy the boundary condition of 0 at the surface.}
1024           sfn\_in\_h = min(max(sfn\_est * gv%Z\_to\_H, -h\_avail\_rsum(i,j,k)), h\_avail\_rsum(i+1,j,k))
1025 
1026           \textcolor{comment}{! The actual transport is limited by the mass available in the two}
1027           \textcolor{comment}{! neighboring grid cells.}
1028           uhd(i,j,k) = max(min((sfn\_in\_h - uhtot(i,j)), h\_avail(i,j,k)), &
1029                            -h\_avail(i+1,j,k))
1030 
1031           \textcolor{keywordflow}{if} (cs%id\_sfn\_x>0) diag\_sfn\_x(i,j,k) = diag\_sfn\_x(i,j,k+1) + uhd(i,j,k)
1032 \textcolor{comment}{!         sfn\_x(I,j,K) = max(min(Sfn\_in\_h, uhtot(I,j)+h\_avail(i,j,k)), &}
1033 \textcolor{comment}{!                            uhtot(I,j)-h\_avail(i+1,j,K))}
1034 \textcolor{comment}{!         sfn\_slope\_x(I,j,K) = max(uhtot(I,j)-h\_avail(i+1,j,k), &}
1035 \textcolor{comment}{!                                  min(uhtot(I,j)+h\_avail(i,j,k), &}
1036 \textcolor{comment}{!               min(h\_avail\_rsum(i+1,j,K), max(-h\_avail\_rsum(i,j,K), &}
1037 \textcolor{comment}{!               (KH\_u(I,j,K)*G%dy\_Cu(I,j)) * ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) )) ))}
1038         \textcolor{keywordflow}{else} \textcolor{comment}{! k <= nk\_linear}
1039           \textcolor{comment}{! Balance the deeper flow with a return flow uniformly distributed}
1040           \textcolor{comment}{! though the remaining near-surface layers.  This is the same as}
1041           \textcolor{comment}{! using Sfn\_safe above.  There is no need to apply the limiters in}
1042           \textcolor{comment}{! this case.}
1043           \textcolor{keywordflow}{if} (uhtot(i,j) <= 0.0) \textcolor{keywordflow}{then}
1044             uhd(i,j,k) = -uhtot(i,j) * h\_frac(i,j,k)
1045           \textcolor{keywordflow}{else} \textcolor{comment}{!  (uhtot(I,j) > 0.0)}
1046             uhd(i,j,k) = -uhtot(i,j) * h\_frac(i+1,j,k)
1047 \textcolor{keywordflow}{          endif}
1048 
1049           \textcolor{keywordflow}{if} (cs%id\_sfn\_x>0) diag\_sfn\_x(i,j,k) = diag\_sfn\_x(i,j,k+1) + uhd(i,j,k)
1050 \textcolor{comment}{!         if (sfn\_slope\_x(I,j,K+1) <= 0.0) then}
1051 \textcolor{comment}{!           sfn\_slope\_x(I,j,K) = sfn\_slope\_x(I,j,K+1) * (1.0 - h\_frac(i,j,k))}
1052 \textcolor{comment}{!         else}
1053 \textcolor{comment}{!           sfn\_slope\_x(I,j,K) = sfn\_slope\_x(I,j,K+1) * (1.0 - h\_frac(i+1,j,k))}
1054 \textcolor{comment}{!         endif}
1055 \textcolor{keywordflow}{        endif}
1056 
1057         uhtot(i,j) = uhtot(i,j) + uhd(i,j,k)
1058 
1059         \textcolor{keywordflow}{if} (find\_work) \textcolor{keywordflow}{then}
1060           \textcolor{comment}{!   This is the energy tendency based on the original profiles, and does}
1061           \textcolor{comment}{! not include any nonlinear terms due to a finite time step (which would}
1062           \textcolor{comment}{! involve interactions between the fluxes through the different faces.}
1063           \textcolor{comment}{!   A second order centered estimate is used for the density transfered}
1064           \textcolor{comment}{! between water columns.}
1065 
1066           work\_u(i,j) = work\_u(i,j) + g\_scale * &
1067             ( uhtot(i,j) * drdkde\_u(i,k) - &
1068               (uhd(i,j,k) * drdi\_u(i,k)) * 0.25 * &
1069               ((e(i,j,k) + e(i,j,k+1)) + (e(i+1,j,k) + e(i+1,j,k+1))) )
1070 \textcolor{keywordflow}{        endif}
1071 
1072 \textcolor{keywordflow}{      enddo}
1073 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! end of k-loop}
1074 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end of j-loop}
1075 
1076     \textcolor{comment}{! Calculate the meridional fluxes and gradients.}
1077     eosdom\_v(:) = eos\_domain(g%HI)
1078 \textcolor{comment}{!$OMP parallel do default(none) shared(nz,is,ie,js,je,find\_work,use\_EOS,G,GV,US,pres,T,S, &}
1079 \textcolor{comment}{!$OMP                                  nk\_linear,IsdB,tv,h,h\_neglect,e,dz\_neglect,  &}
1080 \textcolor{comment}{!$OMP                                  I\_slope\_max2,h\_neglect2,present\_int\_slope\_v, &}
1081 \textcolor{comment}{!$OMP                                  int\_slope\_v,KH\_v,vhtot,h\_frac,h\_avail\_rsum,  &}
1082 \textcolor{comment}{!$OMP                                  vhD,h\_avail,G\_scale,Work\_v,CS,slope\_y,cg1,   &}
1083 \textcolor{comment}{!$OMP                                  diag\_sfn\_y,diag\_sfn\_unlim\_y,N2\_floor,EOSdom\_v,&}
1084 \textcolor{comment}{!$OMP                                  use\_stanley, Tsgs2,                          &}
1085 \textcolor{comment}{!$OMP                                  present\_slope\_y,G\_rho0,Slope\_y\_PE,hN2\_y\_PE)  &}
1086 \textcolor{comment}{!$OMP                          private(drdjA,drdjB,drdkL,drdkR,pres\_v,T\_v,S\_v,      &}
1087 \textcolor{comment}{!$OMP                                  drho\_dT\_v,drho\_dS\_v,hg2A,hg2B,hg2L,hg2R,haA, &}
1088 \textcolor{comment}{!$OMP                                  drho\_dT\_dT\_v,scrap,                          &}
1089 \textcolor{comment}{!$OMP                                  haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,  &}
1090 \textcolor{comment}{!$OMP                                  drdy,mag\_grad2,Slope,slope2\_Ratio\_v,hN2\_v,   &}
1091 \textcolor{comment}{!$OMP                                  Sfn\_unlim\_v,drdj\_v,drdkDe\_v,h\_harm,c2\_h\_v,   &}
1092 \textcolor{comment}{!$OMP                                  Sfn\_safe,Sfn\_est,Sfn\_in\_h,calc\_derivatives)}
1093   \textcolor{keywordflow}{do} j=js-1,je
1094     \textcolor{keywordflow}{do} k=nz,2,-1
1095       \textcolor{keywordflow}{if} (find\_work .and. .not.(use\_eos)) \textcolor{keywordflow}{then}
1096         drdja = 0.0 ; drdjb = 0.0
1097         drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
1098 \textcolor{keywordflow}{      endif}
1099 
1100       calc\_derivatives = use\_eos .and. (k >= nk\_linear) .and. &
1101          (find\_work .or. .not. present\_slope\_y .or. cs%use\_FGNV\_streamfn .or. use\_stanley)
1102 
1103       \textcolor{keywordflow}{if} (calc\_derivatives) \textcolor{keywordflow}{then}
1104         \textcolor{keywordflow}{do} i=is,ie
1105           pres\_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
1106           t\_v(i) = 0.25*((t(i,j,k) + t(i,j+1,k)) + (t(i,j,k-1) + t(i,j+1,k-1)))
1107           s\_v(i) = 0.25*((s(i,j,k) + s(i,j+1,k)) + (s(i,j,k-1) + s(i,j+1,k-1)))
1108 \textcolor{keywordflow}{        enddo}
1109         \textcolor{keyword}{call }calculate\_density\_derivs(t\_v, s\_v, pres\_v, drho\_dt\_v, drho\_ds\_v, &
1110                                       tv%eqn\_of\_state, eosdom\_v)
1111 \textcolor{keywordflow}{      endif}
1112       \textcolor{keywordflow}{if} (use\_stanley) \textcolor{keywordflow}{then}
1113         \textcolor{comment}{! The second line below would correspond to arguments}
1114         \textcolor{comment}{!            drho\_dS\_dS, drho\_dS\_dT, drho\_dT\_dT, drho\_dS\_dP, drho\_dT\_dP, &}
1115         \textcolor{keyword}{call }calculate\_density\_second\_derivs(t\_v, s\_v, pres\_v, &
1116                      scrap, scrap, drho\_dt\_dt\_v, scrap, scrap, &
1117                      is, ie-is+1, tv%eqn\_of\_state)
1118 \textcolor{keywordflow}{      endif}
1119       \textcolor{keywordflow}{do} i=is,ie
1120         \textcolor{keywordflow}{if} (calc\_derivatives) \textcolor{keywordflow}{then}
1121           \textcolor{comment}{! Estimate the horizontal density gradients along layers.}
1122           drdja = drho\_dt\_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
1123                   drho\_ds\_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
1124           drdjb = drho\_dt\_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
1125                   drho\_ds\_v(i) * (s(i,j+1,k)-s(i,j,k))
1126 
1127           \textcolor{comment}{! Estimate the vertical density gradients times the grid spacing.}
1128           drdkl = (drho\_dt\_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
1129                    drho\_ds\_v(i) * (s(i,j,k)-s(i,j,k-1)))
1130           drdkr = (drho\_dt\_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
1131                    drho\_ds\_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
1132           drdkde\_v(i,k) =  drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1133         \textcolor{keywordflow}{elseif} (find\_work) \textcolor{keywordflow}{then} \textcolor{comment}{! This is used in pure stacked SW mode}
1134           drdkde\_v(i,k) =  drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1135 \textcolor{keywordflow}{        endif}
1136         \textcolor{keywordflow}{if} (use\_stanley) \textcolor{keywordflow}{then}
1137           \textcolor{comment}{! Correction to the horizontal density gradient due to nonlinearity in}
1138           \textcolor{comment}{! the EOS rectifying SGS temperature anomalies}
1139           drdja = drdja + drho\_dt\_dt\_v(i) * 0.5 * ( tsgs2(i,j+1,k-1)-tsgs2(i,j,k-1) )
1140           drdjb = drdjb + drho\_dt\_dt\_v(i) * 0.5 * ( tsgs2(i,j+1,k)-tsgs2(i,j,k) )
1141 \textcolor{keywordflow}{        endif}
1142 
1143         \textcolor{keywordflow}{if} (find\_work) drdj\_v(i,k) = drdjb
1144 
1145         \textcolor{keywordflow}{if} (k > nk\_linear) \textcolor{keywordflow}{then}
1146           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1147             \textcolor{keywordflow}{if} (cs%use\_FGNV\_streamfn .or. find\_work .or. .not. present\_slope\_y) \textcolor{keywordflow}{then}
1148               hg2l = h(i,j,k-1)*h(i,j,k) + h\_neglect2
1149               hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h\_neglect2
1150               hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h\_neglect
1151               har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h\_neglect
1152               \textcolor{keywordflow}{if} (gv%Boussinesq) \textcolor{keywordflow}{then}
1153                 dzal = hal * gv%H\_to\_Z ; dzar = har * gv%H\_to\_Z
1154               \textcolor{keywordflow}{else}
1155                 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz\_neglect
1156                 dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz\_neglect
1157 \textcolor{keywordflow}{              endif}
1158               \textcolor{comment}{! Use the harmonic mean thicknesses to weight the horizontal gradients.}
1159               \textcolor{comment}{! These unnormalized weights have been rearranged to minimize divisions.}
1160               wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
1161 
1162               drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
1163               \textcolor{comment}{! The expression for drdz above is mathematically equivalent to:}
1164               \textcolor{comment}{!   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &}
1165               \textcolor{comment}{!          ((hg2L/haL) + (hg2R/haR))}
1166               hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h\_neglect2
1167               hg2b = h(i,j,k)*h(i,j+1,k) + h\_neglect2
1168               haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h\_neglect
1169               hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h\_neglect
1170 
1171               \textcolor{comment}{! hN2\_v is used with the FGNV streamfunction formulation}
1172               hn2\_v(i,k) = (0.5 * gv%H\_to\_Z * ( hg2a / haa + hg2b / hab )) * &
1173                            max(drdz*g\_rho0, n2\_floor)
1174 \textcolor{keywordflow}{            endif}
1175             \textcolor{keywordflow}{if} (present\_slope\_y) \textcolor{keywordflow}{then}
1176               slope = slope\_y(i,j,k)
1177               slope2\_ratio\_v(i,k) = slope**2 * i\_slope\_max2
1178             \textcolor{keywordflow}{else}
1179               \textcolor{comment}{! Use the harmonic mean thicknesses to weight the horizontal gradients.}
1180               \textcolor{comment}{! These unnormalized weights have been rearranged to minimize divisions.}
1181               wta = hg2a*hab ; wtb = hg2b*haa
1182               \textcolor{comment}{! This is the gradient of density along geopotentials.}
1183               drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
1184                       drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
1185 
1186               \textcolor{comment}{! This estimate of slope is accurate for small slopes, but bounded}
1187               \textcolor{comment}{! to be between -1 and 1.}
1188               mag\_grad2 = drdy**2 + (us%L\_to\_Z*drdz)**2
1189               \textcolor{keywordflow}{if} (mag\_grad2 > 0.0) \textcolor{keywordflow}{then}
1190                 slope = drdy / sqrt(mag\_grad2)
1191                 slope2\_ratio\_v(i,k) = slope**2 * i\_slope\_max2
1192               \textcolor{keywordflow}{else} \textcolor{comment}{! Just in case mag\_grad2 = 0 ever.}
1193                 slope = 0.0
1194                 slope2\_ratio\_v(i,k) = 1.0e20  \textcolor{comment}{! Force the use of the safe streamfunction.}
1195 \textcolor{keywordflow}{              endif}
1196 \textcolor{keywordflow}{            endif}
1197 
1198             \textcolor{comment}{! Adjust real slope by weights that bias towards slope of interfaces}
1199             \textcolor{comment}{! that ignore density gradients along layers.}
1200             \textcolor{keywordflow}{if} (present\_int\_slope\_v) \textcolor{keywordflow}{then}
1201               slope = (1.0 - int\_slope\_v(i,j,k)) * slope + &
1202                       int\_slope\_v(i,j,k) * us%Z\_to\_L*((e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j))
1203               slope2\_ratio\_v(i,k) = (1.0 - int\_slope\_v(i,j,k)) * slope2\_ratio\_v(i,k)
1204 \textcolor{keywordflow}{            endif}
1205 
1206             slope\_y\_pe(i,j,k) = min(slope,cs%slope\_max)
1207             hn2\_y\_pe(i,j,k) = hn2\_v(i,k)
1208             \textcolor{keywordflow}{if} (cs%id\_slope\_y > 0) cs%diagSlopeY(i,j,k) = slope
1209 
1210             \textcolor{comment}{! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1].}
1211             sfn\_unlim\_v(i,k) = -((kh\_v(i,j,k)*g%dx\_Cv(i,j))*us%L\_to\_Z*slope)
1212 
1213             \textcolor{comment}{! Avoid moving dense water upslope from below the level of}
1214             \textcolor{comment}{! the bottom on the receiving side.}
1215             \textcolor{keywordflow}{if} (sfn\_unlim\_v(i,k) > 0.0) \textcolor{keywordflow}{then} \textcolor{comment}{! The flow below this interface is positive.}
1216               \textcolor{keywordflow}{if} (e(i,j,k) < e(i,j+1,nz+1)) \textcolor{keywordflow}{then}
1217                 sfn\_unlim\_v(i,k) = 0.0 \textcolor{comment}{! This is not vhtot, because it may compensate for}
1218                                 \textcolor{comment}{! deeper flow in very unusual cases.}
1219               \textcolor{keywordflow}{elseif} (e(i,j+1,nz+1) > e(i,j,k+1)) \textcolor{keywordflow}{then}
1220                 \textcolor{comment}{! Scale the transport with the fraction of the donor layer above}
1221                 \textcolor{comment}{! the bottom on the receiving side.}
1222                 sfn\_unlim\_v(i,k) = sfn\_unlim\_v(i,k) * ((e(i,j,k) - e(i,j+1,nz+1)) / &
1223                                          ((e(i,j,k) - e(i,j,k+1)) + dz\_neglect))
1224 \textcolor{keywordflow}{              endif}
1225             \textcolor{keywordflow}{else}
1226               \textcolor{keywordflow}{if} (e(i,j+1,k) < e(i,j,nz+1)) \textcolor{keywordflow}{then} ; sfn\_unlim\_v(i,k) = 0.0
1227               \textcolor{keywordflow}{elseif} (e(i,j,nz+1) > e(i,j+1,k+1)) \textcolor{keywordflow}{then}
1228                 sfn\_unlim\_v(i,k) = sfn\_unlim\_v(i,k) * ((e(i,j+1,k) - e(i,j,nz+1)) / &
1229                                        ((e(i,j+1,k) - e(i,j+1,k+1)) + dz\_neglect))
1230 \textcolor{keywordflow}{              endif}
1231 \textcolor{keywordflow}{            endif}
1232 
1233           \textcolor{keywordflow}{else} \textcolor{comment}{! .not. use\_EOS}
1234             \textcolor{keywordflow}{if} (present\_slope\_y) \textcolor{keywordflow}{then}
1235               slope = slope\_y(i,j,k)
1236             \textcolor{keywordflow}{else}
1237               slope = us%Z\_to\_L*((e(i,j,k)-e(i,j+1,k))*g%IdyCv(i,j)) * g%mask2dCv(i,j)
1238 \textcolor{keywordflow}{            endif}
1239             \textcolor{keywordflow}{if} (cs%id\_slope\_y > 0) cs%diagSlopeY(i,j,k) = slope
1240             sfn\_unlim\_v(i,k) = ((kh\_v(i,j,k)*g%dx\_Cv(i,j))*us%L\_to\_Z*slope)
1241             hn2\_v(i,k) = gv%g\_prime(k)
1242 \textcolor{keywordflow}{          endif} \textcolor{comment}{! if (use\_EOS)}
1243         \textcolor{keywordflow}{else} \textcolor{comment}{! if (k > nk\_linear)}
1244           hn2\_v(i,k) = n2\_floor * dz\_neglect
1245           sfn\_unlim\_v(i,k) = 0.
1246 \textcolor{keywordflow}{        endif} \textcolor{comment}{! if (k > nk\_linear)}
1247         \textcolor{keywordflow}{if} (cs%id\_sfn\_unlim\_y>0) diag\_sfn\_unlim\_y(i,j,k) = sfn\_unlim\_v(i,k)
1248 \textcolor{keywordflow}{      enddo} \textcolor{comment}{! i-loop}
1249 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! k-loop}
1250 
1251     \textcolor{keywordflow}{if} (cs%use\_FGNV\_streamfn) \textcolor{keywordflow}{then}
1252       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (g%mask2dCv(i,j)>0.) \textcolor{keywordflow}{then}
1253         h\_harm = max( h\_neglect, &
1254               2. * h(i,j,k) * h(i,j+1,k) / ( ( h(i,j,k) + h(i,j+1,k) ) + h\_neglect ) )
1255         c2\_h\_v(i,k) = cs%FGNV\_scale * &
1256             ( 0.5*( cg1(i,j) + cg1(i,j+1) ) )**2 / (gv%H\_to\_Z*h\_harm)
1257 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1258 
1259       \textcolor{comment}{! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.}
1260       \textcolor{keywordflow}{do} i=is,ie
1261         \textcolor{keywordflow}{if} (g%mask2dCv(i,j)>0.) \textcolor{keywordflow}{then}
1262           sfn\_unlim\_v(i,:) = ( 1. + cs%FGNV\_scale ) * sfn\_unlim\_v(i,:)
1263           \textcolor{keyword}{call }streamfn\_solver(nz, c2\_h\_v(i,:), hn2\_v(i,:), sfn\_unlim\_v(i,:))
1264         \textcolor{keywordflow}{else}
1265           sfn\_unlim\_v(i,:) = 0.
1266 \textcolor{keywordflow}{        endif}
1267 \textcolor{keywordflow}{      enddo}
1268 \textcolor{keywordflow}{    endif}
1269 
1270     \textcolor{keywordflow}{do} k=nz,2,-1
1271       \textcolor{keywordflow}{do} i=is,ie
1272         \textcolor{keywordflow}{if} (k > nk\_linear) \textcolor{keywordflow}{then}
1273           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1274 
1275             \textcolor{keywordflow}{if} (vhtot(i,j) <= 0.0) \textcolor{keywordflow}{then}
1276               \textcolor{comment}{! The transport that must balance the transport below is positive.}
1277               sfn\_safe = vhtot(i,j) * (1.0 - h\_frac(i,j,k)) * gv%H\_to\_Z
1278             \textcolor{keywordflow}{else} \textcolor{comment}{!  (vhtot(I,j) > 0.0)}
1279               sfn\_safe = vhtot(i,j) * (1.0 - h\_frac(i,j+1,k)) * gv%H\_to\_Z
1280 \textcolor{keywordflow}{            endif}
1281 
1282             \textcolor{comment}{! The actual streamfunction at each interface.}
1283             sfn\_est = (sfn\_unlim\_v(i,k) + slope2\_ratio\_v(i,k)*sfn\_safe) / (1.0 + slope2\_ratio\_v(i,k))
1284           \textcolor{keywordflow}{else}      \textcolor{comment}{! With .not.use\_EOS, the layers are constant density.}
1285             sfn\_est = sfn\_unlim\_v(i,k)
1286 \textcolor{keywordflow}{          endif}
1287 
1288           \textcolor{comment}{! Make sure that there is enough mass above to allow the streamfunction}
1289           \textcolor{comment}{! to satisfy the boundary condition of 0 at the surface.}
1290           sfn\_in\_h = min(max(sfn\_est * gv%Z\_to\_H, -h\_avail\_rsum(i,j,k)), h\_avail\_rsum(i,j+1,k))
1291 
1292           \textcolor{comment}{! The actual transport is limited by the mass available in the two}
1293           \textcolor{comment}{! neighboring grid cells.}
1294           vhd(i,j,k) = max(min((sfn\_in\_h - vhtot(i,j)), h\_avail(i,j,k)), -h\_avail(i,j+1,k))
1295 
1296           \textcolor{keywordflow}{if} (cs%id\_sfn\_y>0) diag\_sfn\_y(i,j,k) = diag\_sfn\_y(i,j,k+1) + vhd(i,j,k)
1297 \textcolor{comment}{!         sfn\_y(i,J,K) = max(min(Sfn\_in\_h, vhtot(i,J)+h\_avail(i,j,k)), &}
1298 \textcolor{comment}{!                            vhtot(i,J)-h\_avail(i,j+1,k))}
1299 \textcolor{comment}{!         sfn\_slope\_y(i,J,K) = max(vhtot(i,J)-h\_avail(i,j+1,k), &}
1300 \textcolor{comment}{!                                  min(vhtot(i,J)+h\_avail(i,j,k), &}
1301 \textcolor{comment}{!               min(h\_avail\_rsum(i,j+1,K), max(-h\_avail\_rsum(i,j,K), &}
1302 \textcolor{comment}{!               (KH\_v(i,J,K)*G%dx\_Cv(i,J)) * ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) )) ))}
1303         \textcolor{keywordflow}{else} \textcolor{comment}{! k <= nk\_linear}
1304           \textcolor{comment}{! Balance the deeper flow with a return flow uniformly distributed}
1305           \textcolor{comment}{! though the remaining near-surface layers.  This is the same as}
1306           \textcolor{comment}{! using Sfn\_safe above.  There is no need to apply the limiters in}
1307           \textcolor{comment}{! this case.}
1308           \textcolor{keywordflow}{if} (vhtot(i,j) <= 0.0) \textcolor{keywordflow}{then}
1309             vhd(i,j,k) = -vhtot(i,j) * h\_frac(i,j,k)
1310           \textcolor{keywordflow}{else} \textcolor{comment}{!  (vhtot(i,J) > 0.0)}
1311             vhd(i,j,k) = -vhtot(i,j) * h\_frac(i,j+1,k)
1312 \textcolor{keywordflow}{          endif}
1313 
1314           \textcolor{keywordflow}{if} (cs%id\_sfn\_y>0) diag\_sfn\_y(i,j,k) = diag\_sfn\_y(i,j,k+1) + vhd(i,j,k)
1315 \textcolor{comment}{!         if (sfn\_slope\_y(i,J,K+1) <= 0.0) then}
1316 \textcolor{comment}{!           sfn\_slope\_y(i,J,K) = sfn\_slope\_y(i,J,K+1) * (1.0 - h\_frac(i,j,k))}
1317 \textcolor{comment}{!         else}
1318 \textcolor{comment}{!           sfn\_slope\_y(i,J,K) = sfn\_slope\_y(i,J,K+1) * (1.0 - h\_frac(i,j+1,k))}
1319 \textcolor{comment}{!         endif}
1320 \textcolor{keywordflow}{        endif}
1321 
1322         vhtot(i,j) = vhtot(i,j)  + vhd(i,j,k)
1323 
1324         \textcolor{keywordflow}{if} (find\_work) \textcolor{keywordflow}{then}
1325           \textcolor{comment}{!   This is the energy tendency based on the original profiles, and does}
1326           \textcolor{comment}{! not include any nonlinear terms due to a finite time step (which would}
1327           \textcolor{comment}{! involve interactions between the fluxes through the different faces.}
1328           \textcolor{comment}{!   A second order centered estimate is used for the density transfered}
1329           \textcolor{comment}{! between water columns.}
1330 
1331           work\_v(i,j) = work\_v(i,j) + g\_scale * &
1332             ( vhtot(i,j) * drdkde\_v(i,k) - &
1333              (vhd(i,j,k) * drdj\_v(i,k)) * 0.25 * &
1334              ((e(i,j,k) + e(i,j,k+1)) + (e(i,j+1,k) + e(i,j+1,k+1))) )
1335 \textcolor{keywordflow}{        endif}
1336 
1337 \textcolor{keywordflow}{      enddo}
1338 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! end of k-loop}
1339 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end of j-loop}
1340 
1341   \textcolor{comment}{! In layer 1, enforce the boundary conditions that Sfn(z=0) = 0.0}
1342   \textcolor{keywordflow}{if} (.not.find\_work .or. .not.(use\_eos)) \textcolor{keywordflow}{then}
1343     \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is-1,ie ; uhd(i,j,1) = -uhtot(i,j) ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1344     \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is,ie ; vhd(i,j,1) = -vhtot(i,j) ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1345   \textcolor{keywordflow}{else}
1346     eosdom\_u(1) = (is-1) - (g%IsdB-1) ; eosdom\_u(2) = ie - (g%IsdB-1)
1347     \textcolor{comment}{!$OMP parallel do default(shared) private(pres\_u,T\_u,S\_u,drho\_dT\_u,drho\_dS\_u,drdiB)}
1348     \textcolor{keywordflow}{do} j=js,je
1349       \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1350         \textcolor{keywordflow}{do} i=is-1,ie
1351           pres\_u(i) = 0.5*(pres(i,j,1) + pres(i+1,j,1))
1352           t\_u(i) = 0.5*(t(i,j,1) + t(i+1,j,1))
1353           s\_u(i) = 0.5*(s(i,j,1) + s(i+1,j,1))
1354 \textcolor{keywordflow}{        enddo}
1355         \textcolor{keyword}{call }calculate\_density\_derivs(t\_u, s\_u, pres\_u, drho\_dt\_u, drho\_ds\_u, &
1356                                       tv%eqn\_of\_state, eosdom\_u )
1357 \textcolor{keywordflow}{      endif}
1358       \textcolor{keywordflow}{do} i=is-1,ie
1359         uhd(i,j,1) = -uhtot(i,j)
1360 
1361         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1362           drdib = drho\_dt\_u(i) * (t(i+1,j,1)-t(i,j,1)) + &
1363                   drho\_ds\_u(i) * (s(i+1,j,1)-s(i,j,1))
1364 \textcolor{keywordflow}{        endif}
1365         \textcolor{keywordflow}{if} (cs%use\_GM\_work\_bug) \textcolor{keywordflow}{then}
1366           work\_u(i,j) = work\_u(i,j) + g\_scale * &
1367               ( (uhd(i,j,1) * drdib) * 0.25 * &
1368                 ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1369         \textcolor{keywordflow}{else}
1370           work\_u(i,j) = work\_u(i,j) - g\_scale * &
1371               ( (uhd(i,j,1) * drdib) * 0.25 * &
1372                 ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1373 \textcolor{keywordflow}{        endif}
1374 \textcolor{keywordflow}{      enddo}
1375 \textcolor{keywordflow}{    enddo}
1376 
1377     eosdom\_v(:) = eos\_domain(g%HI)
1378     \textcolor{comment}{!$OMP parallel do default(shared) private(pres\_v,T\_v,S\_v,drho\_dT\_v,drho\_dS\_v,drdjB)}
1379     \textcolor{keywordflow}{do} j=js-1,je
1380       \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1381         \textcolor{keywordflow}{do} i=is,ie
1382           pres\_v(i) = 0.5*(pres(i,j,1) + pres(i,j+1,1))
1383           t\_v(i) = 0.5*(t(i,j,1) + t(i,j+1,1))
1384           s\_v(i) = 0.5*(s(i,j,1) + s(i,j+1,1))
1385 \textcolor{keywordflow}{        enddo}
1386         \textcolor{keyword}{call }calculate\_density\_derivs(t\_v, s\_v, pres\_v, drho\_dt\_v, drho\_ds\_v, &
1387                                       tv%eqn\_of\_state, eosdom\_v)
1388 \textcolor{keywordflow}{      endif}
1389       \textcolor{keywordflow}{do} i=is,ie
1390         vhd(i,j,1) = -vhtot(i,j)
1391 
1392         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1393           drdjb = drho\_dt\_v(i) * (t(i,j+1,1)-t(i,j,1)) + &
1394                   drho\_ds\_v(i) * (s(i,j+1,1)-s(i,j,1))
1395 \textcolor{keywordflow}{        endif}
1396         work\_v(i,j) = work\_v(i,j) - g\_scale * &
1397             ( (vhd(i,j,1) * drdjb) * 0.25 * &
1398               ((e(i,j,1) + e(i,j,2)) + (e(i,j+1,1) + e(i,j+1,2))) )
1399 \textcolor{keywordflow}{      enddo}
1400 \textcolor{keywordflow}{    enddo}
1401 \textcolor{keywordflow}{  endif}
1402 
1403   \textcolor{keywordflow}{if} (find\_work) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
1404     \textcolor{comment}{! Note that the units of Work\_v and Work\_u are W, while Work\_h is W m-2.}
1405     work\_h = 0.5 * g%IareaT(i,j) * &
1406       ((work\_u(i-1,j) + work\_u(i,j)) + (work\_v(i,j-1) + work\_v(i,j)))
1407     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs%GMwork)) cs%GMwork(i,j) = work\_h
1408     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke) .and. .not.cs%GM\_src\_alt) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke%GM\_src)) \textcolor{keywordflow}{then}
1409       meke%GM\_src(i,j) = meke%GM\_src(i,j) + work\_h
1410 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
1411 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
1412 
1413   \textcolor{keywordflow}{if} (find\_work .and. cs%GM\_src\_alt .and. \textcolor{keyword}{associated}(meke)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(meke%GM\_src)) \textcolor{keywordflow}{then}
1414     \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{do} k=nz,1,-1
1415       pe\_release\_h = -0.25*(kh\_u(i,j,k)*(slope\_x\_pe(i,j,k)**2) * hn2\_x\_pe(i,j,k) + &
1416                             kh\_u(i-1,j,k)*(slope\_x\_pe(i-1,j,k)**2) * hn2\_x\_pe(i-1,j,k) + &
1417                             kh\_v(i,j,k)*(slope\_y\_pe(i,j,k)**2) * hn2\_y\_pe(i,j,k) + &
1418                             kh\_v(i,j-1,k)*(slope\_y\_pe(i,j-1,k)**2) * hn2\_y\_pe(i,j-1,k))
1419       meke%GM\_src(i,j) = meke%GM\_src(i,j) + us%L\_to\_Z**2 * gv%Rho0 * pe\_release\_h
1420 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1421 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
1422 
1423   \textcolor{keywordflow}{if} (cs%id\_slope\_x > 0) \textcolor{keyword}{call }post\_data(cs%id\_slope\_x, cs%diagSlopeX, cs%diag)
1424   \textcolor{keywordflow}{if} (cs%id\_slope\_y > 0) \textcolor{keyword}{call }post\_data(cs%id\_slope\_y, cs%diagSlopeY, cs%diag)
1425   \textcolor{keywordflow}{if} (cs%id\_sfn\_x > 0) \textcolor{keyword}{call }post\_data(cs%id\_sfn\_x, diag\_sfn\_x, cs%diag)
1426   \textcolor{keywordflow}{if} (cs%id\_sfn\_y > 0) \textcolor{keyword}{call }post\_data(cs%id\_sfn\_y, diag\_sfn\_y, cs%diag)
1427   \textcolor{keywordflow}{if} (cs%id\_sfn\_unlim\_x > 0) \textcolor{keyword}{call }post\_data(cs%id\_sfn\_unlim\_x, diag\_sfn\_unlim\_x, cs%diag)
1428   \textcolor{keywordflow}{if} (cs%id\_sfn\_unlim\_y > 0) \textcolor{keyword}{call }post\_data(cs%id\_sfn\_unlim\_y, diag\_sfn\_unlim\_y, cs%diag)
1429 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__thickness__diffuse_a96e66e8491141614027c1583aeaecd7a}\label{namespacemom__thickness__diffuse_a96e66e8491141614027c1583aeaecd7a}} 
\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}!thickness\+\_\+diffuse\+\_\+get\+\_\+kh@{thickness\+\_\+diffuse\+\_\+get\+\_\+kh}}
\index{thickness\+\_\+diffuse\+\_\+get\+\_\+kh@{thickness\+\_\+diffuse\+\_\+get\+\_\+kh}!mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}
\subsubsection{\texorpdfstring{thickness\+\_\+diffuse\+\_\+get\+\_\+kh()}{thickness\_diffuse\_get\_kh()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+thickness\+\_\+diffuse\+::thickness\+\_\+diffuse\+\_\+get\+\_\+kh (\begin{DoxyParamCaption}\item[{type(\hyperlink{structmom__thickness__diffuse_1_1thickness__diffuse__cs}{thickness\+\_\+diffuse\+\_\+cs}), pointer}]{CS,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)+1), intent(inout)}]{K\+H\+\_\+u\+\_\+\+G\+ME,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)+1), intent(inout)}]{K\+H\+\_\+v\+\_\+\+G\+ME,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G }\end{DoxyParamCaption})}



Copies ubtav and vbtav from private type into arrays. 


\begin{DoxyParams}[1]{Parameters}
 & {\em cs} & Control structure for this module\\
\hline
\mbox{\tt in}  & {\em g} & Grid structure\\
\hline
\mbox{\tt in,out}  & {\em kh\+\_\+u\+\_\+gme} & interface height diffusivities at u-\/faces \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em kh\+\_\+v\+\_\+gme} & interface height diffusivities at v-\/faces \mbox{[}L2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 2091 of file M\+O\+M\+\_\+thickness\+\_\+diffuse.\+F90.


\begin{DoxyCode}
2091   \textcolor{keywordtype}{type}(thickness\_diffuse\_cs),          \textcolor{keywordtype}{pointer}     :: cs\textcolor{comment}{   !< Control structure for}
2092 \textcolor{comment}{                                                   !! this module}
2093   \textcolor{keywordtype}{type}(ocean\_grid\_type),               \textcolor{keywordtype}{intent(in)}  :: g\textcolor{comment}{    !< Grid structure}
2094   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(inout)} :: kh\_u\_gme\textcolor{comment}{ !< interface height}
2095 \textcolor{comment}{                                                   !! diffusivities at u-faces [L2 T-1 ~> m2 s-1]}
2096   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G)+1)}, \textcolor{keywordtype}{intent(inout)} :: kh\_v\_gme\textcolor{comment}{ !< interface height}
2097 \textcolor{comment}{                                                   !! diffusivities at v-faces [L2 T-1 ~> m2 s-1]}
2098   \textcolor{comment}{! Local variables}
2099   \textcolor{keywordtype}{integer} :: i,j,k
2100 
2101   \textcolor{keywordflow}{do} k=1,g%ke+1 ; \textcolor{keywordflow}{do} j = g%jsc, g%jec ; \textcolor{keywordflow}{do} i = g%isc-1, g%iec
2102     kh\_u\_gme(i,j,k) = cs%KH\_u\_GME(i,j,k)
2103 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
2104 
2105   \textcolor{keywordflow}{do} k=1,g%ke+1 ; \textcolor{keywordflow}{do} j = g%jsc-1, g%jec ; \textcolor{keywordflow}{do} i = g%isc, g%iec
2106     kh\_v\_gme(i,j,k) = cs%KH\_v\_GME(i,j,k)
2107 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
2108 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__thickness__diffuse_affd209ab72d54799680c78d8ba1ad779}\label{namespacemom__thickness__diffuse_affd209ab72d54799680c78d8ba1ad779}} 
\index{mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}!thickness\+\_\+diffuse\+\_\+init@{thickness\+\_\+diffuse\+\_\+init}}
\index{thickness\+\_\+diffuse\+\_\+init@{thickness\+\_\+diffuse\+\_\+init}!mom\+\_\+thickness\+\_\+diffuse@{mom\+\_\+thickness\+\_\+diffuse}}
\subsubsection{\texorpdfstring{thickness\+\_\+diffuse\+\_\+init()}{thickness\_diffuse\_init()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+thickness\+\_\+diffuse\+::thickness\+\_\+diffuse\+\_\+init (\begin{DoxyParamCaption}\item[{type(time\+\_\+type), intent(in)}]{Time,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(param\+\_\+file\+\_\+type), intent(in)}]{param\+\_\+file,  }\item[{type(diag\+\_\+ctrl), intent(inout), target}]{diag,  }\item[{type(cont\+\_\+diag\+\_\+ptrs), intent(inout)}]{C\+Dp,  }\item[{type(\hyperlink{structmom__thickness__diffuse_1_1thickness__diffuse__cs}{thickness\+\_\+diffuse\+\_\+cs}), pointer}]{CS }\end{DoxyParamCaption})}



Initialize the thickness diffusion module/structure. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em time} & Current model time\\
\hline
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em param\+\_\+file} & Parameter file handles\\
\hline
\mbox{\tt in,out}  & {\em diag} & Diagnostics control structure\\
\hline
\mbox{\tt in,out}  & {\em cdp} & Continuity equation diagnostics\\
\hline
 & {\em cs} & Control structure for thickness diffusion \\
\hline
\end{DoxyParams}


Definition at line 1882 of file M\+O\+M\+\_\+thickness\+\_\+diffuse.\+F90.


\begin{DoxyCode}
1882   \textcolor{keywordtype}{type}(time\_type),         \textcolor{keywordtype}{intent(in)} :: time\textcolor{comment}{    !< Current model time}
1883   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{       !< Ocean grid structure}
1884   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)} :: gv\textcolor{comment}{      !< Vertical grid structure}
1885   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)} :: us\textcolor{comment}{      !< A dimensional unit scaling type}
1886   \textcolor{keywordtype}{type}(param\_file\_type),   \textcolor{keywordtype}{intent(in)} :: param\_file\textcolor{comment}{ !< Parameter file handles}
1887   \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(inout)} :: diag\textcolor{comment}{ !< Diagnostics control structure}
1888   \textcolor{keywordtype}{type}(cont\_diag\_ptrs),    \textcolor{keywordtype}{intent(inout)} :: cdp\textcolor{comment}{  !< Continuity equation diagnostics}
1889   \textcolor{keywordtype}{type}(thickness\_diffuse\_cs), \textcolor{keywordtype}{pointer}    :: cs\textcolor{comment}{   !< Control structure for thickness diffusion}
1890 
1891 \textcolor{comment}{! This include declares and sets the variable "version".}
1892 \textcolor{preprocessor}{#include "version\_variable.h"}
1893 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_thickness\_diffuse"} \textcolor{comment}{! This module's name.}
1894   \textcolor{keywordtype}{real} :: omega        \textcolor{comment}{! The Earth's rotation rate [T-1 ~> s-1]}
1895   \textcolor{keywordtype}{real} :: strat\_floor  \textcolor{comment}{! A floor for Brunt-Vasaila frequency in the Ferrari et al. 2010,}
1896                        \textcolor{comment}{! streamfunction formulation, expressed as a fraction of planetary}
1897                        \textcolor{comment}{! rotation [nondim].}
1898   \textcolor{keywordtype}{logical} :: default\_2018\_answers \textcolor{comment}{! The default setting for the various 2018\_ANSWERS flags.}
1899 
1900   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}
1901     \textcolor{keyword}{call }mom\_error(warning, &
1902       \textcolor{stringliteral}{"Thickness\_diffuse\_init called with an associated control structure."})
1903     \textcolor{keywordflow}{return}
1904   \textcolor{keywordflow}{else} ; \textcolor{keyword}{allocate}(cs) ;\textcolor{keywordflow}{ endif}
1905 
1906   cs%diag => diag
1907 
1908   \textcolor{comment}{! Read all relevant parameters and write them to the model log.}
1909   \textcolor{keyword}{call }log\_version(param\_file, mdl, version, \textcolor{stringliteral}{""})
1910   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"THICKNESSDIFFUSE"}, cs%thickness\_diffuse, &
1911                  \textcolor{stringliteral}{"If true, interface heights are diffused with a "}//&
1912                  \textcolor{stringliteral}{"coefficient of KHTH."}, default=.false.)
1913   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KHTH"}, cs%Khth, &
1914                  \textcolor{stringliteral}{"The background horizontal thickness diffusivity."}, &
1915                  default=0.0, units=\textcolor{stringliteral}{"m2 s-1"}, scale=us%m\_to\_L**2*us%T\_to\_s)
1916   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KHTH\_SLOPE\_CFF"}, cs%KHTH\_Slope\_Cff, &
1917                  \textcolor{stringliteral}{"The nondimensional coefficient in the Visbeck formula "}//&
1918                  \textcolor{stringliteral}{"for the interface depth diffusivity"}, units=\textcolor{stringliteral}{"nondim"}, &
1919                  default=0.0)
1920   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KHTH\_MIN"}, cs%KHTH\_Min, &
1921                  \textcolor{stringliteral}{"The minimum horizontal thickness diffusivity."}, &
1922                  default=0.0, units=\textcolor{stringliteral}{"m2 s-1"}, scale=us%m\_to\_L**2*us%T\_to\_s)
1923   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KHTH\_MAX"}, cs%KHTH\_Max, &
1924                  \textcolor{stringliteral}{"The maximum horizontal thickness diffusivity."}, &
1925                  default=0.0, units=\textcolor{stringliteral}{"m2 s-1"}, scale=us%m\_to\_L**2*us%T\_to\_s)
1926   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KHTH\_MAX\_CFL"}, cs%max\_Khth\_CFL, &
1927                  \textcolor{stringliteral}{"The maximum value of the local diffusive CFL ratio that "}//&
1928                  \textcolor{stringliteral}{"is permitted for the thickness diffusivity. 1.0 is the "}//&
1929                  \textcolor{stringliteral}{"marginally unstable value in a pure layered model, but "}//&
1930                  \textcolor{stringliteral}{"much smaller numbers (e.g. 0.1) seem to work better for "}//&
1931                  \textcolor{stringliteral}{"ALE-based models."}, units = \textcolor{stringliteral}{"nondimensional"}, default=0.8)
1932 
1933 \textcolor{comment}{!  call get\_param(param\_file, mdl, "USE\_QG\_LEITH\_GM", CS%QG\_Leith\_GM, &}
1934 \textcolor{comment}{!               "If true, use the QG Leith viscosity as the GM coefficient.", &}
1935 \textcolor{comment}{!               default=.false.)}
1936 
1937   \textcolor{keywordflow}{if} (cs%max\_Khth\_CFL < 0.0) cs%max\_Khth\_CFL = 0.0
1938   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DETANGLE\_INTERFACES"}, cs%detangle\_interfaces, &
1939                  \textcolor{stringliteral}{"If defined add 3-d structured enhanced interface height "}//&
1940                  \textcolor{stringliteral}{"diffusivities to horizontally smooth jagged layers."}, &
1941                  default=.false.)
1942   cs%detangle\_time = 0.0
1943   \textcolor{keywordflow}{if} (cs%detangle\_interfaces) &
1944     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DETANGLE\_TIMESCALE"}, cs%detangle\_time, &
1945                  \textcolor{stringliteral}{"A timescale over which maximally jagged grid-scale "}//&
1946                  \textcolor{stringliteral}{"thickness variations are suppressed.  This must be "}//&
1947                  \textcolor{stringliteral}{"longer than DT, or 0 to use DT."}, units=\textcolor{stringliteral}{"s"}, default=0.0, scale=us%s\_to\_T)
1948   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KHTH\_SLOPE\_MAX"}, cs%slope\_max, &
1949                  \textcolor{stringliteral}{"A slope beyond which the calculated isopycnal slope is "}//&
1950                  \textcolor{stringliteral}{"not reliable and is scaled away."}, units=\textcolor{stringliteral}{"nondim"}, default=0.01)
1951   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KD\_SMOOTH"}, cs%kappa\_smooth, &
1952                  \textcolor{stringliteral}{"A diapycnal diffusivity that is used to interpolate "}//&
1953                  \textcolor{stringliteral}{"more sensible values of T & S into thin layers."}, &
1954                  units=\textcolor{stringliteral}{"m2 s-1"}, default=1.0e-6, scale=us%m\_to\_Z**2*us%T\_to\_s)
1955   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KHTH\_USE\_FGNV\_STREAMFUNCTION"}, cs%use\_FGNV\_streamfn, &
1956                  \textcolor{stringliteral}{"If true, use the streamfunction formulation of "}//&
1957                  \textcolor{stringliteral}{"Ferrari et al., 2010, which effectively emphasizes "}//&
1958                  \textcolor{stringliteral}{"graver vertical modes by smoothing in the vertical."},  &
1959                  default=.false.)
1960   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"FGNV\_FILTER\_SCALE"}, cs%FGNV\_scale, &
1961                  \textcolor{stringliteral}{"A coefficient scaling the vertical smoothing term in the "}//&
1962                  \textcolor{stringliteral}{"Ferrari et al., 2010, streamfunction formulation."}, &
1963                  units=\textcolor{stringliteral}{"nondim"}, default=1., do\_not\_log=.not.cs%use\_FGNV\_streamfn)
1964   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"FGNV\_C\_MIN"}, cs%FGNV\_c\_min, &
1965                  \textcolor{stringliteral}{"A minium wave speed used in the Ferrari et al., 2010, "}//&
1966                  \textcolor{stringliteral}{"streamfunction formulation."}, &
1967                  default=0., units=\textcolor{stringliteral}{"m s-1"}, scale=us%m\_s\_to\_L\_T, do\_not\_log=.not.cs%use\_FGNV\_streamfn)
1968   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"FGNV\_STRAT\_FLOOR"}, strat\_floor, &
1969                  \textcolor{stringliteral}{"A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010, "}//&
1970                  \textcolor{stringliteral}{"streamfunction formulation, expressed as a fraction of planetary "}//&
1971                  \textcolor{stringliteral}{"rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy."}, &
1972                  default=1.e-15, units=\textcolor{stringliteral}{"nondim"}, do\_not\_log=.not.cs%use\_FGNV\_streamfn)
1973   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"STANLEY\_PRM\_DET\_COEFF"}, cs%Stanley\_det\_coeff, &
1974                  \textcolor{stringliteral}{"The coefficient correlating SGS temperature variance with the mean "}//&
1975                  \textcolor{stringliteral}{"temperature gradient in the deterministic part of the Stanley parameterization. "}//&
1976                  \textcolor{stringliteral}{"Negative values disable the scheme."}, units=\textcolor{stringliteral}{"nondim"}, default=-1.0)
1977   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"OMEGA"}, omega, &
1978                  \textcolor{stringliteral}{"The rotation rate of the earth."}, &
1979                  default=7.2921e-5, units=\textcolor{stringliteral}{"s-1"}, scale=us%T\_to\_s, do\_not\_log=.not.cs%use\_FGNV\_streamfn)
1980   \textcolor{keywordflow}{if} (cs%use\_FGNV\_streamfn) cs%N2\_floor = (strat\_floor*omega)**2
1981   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DEBUG"}, cs%debug, &
1982                  \textcolor{stringliteral}{"If true, write out verbose debugging data."}, &
1983                  default=.false., debuggingparam=.true.)
1984 
1985   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MEKE\_GM\_SRC\_ALT"}, cs%GM\_src\_alt, &
1986                  \textcolor{stringliteral}{"If true, use the GM energy conversion form S^2*N^2*kappa rather "}//&
1987                  \textcolor{stringliteral}{"than the streamfunction for the GM source term."}, default=.false.)
1988   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MEKE\_GEOMETRIC"}, cs%MEKE\_GEOMETRIC, &
1989                  \textcolor{stringliteral}{"If true, uses the GM coefficient formulation from the GEOMETRIC "}//&
1990                  \textcolor{stringliteral}{"framework (Marshall et al., 2012)."}, default=.false.)
1991   \textcolor{keywordflow}{if} (cs%MEKE\_GEOMETRIC) \textcolor{keywordflow}{then}
1992     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MEKE\_GEOMETRIC\_EPSILON"}, cs%MEKE\_GEOMETRIC\_epsilon, &
1993                  \textcolor{stringliteral}{"Minimum Eady growth rate used in the calculation of GEOMETRIC "}//&
1994                  \textcolor{stringliteral}{"thickness diffusivity."}, units=\textcolor{stringliteral}{"s-1"}, default=1.0e-7, scale=us%T\_to\_s)
1995     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MEKE\_GEOMETRIC\_ALPHA"}, cs%MEKE\_GEOMETRIC\_alpha, &
1996                  \textcolor{stringliteral}{"The nondimensional coefficient governing the efficiency of the GEOMETRIC "}//&
1997                  \textcolor{stringliteral}{"thickness diffusion."}, units=\textcolor{stringliteral}{"nondim"}, default=0.05)
1998 
1999     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DEFAULT\_2018\_ANSWERS"}, default\_2018\_answers, &
2000                  \textcolor{stringliteral}{"This sets the default value for the various \_2018\_ANSWERS parameters."}, &
2001                  default=.false.)
2002     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MEKE\_GEOMETRIC\_2018\_ANSWERS"}, cs%MEKE\_GEOM\_answers\_2018, &
2003                  \textcolor{stringliteral}{"If true, use expressions in the MEKE\_GEOMETRIC calculation that recover the "}//&
2004                  \textcolor{stringliteral}{"answers from the original implementation.  Otherwise, use expressions that "}//&
2005                  \textcolor{stringliteral}{"satisfy rotational symmetry."}, default=default\_2018\_answers)
2006 \textcolor{keywordflow}{  endif}
2007 
2008   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_KH\_IN\_MEKE"}, cs%Use\_KH\_in\_MEKE, &
2009                  \textcolor{stringliteral}{"If true, uses the thickness diffusivity calculated here to diffuse MEKE."}, &
2010                  default=.false.)
2011 
2012   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_GME"}, cs%use\_GME\_thickness\_diffuse, &
2013                  \textcolor{stringliteral}{"If true, use the GM+E backscatter scheme in association "}//&
2014                  \textcolor{stringliteral}{"with the Gent and McWilliams parameterization."}, default=.false.)
2015 
2016   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_GM\_WORK\_BUG"}, cs%use\_GM\_work\_bug, &
2017                  \textcolor{stringliteral}{"If true, compute the top-layer work tendency on the u-grid "}//&
2018                  \textcolor{stringliteral}{"with the incorrect sign, for legacy reproducibility."}, &
2019                  default=.false.)
2020 
2021   \textcolor{keywordflow}{if} (cs%use\_GME\_thickness\_diffuse) \textcolor{keywordflow}{then}
2022     \textcolor{keyword}{call }safe\_alloc\_ptr(cs%KH\_u\_GME,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
2023     \textcolor{keyword}{call }safe\_alloc\_ptr(cs%KH\_v\_GME,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
2024 \textcolor{keywordflow}{  endif}
2025 
2026   cs%id\_uhGM = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'uhGM'}, diag%axesCuL, time, &
2027            \textcolor{stringliteral}{'Time Mean Diffusive Zonal Thickness Flux'}, &
2028            \textcolor{stringliteral}{'kg s-1'}, conversion=gv%H\_to\_kg\_m2*us%L\_to\_m**2*us%s\_to\_T, &
2029            y\_cell\_method=\textcolor{stringliteral}{'sum'}, v\_extensive=.true.)
2030   \textcolor{keywordflow}{if} (cs%id\_uhGM > 0) \textcolor{keyword}{call }safe\_alloc\_ptr(cdp%uhGM,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
2031   cs%id\_vhGM = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'vhGM'}, diag%axesCvL, time, &
2032            \textcolor{stringliteral}{'Time Mean Diffusive Meridional Thickness Flux'}, &
2033            \textcolor{stringliteral}{'kg s-1'}, conversion=gv%H\_to\_kg\_m2*us%L\_to\_m**2*us%s\_to\_T, &
2034            x\_cell\_method=\textcolor{stringliteral}{'sum'}, v\_extensive=.true.)
2035   \textcolor{keywordflow}{if} (cs%id\_vhGM > 0) \textcolor{keyword}{call }safe\_alloc\_ptr(cdp%vhGM,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
2036 
2037   cs%id\_GMwork = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'GMwork'}, diag%axesT1, time, &
2038           \textcolor{stringliteral}{'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection'}, &
2039           \textcolor{stringliteral}{'W m-2'}, conversion=us%RZ3\_T3\_to\_W\_m2*us%L\_to\_Z**2, cmor\_field\_name=\textcolor{stringliteral}{'tnkebto'}, &
2040           cmor\_long\_name=\textcolor{stringliteral}{'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection'}
      , &
2041           cmor\_standard\_name=\textcolor{stringliteral}{
      'tendency\_of\_ocean\_eddy\_kinetic\_energy\_content\_due\_to\_parameterized\_eddy\_advection'})
2042   \textcolor{keywordflow}{if} (cs%id\_GMwork > 0) \textcolor{keyword}{call }safe\_alloc\_ptr(cs%GMwork,g%isd,g%ied,g%jsd,g%jed)
2043 
2044   cs%id\_KH\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'KHTH\_u'}, diag%axesCui, time, &
2045            \textcolor{stringliteral}{'Parameterized mesoscale eddy advection diffusivity at U-point'}, &
2046            \textcolor{stringliteral}{'m2 s-1'}, conversion=us%L\_to\_m**2*us%s\_to\_T)
2047   cs%id\_KH\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'KHTH\_v'}, diag%axesCvi, time, &
2048            \textcolor{stringliteral}{'Parameterized mesoscale eddy advection diffusivity at V-point'}, &
2049            \textcolor{stringliteral}{'m2 s-1'}, conversion=us%L\_to\_m**2*us%s\_to\_T)
2050   cs%id\_KH\_t = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'KHTH\_t'}, diag%axesTL, time, &
2051           \textcolor{stringliteral}{'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection'}, &
2052           \textcolor{stringliteral}{'m2 s-1'}, conversion=us%L\_to\_m**2*us%s\_to\_T, &
2053           cmor\_field\_name=\textcolor{stringliteral}{'diftrblo'}, &
2054           cmor\_long\_name=\textcolor{stringliteral}{'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection'}, &
2055           cmor\_units=\textcolor{stringliteral}{'m2 s-1'}, &
2056           cmor\_standard\_name=\textcolor{stringliteral}{'ocean\_tracer\_diffusivity\_due\_to\_parameterized\_mesoscale\_advection'})
2057 
2058   cs%id\_KH\_u1 = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'KHTH\_u1'}, diag%axesCu1, time,         &
2059            \textcolor{stringliteral}{'Parameterized mesoscale eddy advection diffusivity at U-points (2-D)'}, &
2060            \textcolor{stringliteral}{'m2 s-1'}, conversion=us%L\_to\_m**2*us%s\_to\_T)
2061   cs%id\_KH\_v1 = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'KHTH\_v1'}, diag%axesCv1, time,         &
2062            \textcolor{stringliteral}{'Parameterized mesoscale eddy advection diffusivity at V-points (2-D)'}, &
2063            \textcolor{stringliteral}{'m2 s-1'}, conversion=us%L\_to\_m**2*us%s\_to\_T)
2064   cs%id\_KH\_t1 = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'KHTH\_t1'}, diag%axesT1, time, &
2065            \textcolor{stringliteral}{'Parameterized mesoscale eddy advection diffusivity at T-points (2-D)'}, &
2066            \textcolor{stringliteral}{'m2 s-1'}, conversion=us%L\_to\_m**2*us%s\_to\_T)
2067 
2068   cs%id\_slope\_x =  register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'neutral\_slope\_x'}, diag%axesCui, time, &
2069            \textcolor{stringliteral}{'Zonal slope of neutral surface'}, \textcolor{stringliteral}{'nondim'})
2070   \textcolor{keywordflow}{if} (cs%id\_slope\_x > 0) \textcolor{keyword}{call }safe\_alloc\_ptr(cs%diagSlopeX,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
2071   cs%id\_slope\_y =  register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'neutral\_slope\_y'}, diag%axesCvi, time, &
2072            \textcolor{stringliteral}{'Meridional slope of neutral surface'}, \textcolor{stringliteral}{'nondim'})
2073   \textcolor{keywordflow}{if} (cs%id\_slope\_y > 0) \textcolor{keyword}{call }safe\_alloc\_ptr(cs%diagSlopeY,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
2074   cs%id\_sfn\_x =  register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'GM\_sfn\_x'}, diag%axesCui, time, &
2075            \textcolor{stringliteral}{'Parameterized Zonal Overturning Streamfunction'}, &
2076            \textcolor{stringliteral}{'m3 s-1'}, conversion=gv%H\_to\_m*us%L\_to\_m**2*us%s\_to\_T)
2077   cs%id\_sfn\_y =  register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'GM\_sfn\_y'}, diag%axesCvi, time, &
2078            \textcolor{stringliteral}{'Parameterized Meridional Overturning Streamfunction'}, &
2079            \textcolor{stringliteral}{'m3 s-1'}, conversion=gv%H\_to\_m*us%L\_to\_m**2*us%s\_to\_T)
2080   cs%id\_sfn\_unlim\_x =  register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'GM\_sfn\_unlim\_x'}, diag%axesCui, time, &
2081            \textcolor{stringliteral}{'Parameterized Zonal Overturning Streamfunction before limiting/smoothing'}, &
2082            \textcolor{stringliteral}{'m3 s-1'}, conversion=us%Z\_to\_m*us%L\_to\_m**2*us%s\_to\_T)
2083   cs%id\_sfn\_unlim\_y =  register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'GM\_sfn\_unlim\_y'}, diag%axesCvi, time, &
2084            \textcolor{stringliteral}{'Parameterized Meridional Overturning Streamfunction before limiting/smoothing'}, &
2085            \textcolor{stringliteral}{'m3 s-1'}, conversion=us%Z\_to\_m*us%L\_to\_m**2*us%s\_to\_T)
2086 
\end{DoxyCode}