\hypertarget{namespacemom__set__visc}{}\section{mom\+\_\+set\+\_\+visc Module Reference}
\label{namespacemom__set__visc}\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}


\subsection{Detailed Description}
Calculates various values related to the bottom boundary layer, such as the viscosity and thickness of the B\+BL (set\+\_\+viscous\+\_\+\+B\+BL). 

This would also be the module in which other viscous quantities that are flow-\/independent might be set. This information is transmitted to other modules via a vertvisc type structure.

The same code is used for the two velocity components, by indirectly referencing the velocities and defining a handful of direction-\/specific defined variables. \subsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \hyperlink{structmom__set__visc_1_1set__visc__cs}{set\+\_\+visc\+\_\+cs}
\begin{DoxyCompactList}\small\item\em Control structure for M\+O\+M\+\_\+set\+\_\+visc. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \hyperlink{namespacemom__set__visc_a9865113fe07928e7a240c2868ed45e5f}{set\+\_\+viscous\+\_\+bbl} (u, v, h, tv, visc, G, GV, US, CS, symmetrize)
\begin{DoxyCompactList}\small\item\em Calculates the thickness of the bottom boundary layer and the viscosity within that layer. \end{DoxyCompactList}\item 
real function \hyperlink{namespacemom__set__visc_a0356a4e81cca9f7f31bbf87c717a6600}{set\+\_\+v\+\_\+at\+\_\+u} (v, h, G, i, j, k, mask2d\+Cv, O\+BC)
\begin{DoxyCompactList}\small\item\em This subroutine finds a thickness-\/weighted value of v at the u-\/points. \end{DoxyCompactList}\item 
real function \hyperlink{namespacemom__set__visc_a46583b82467e74d8654c3c0a037a25cd}{set\+\_\+u\+\_\+at\+\_\+v} (u, h, G, i, j, k, mask2d\+Cu, O\+BC)
\begin{DoxyCompactList}\small\item\em This subroutine finds a thickness-\/weighted value of u at the v-\/points. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__set__visc_aba41cd4f8baa1cda9036d97087ce8a22}{set\+\_\+viscous\+\_\+ml} (u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetrize)
\begin{DoxyCompactList}\small\item\em Calculates the thickness of the surface boundary layer for applying an elevated viscosity. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__set__visc_ae2d9d9f74c1e9aec56257cfad372b0fd}{set\+\_\+visc\+\_\+register\+\_\+restarts} (HI, GV, param\+\_\+file, visc, restart\+\_\+\+CS)
\begin{DoxyCompactList}\small\item\em Register any fields associated with the vertvisc\+\_\+type. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__set__visc_a2ec9be1a61c7d4a062aab4ff03e3ca29}{set\+\_\+visc\+\_\+init} (Time, G, GV, US, param\+\_\+file, diag, visc, CS, restart\+\_\+\+CS, O\+BC)
\begin{DoxyCompactList}\small\item\em Initializes the M\+O\+M\+\_\+set\+\_\+visc control structure. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespacemom__set__visc_abbac7a7410557759c74444787de415a9}{set\+\_\+visc\+\_\+end} (visc, CS)
\begin{DoxyCompactList}\small\item\em This subroutine dellocates any memory in the set\+\_\+visc control structure. \end{DoxyCompactList}\end{DoxyCompactItemize}


\subsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespacemom__set__visc_a46583b82467e74d8654c3c0a037a25cd}\label{namespacemom__set__visc_a46583b82467e74d8654c3c0a037a25cd}} 
\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}!set\+\_\+u\+\_\+at\+\_\+v@{set\+\_\+u\+\_\+at\+\_\+v}}
\index{set\+\_\+u\+\_\+at\+\_\+v@{set\+\_\+u\+\_\+at\+\_\+v}!mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}
\subsubsection{\texorpdfstring{set\+\_\+u\+\_\+at\+\_\+v()}{set\_u\_at\_v()}}
{\footnotesize\ttfamily real function mom\+\_\+set\+\_\+visc\+::set\+\_\+u\+\_\+at\+\_\+v (\begin{DoxyParamCaption}\item[{real, dimension( g \%isdb\+: g \%iedb, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{u,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{integer, intent(in)}]{i,  }\item[{integer, intent(in)}]{j,  }\item[{integer, intent(in)}]{k,  }\item[{real, dimension( g \%isdb\+: g \%iedb, g \%jsd\+: g \%jed), intent(in)}]{mask2d\+Cu,  }\item[{type(ocean\+\_\+obc\+\_\+type), pointer}]{O\+BC }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



This subroutine finds a thickness-\/weighted value of u at the v-\/points. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em u} & The zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]} or other units.\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em i} & The i-\/index of the u-\/location to work on.\\
\hline
\mbox{\tt in}  & {\em j} & The j-\/index of the u-\/location to work on.\\
\hline
\mbox{\tt in}  & {\em k} & The k-\/index of the u-\/location to work on.\\
\hline
\mbox{\tt in}  & {\em mask2dcu} & A multiplicative mask of the u-\/points\\
\hline
 & {\em obc} & A pointer to an open boundary condition structure\\
\hline
\end{DoxyParams}
\begin{DoxyReturn}{Returns}
The return value of u at v points in the same units as u, i.\+e. \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]} or other units. 
\end{DoxyReturn}


Definition at line 1160 of file M\+O\+M\+\_\+set\+\_\+viscosity.\+F90.


\begin{DoxyCode}
1160   \textcolor{keywordtype}{type}(ocean\_grid\_type),                     \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{    !< The ocean's grid structure}
1161   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))}, &
1162                          \textcolor{keywordtype}{intent(in)} :: u\textcolor{comment}{    !< The zonal velocity [L T-1 ~> m s-1] or other units.}
1163   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
1164                          \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2]}
1165   \textcolor{keywordtype}{integer},               \textcolor{keywordtype}{intent(in)} :: i\textcolor{comment}{    !< The i-index of the u-location to work on.}
1166   \textcolor{keywordtype}{integer},               \textcolor{keywordtype}{intent(in)} :: j\textcolor{comment}{    !< The j-index of the u-location to work on.}
1167   \textcolor{keywordtype}{integer},               \textcolor{keywordtype}{intent(in)} :: k\textcolor{comment}{    !< The k-index of the u-location to work on.}
1168   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G))}, &
1169                          \textcolor{keywordtype}{intent(in)} :: mask2dcu\textcolor{comment}{ !< A multiplicative mask of the u-points}
1170   \textcolor{keywordtype}{type}(ocean\_obc\_type),  \textcolor{keywordtype}{pointer}    :: obc\textcolor{comment}{  !< A pointer to an open boundary condition structure}
1171   \textcolor{keywordtype}{real}                              :: set\_u\_at\_v\textcolor{comment}{ !< The return value of u at v points in the}
1172 \textcolor{comment}{                                            !! same units as u, i.e. [L T-1 ~> m s-1] or other units.}
1173 
1174   \textcolor{comment}{! This subroutine finds a thickness-weighted value of u at the v-points.}
1175   \textcolor{keywordtype}{real} :: hwt(-1:0,0:1)    \textcolor{comment}{! Masked weights used to average u onto v [H ~> m or kg m-2].}
1176   \textcolor{keywordtype}{real} :: hwt\_tot          \textcolor{comment}{! The sum of the masked thicknesses [H ~> m or kg m-2].}
1177   \textcolor{keywordtype}{integer} :: i0, j0, i1, j1
1178 
1179   \textcolor{keywordflow}{do} j0 = 0,1 ; \textcolor{keywordflow}{do} i0 = -1,0 ; i1 = i+i0 ; j1 = j+j0
1180     hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dcu(i1,j1)
1181 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
1182 
1183   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (obc%number\_of\_segments > 0) \textcolor{keywordflow}{then}
1184     \textcolor{keywordflow}{do} j0 = 0,1 ; \textcolor{keywordflow}{do} i0 = -1,0 ; \textcolor{keywordflow}{if} ((obc%segnum\_u(i+i0,j+j0) /= obc\_none)) \textcolor{keywordflow}{then}
1185       i1 = i+i0 ; j1 = j+j0
1186       \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_u(i1,j1))%direction == obc\_direction\_e) \textcolor{keywordflow}{then}
1187         hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcu(i1,j1)
1188       \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_u(i1,j1))%direction == obc\_direction\_w) \textcolor{keywordflow}{then}
1189         hwt(i0,j0) = 2.0 * h(i1+1,j1,k) * mask2dcu(i1,j1)
1190 \textcolor{keywordflow}{      endif}
1191 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1192 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
1193 
1194   hwt\_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
1195   set\_u\_at\_v = 0.0
1196   \textcolor{keywordflow}{if} (hwt\_tot > 0.0) set\_u\_at\_v = &
1197           ((hwt(0,0) * u(i,j,k) + hwt(-1,1) * u(i-1,j+1,k)) + &
1198            (hwt(-1,0) * u(i-1,j,k) + hwt(0,1) * u(i,j+1,k))) / hwt\_tot
1199 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__set__visc_a0356a4e81cca9f7f31bbf87c717a6600}\label{namespacemom__set__visc_a0356a4e81cca9f7f31bbf87c717a6600}} 
\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}!set\+\_\+v\+\_\+at\+\_\+u@{set\+\_\+v\+\_\+at\+\_\+u}}
\index{set\+\_\+v\+\_\+at\+\_\+u@{set\+\_\+v\+\_\+at\+\_\+u}!mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}
\subsubsection{\texorpdfstring{set\+\_\+v\+\_\+at\+\_\+u()}{set\_v\_at\_u()}}
{\footnotesize\ttfamily real function mom\+\_\+set\+\_\+visc\+::set\+\_\+v\+\_\+at\+\_\+u (\begin{DoxyParamCaption}\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)), intent(in)}]{v,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{integer, intent(in)}]{i,  }\item[{integer, intent(in)}]{j,  }\item[{integer, intent(in)}]{k,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g)), intent(in)}]{mask2d\+Cv,  }\item[{type(ocean\+\_\+obc\+\_\+type), pointer}]{O\+BC }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



This subroutine finds a thickness-\/weighted value of v at the u-\/points. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & The ocean\textquotesingle{}s grid structure\\
\hline
\mbox{\tt in}  & {\em v} & The meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em i} & The i-\/index of the u-\/location to work on.\\
\hline
\mbox{\tt in}  & {\em j} & The j-\/index of the u-\/location to work on.\\
\hline
\mbox{\tt in}  & {\em k} & The k-\/index of the u-\/location to work on.\\
\hline
\mbox{\tt in}  & {\em mask2dcv} & A multiplicative mask of the v-\/points\\
\hline
 & {\em obc} & A pointer to an open boundary condition structure\\
\hline
\end{DoxyParams}
\begin{DoxyReturn}{Returns}
The return value of v at u points points in the same units as u, i.\+e. \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]} or other units. 
\end{DoxyReturn}


Definition at line 1116 of file M\+O\+M\+\_\+set\+\_\+viscosity.\+F90.


\begin{DoxyCode}
1116   \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)} :: g\textcolor{comment}{    !< The ocean's grid structure}
1117   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))}, &
1118                          \textcolor{keywordtype}{intent(in)} :: v\textcolor{comment}{    !< The meridional velocity [L T-1 ~> m s-1]}
1119   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))}, &
1120                          \textcolor{keywordtype}{intent(in)} :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2]}
1121   \textcolor{keywordtype}{integer},               \textcolor{keywordtype}{intent(in)} :: i\textcolor{comment}{    !< The i-index of the u-location to work on.}
1122   \textcolor{keywordtype}{integer},               \textcolor{keywordtype}{intent(in)} :: j\textcolor{comment}{    !< The j-index of the u-location to work on.}
1123   \textcolor{keywordtype}{integer},               \textcolor{keywordtype}{intent(in)} :: k\textcolor{comment}{    !< The k-index of the u-location to work on.}
1124   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G))},&
1125                          \textcolor{keywordtype}{intent(in)} :: mask2dcv\textcolor{comment}{ !< A multiplicative mask of the v-points}
1126   \textcolor{keywordtype}{type}(ocean\_obc\_type),  \textcolor{keywordtype}{pointer}    :: obc\textcolor{comment}{  !< A pointer to an open boundary condition structure}
1127   \textcolor{keywordtype}{real}                              :: set\_v\_at\_u\textcolor{comment}{ !< The return value of v at u points points in the}
1128 \textcolor{comment}{                                            !! same units as u, i.e. [L T-1 ~> m s-1] or other units.}
1129 
1130   \textcolor{comment}{! This subroutine finds a thickness-weighted value of v at the u-points.}
1131   \textcolor{keywordtype}{real} :: hwt(0:1,-1:0)    \textcolor{comment}{! Masked weights used to average u onto v [H ~> m or kg m-2].}
1132   \textcolor{keywordtype}{real} :: hwt\_tot          \textcolor{comment}{! The sum of the masked thicknesses [H ~> m or kg m-2].}
1133   \textcolor{keywordtype}{integer} :: i0, j0, i1, j1
1134 
1135   \textcolor{keywordflow}{do} j0 = -1,0 ; \textcolor{keywordflow}{do} i0 = 0,1 ; i1 = i+i0 ; j1 = j+j0
1136     hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dcv(i1,j1)
1137 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
1138 
1139   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (obc%number\_of\_segments > 0) \textcolor{keywordflow}{then}
1140     \textcolor{keywordflow}{do} j0 = -1,0 ; \textcolor{keywordflow}{do} i0 = 0,1 ; \textcolor{keywordflow}{if} ((obc%segnum\_v(i+i0,j+j0) /= obc\_none)) \textcolor{keywordflow}{then}
1141       i1 = i+i0 ; j1 = j+j0
1142       \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_v(i1,j1))%direction == obc\_direction\_n) \textcolor{keywordflow}{then}
1143         hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcv(i1,j1)
1144       \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_v(i1,j1))%direction == obc\_direction\_s) \textcolor{keywordflow}{then}
1145         hwt(i0,j0) = 2.0 * h(i1,j1+1,k) * mask2dcv(i1,j1)
1146 \textcolor{keywordflow}{      endif}
1147 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1148 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
1149 
1150   hwt\_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
1151   set\_v\_at\_u = 0.0
1152   \textcolor{keywordflow}{if} (hwt\_tot > 0.0) set\_v\_at\_u = &
1153           ((hwt(0,0) * v(i,j,k) + hwt(1,-1) * v(i+1,j-1,k)) + &
1154            (hwt(1,0) * v(i+1,j,k) + hwt(0,-1) * v(i,j-1,k))) / hwt\_tot
1155 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__set__visc_abbac7a7410557759c74444787de415a9}\label{namespacemom__set__visc_abbac7a7410557759c74444787de415a9}} 
\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}!set\+\_\+visc\+\_\+end@{set\+\_\+visc\+\_\+end}}
\index{set\+\_\+visc\+\_\+end@{set\+\_\+visc\+\_\+end}!mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}
\subsubsection{\texorpdfstring{set\+\_\+visc\+\_\+end()}{set\_visc\_end()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+set\+\_\+visc\+::set\+\_\+visc\+\_\+end (\begin{DoxyParamCaption}\item[{type(vertvisc\+\_\+type), intent(inout)}]{visc,  }\item[{type(\hyperlink{structmom__set__visc_1_1set__visc__cs}{set\+\_\+visc\+\_\+cs}), pointer}]{CS }\end{DoxyParamCaption})}



This subroutine dellocates any memory in the set\+\_\+visc control structure. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in,out}  & {\em visc} & A structure containing vertical viscosities and related fields. Elements are deallocated here.\\
\hline
 & {\em cs} & The control structure returned by a previous call to set\+\_\+visc\+\_\+init. \\
\hline
\end{DoxyParams}


Definition at line 2315 of file M\+O\+M\+\_\+set\+\_\+viscosity.\+F90.


\begin{DoxyCode}
2315   \textcolor{keywordtype}{type}(vertvisc\_type), \textcolor{keywordtype}{intent(inout)} :: visc\textcolor{comment}{ !< A structure containing vertical viscosities and}
2316 \textcolor{comment}{                                             !! related fields.  Elements are deallocated here.}
2317   \textcolor{keywordtype}{type}(set\_visc\_cs),   \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{   !< The control structure returned by a previous}
2318 \textcolor{comment}{                                             !! call to set\_visc\_init.}
2319   \textcolor{keywordflow}{if} (cs%bottomdraglaw) \textcolor{keywordflow}{then}
2320     \textcolor{keyword}{deallocate}(visc%bbl\_thick\_u) ; \textcolor{keyword}{deallocate}(visc%bbl\_thick\_v)
2321     \textcolor{keyword}{deallocate}(visc%kv\_bbl\_u) ; \textcolor{keyword}{deallocate}(visc%kv\_bbl\_v)
2322     \textcolor{keywordflow}{if} (\textcolor{keyword}{allocated}(cs%bbl\_u)) \textcolor{keyword}{deallocate}(cs%bbl\_u)
2323     \textcolor{keywordflow}{if} (\textcolor{keyword}{allocated}(cs%bbl\_v)) \textcolor{keyword}{deallocate}(cs%bbl\_v)
2324 \textcolor{keywordflow}{  endif}
2325   \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then}
2326     \textcolor{keyword}{deallocate}(visc%Ray\_u) ; \textcolor{keyword}{deallocate}(visc%Ray\_v)
2327 \textcolor{keywordflow}{  endif}
2328   \textcolor{keywordflow}{if} (cs%dynamic\_viscous\_ML) \textcolor{keywordflow}{then}
2329     \textcolor{keyword}{deallocate}(visc%nkml\_visc\_u) ; \textcolor{keyword}{deallocate}(visc%nkml\_visc\_v)
2330 \textcolor{keywordflow}{  endif}
2331   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kd\_shear)) \textcolor{keyword}{deallocate}(visc%Kd\_shear)
2332   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_slow)) \textcolor{keyword}{deallocate}(visc%Kv\_slow)
2333   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%TKE\_turb)) \textcolor{keyword}{deallocate}(visc%TKE\_turb)
2334   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_shear)) \textcolor{keyword}{deallocate}(visc%Kv\_shear)
2335   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_shear\_Bu)) \textcolor{keyword}{deallocate}(visc%Kv\_shear\_Bu)
2336   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%ustar\_bbl)) \textcolor{keyword}{deallocate}(visc%ustar\_bbl)
2337   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%TKE\_bbl)) \textcolor{keyword}{deallocate}(visc%TKE\_bbl)
2338   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%taux\_shelf)) \textcolor{keyword}{deallocate}(visc%taux\_shelf)
2339   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%tauy\_shelf)) \textcolor{keyword}{deallocate}(visc%tauy\_shelf)
2340   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%tbl\_thick\_shelf\_u)) \textcolor{keyword}{deallocate}(visc%tbl\_thick\_shelf\_u)
2341   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%tbl\_thick\_shelf\_v)) \textcolor{keyword}{deallocate}(visc%tbl\_thick\_shelf\_v)
2342   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%kv\_tbl\_shelf\_u)) \textcolor{keyword}{deallocate}(visc%kv\_tbl\_shelf\_u)
2343   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%kv\_tbl\_shelf\_v)) \textcolor{keyword}{deallocate}(visc%kv\_tbl\_shelf\_v)
2344 
2345   \textcolor{keyword}{deallocate}(cs)
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__set__visc_a2ec9be1a61c7d4a062aab4ff03e3ca29}\label{namespacemom__set__visc_a2ec9be1a61c7d4a062aab4ff03e3ca29}} 
\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}!set\+\_\+visc\+\_\+init@{set\+\_\+visc\+\_\+init}}
\index{set\+\_\+visc\+\_\+init@{set\+\_\+visc\+\_\+init}!mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}
\subsubsection{\texorpdfstring{set\+\_\+visc\+\_\+init()}{set\_visc\_init()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+set\+\_\+visc\+::set\+\_\+visc\+\_\+init (\begin{DoxyParamCaption}\item[{type(time\+\_\+type), intent(in), target}]{Time,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(inout)}]{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(vertvisc\+\_\+type), intent(inout)}]{visc,  }\item[{type(\hyperlink{structmom__set__visc_1_1set__visc__cs}{set\+\_\+visc\+\_\+cs}), pointer}]{CS,  }\item[{type(mom\+\_\+restart\+\_\+cs), pointer}]{restart\+\_\+\+CS,  }\item[{type(ocean\+\_\+obc\+\_\+type), pointer}]{O\+BC }\end{DoxyParamCaption})}



Initializes the M\+O\+M\+\_\+set\+\_\+visc control structure. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em time} & The current model time.\\
\hline
\mbox{\tt in,out}  & {\em g} & The ocean\textquotesingle{}s grid structure.\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure.\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em param\+\_\+file} & A structure to parse for run-\/time parameters.\\
\hline
\mbox{\tt in,out}  & {\em diag} & A structure that is used to regulate diagnostic output.\\
\hline
\mbox{\tt in,out}  & {\em visc} & A structure containing vertical viscosities and related fields. Allocated here.\\
\hline
 & {\em cs} & A pointer that is set to point to the control structure for this module\\
\hline
 & {\em restart\+\_\+cs} & A pointer to the restart control structure.\\
\hline
 & {\em obc} & A pointer to an open boundary condition structure \\
\hline
\end{DoxyParams}


Definition at line 1973 of file M\+O\+M\+\_\+set\+\_\+viscosity.\+F90.


\begin{DoxyCode}
1973   \textcolor{keywordtype}{type}(time\_type), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)}    :: time\textcolor{comment}{ !< The current model time.}
1974   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(inout)} :: g\textcolor{comment}{    !< The ocean's grid structure.}
1975   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure.}
1976   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{   !< A dimensional unit scaling type}
1977   \textcolor{keywordtype}{type}(param\_file\_type),   \textcolor{keywordtype}{intent(in)}    :: param\_file\textcolor{comment}{ !< A structure to parse for run-time}
1978 \textcolor{comment}{                                                 !! parameters.}
1979   \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(inout)} :: diag\textcolor{comment}{ !< A structure that is used to regulate diagnostic}
1980 \textcolor{comment}{                                                 !! output.}
1981   \textcolor{keywordtype}{type}(vertvisc\_type),     \textcolor{keywordtype}{intent(inout)} :: visc\textcolor{comment}{ !< A structure containing vertical viscosities and}
1982 \textcolor{comment}{                                                 !! related fields.  Allocated here.}
1983   \textcolor{keywordtype}{type}(set\_visc\_cs),       \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{   !< A pointer that is set to point to the control}
1984 \textcolor{comment}{                                                 !! structure for this module}
1985   \textcolor{keywordtype}{type}(mom\_restart\_cs),    \textcolor{keywordtype}{pointer}       :: restart\_cs\textcolor{comment}{ !< A pointer to the restart control structure.}
1986   \textcolor{keywordtype}{type}(ocean\_obc\_type),    \textcolor{keywordtype}{pointer}       :: obc\textcolor{comment}{  !< A pointer to an open boundary condition structure}
1987 
1988   \textcolor{comment}{! Local variables}
1989   \textcolor{keywordtype}{real}    :: csmag\_chan\_dflt, smag\_const1, tke\_decay\_dflt, bulk\_ri\_ml\_dflt
1990   \textcolor{keywordtype}{real}    :: kv\_background
1991   \textcolor{keywordtype}{real}    :: omega\_frac\_dflt
1992   \textcolor{keywordtype}{real}    :: z\_rescale     \textcolor{comment}{! A rescaling factor for heights from the representation in}
1993                            \textcolor{comment}{! a restart file to the internal representation in this run.}
1994   \textcolor{keywordtype}{real}    :: i\_t\_rescale   \textcolor{comment}{! A rescaling factor for time from the internal representation in this run}
1995                            \textcolor{comment}{! to the representation in a restart file.}
1996   \textcolor{keywordtype}{real}    :: z2\_t\_rescale  \textcolor{comment}{! A rescaling factor for vertical diffusivities and viscosities from the}
1997                            \textcolor{comment}{! representation in a restart file to the internal representation in this run.}
1998   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, n
1999   \textcolor{keywordtype}{integer} :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
2000   \textcolor{keywordtype}{logical} :: default\_2018\_answers
2001   \textcolor{keywordtype}{logical} :: use\_kappa\_shear, adiabatic, use\_omega, mle\_use\_pbl\_mld
2002   \textcolor{keywordtype}{logical} :: use\_kpp
2003   \textcolor{keywordtype}{logical} :: use\_regridding
2004   \textcolor{keywordtype}{character(len=200)} :: filename, tideamp\_file
2005   \textcolor{keywordtype}{type}(obc\_segment\_type), \textcolor{keywordtype}{pointer} :: segment => null() \textcolor{comment}{! pointer to OBC segment type}
2006   \textcolor{comment}{! This include declares and sets the variable "version".}
2007 \textcolor{preprocessor}{# include "version\_variable.h"}
2008 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_set\_visc"}  \textcolor{comment}{! This module's name.}
2009 
2010   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}
2011     \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"set\_visc\_init called with an associated "}// &
2012                             \textcolor{stringliteral}{"control structure."})
2013     \textcolor{keywordflow}{return}
2014 \textcolor{keywordflow}{  endif}
2015   \textcolor{keyword}{allocate}(cs)
2016 
2017   cs%OBC => obc
2018 
2019   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2020   isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
2021   isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2022 
2023   cs%diag => diag
2024 
2025   \textcolor{comment}{! Set default, read and log parameters}
2026   \textcolor{keyword}{call }log\_version(param\_file, mdl, version, \textcolor{stringliteral}{""})
2027   cs%RiNo\_mix = .false.
2028   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"INPUTDIR"}, cs%inputdir, default=\textcolor{stringliteral}{"."})
2029   cs%inputdir = slasher(cs%inputdir)
2030   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DEFAULT\_2018\_ANSWERS"}, default\_2018\_answers, &
2031                  \textcolor{stringliteral}{"This sets the default value for the various \_2018\_ANSWERS parameters."}, &
2032                  default=.false.)
2033   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"SET\_VISC\_2018\_ANSWERS"}, cs%answers\_2018, &
2034                  \textcolor{stringliteral}{"If true, use the order of arithmetic and expressions that recover the "}//&
2035                  \textcolor{stringliteral}{"answers from the end of 2018.  Otherwise, use updated and more robust "}//&
2036                  \textcolor{stringliteral}{"forms of the same expressions."}, default=default\_2018\_answers)
2037   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"BOTTOMDRAGLAW"}, cs%bottomdraglaw, &
2038                  \textcolor{stringliteral}{"If true, the bottom stress is calculated with a drag "}//&
2039                  \textcolor{stringliteral}{"law of the form c\_drag*|u|*u. The velocity magnitude "}//&
2040                  \textcolor{stringliteral}{"may be an assumed value or it may be based on the "}//&
2041                  \textcolor{stringliteral}{"actual velocity in the bottommost HBBL, depending on "}//&
2042                  \textcolor{stringliteral}{"LINEAR\_DRAG."}, default=.true.)
2043   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHANNEL\_DRAG"}, cs%Channel\_drag, &
2044                  \textcolor{stringliteral}{"If true, the bottom drag is exerted directly on each "}//&
2045                  \textcolor{stringliteral}{"layer proportional to the fraction of the bottom it "}//&
2046                  \textcolor{stringliteral}{"overlies."}, default=.false.)
2047   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"LINEAR\_DRAG"}, cs%linear\_drag, &
2048                  \textcolor{stringliteral}{"If LINEAR\_DRAG and BOTTOMDRAGLAW are defined the drag "}//&
2049                  \textcolor{stringliteral}{"law is cdrag*DRAG\_BG\_VEL*u."}, default=.false.)
2050   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"ADIABATIC"}, adiabatic, default=.false., &
2051                  do\_not\_log=.true.)
2052   \textcolor{keywordflow}{if} (adiabatic) \textcolor{keywordflow}{then}
2053     \textcolor{keyword}{call }log\_param(param\_file, mdl, \textcolor{stringliteral}{"ADIABATIC"},adiabatic, &
2054                  \textcolor{stringliteral}{"There are no diapycnal mass fluxes if ADIABATIC is "}//&
2055                  \textcolor{stringliteral}{"true. This assumes that KD = KDML = 0.0 and that "}//&
2056                  \textcolor{stringliteral}{"there is no buoyancy forcing, but makes the model "}//&
2057                  \textcolor{stringliteral}{"faster by eliminating subroutine calls."}, default=.false.)
2058 \textcolor{keywordflow}{  endif}
2059 
2060   \textcolor{keywordflow}{if} (.not.adiabatic) \textcolor{keywordflow}{then}
2061     cs%RiNo\_mix = kappa\_shear\_is\_used(param\_file)
2062 \textcolor{keywordflow}{  endif}
2063 
2064   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"PRANDTL\_TURB"}, visc%Prandtl\_turb, &
2065                  \textcolor{stringliteral}{"The turbulent Prandtl number applied to shear "}//&
2066                  \textcolor{stringliteral}{"instability."}, units=\textcolor{stringliteral}{"nondim"}, default=1.0)
2067   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DEBUG"}, cs%debug, default=.false.)
2068 
2069   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DYNAMIC\_VISCOUS\_ML"}, cs%dynamic\_viscous\_ML, &
2070                  \textcolor{stringliteral}{"If true, use a bulk Richardson number criterion to "}//&
2071                  \textcolor{stringliteral}{"determine the mixed layer thickness for viscosity."}, &
2072                  default=.false.)
2073   \textcolor{keywordflow}{if} (cs%dynamic\_viscous\_ML) \textcolor{keywordflow}{then}
2074     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"BULK\_RI\_ML"}, bulk\_ri\_ml\_dflt, default=0.0)
2075     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"BULK\_RI\_ML\_VISC"}, cs%bulk\_Ri\_ML, &
2076                  \textcolor{stringliteral}{"The efficiency with which mean kinetic energy released "}//&
2077                  \textcolor{stringliteral}{"by mechanically forced entrainment of the mixed layer "}//&
2078                  \textcolor{stringliteral}{"is converted to turbulent kinetic energy.  By default, "}//&
2079                  \textcolor{stringliteral}{"BULK\_RI\_ML\_VISC = BULK\_RI\_ML or 0."}, units=\textcolor{stringliteral}{"nondim"}, &
2080                  default=bulk\_ri\_ml\_dflt)
2081     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"TKE\_DECAY"}, tke\_decay\_dflt, default=0.0)
2082     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"TKE\_DECAY\_VISC"}, cs%TKE\_decay, &
2083                  \textcolor{stringliteral}{"TKE\_DECAY\_VISC relates the vertical rate of decay of "}//&
2084                  \textcolor{stringliteral}{"the TKE available for mechanical entrainment to the "}//&
2085                  \textcolor{stringliteral}{"natural Ekman depth for use in calculating the dynamic "}//&
2086                  \textcolor{stringliteral}{"mixed layer viscosity.  By default, "}//&
2087                  \textcolor{stringliteral}{"TKE\_DECAY\_VISC = TKE\_DECAY or 0."}, units=\textcolor{stringliteral}{"nondim"}, &
2088                  default=tke\_decay\_dflt)
2089     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"ML\_USE\_OMEGA"}, use\_omega, &
2090                  \textcolor{stringliteral}{"If true, use the absolute rotation rate instead of the "}//&
2091                  \textcolor{stringliteral}{"vertical component of rotation when setting the decay "}//&
2092                    \textcolor{stringliteral}{"scale for turbulence."}, default=.false., do\_not\_log=.true.)
2093     omega\_frac\_dflt = 0.0
2094     \textcolor{keywordflow}{if} (use\_omega) \textcolor{keywordflow}{then}
2095       \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"ML\_USE\_OMEGA is depricated; use ML\_OMEGA\_FRAC=1.0 instead."})
2096       omega\_frac\_dflt = 1.0
2097 \textcolor{keywordflow}{    endif}
2098     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"ML\_OMEGA\_FRAC"}, cs%omega\_frac, &
2099                    \textcolor{stringliteral}{"When setting the decay scale for turbulence, use this "}//&
2100                    \textcolor{stringliteral}{"fraction of the absolute rotation rate blended with the "}//&
2101                    \textcolor{stringliteral}{"local value of f, as sqrt((1-of)*f^2 + of*4*omega^2)."}, &
2102                    units=\textcolor{stringliteral}{"nondim"}, default=omega\_frac\_dflt)
2103     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"OMEGA"}, cs%omega, &
2104                  \textcolor{stringliteral}{"The rotation rate of the earth."}, units=\textcolor{stringliteral}{"s-1"}, &
2105                  default=7.2921e-5, scale=us%T\_to\_s)
2106     \textcolor{comment}{! This give a minimum decay scale that is typically much less than Angstrom.}
2107     cs%ustar\_min = 2e-4*cs%omega*(gv%Angstrom\_Z + gv%H\_to\_Z*gv%H\_subroundoff)
2108   \textcolor{keywordflow}{else}
2109     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"OMEGA"}, cs%omega, &
2110                  \textcolor{stringliteral}{"The rotation rate of the earth."}, units=\textcolor{stringliteral}{"s-1"}, &
2111                  default=7.2921e-5, scale=us%T\_to\_s)
2112 \textcolor{keywordflow}{  endif}
2113 
2114   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HBBL"}, cs%Hbbl, &
2115                  \textcolor{stringliteral}{"The thickness of a bottom boundary layer with a "}//&
2116                  \textcolor{stringliteral}{"viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "}//&
2117                  \textcolor{stringliteral}{"the thickness over which near-bottom velocities are "}//&
2118                  \textcolor{stringliteral}{"averaged for the drag law if BOTTOMDRAGLAW is defined "}//&
2119                  \textcolor{stringliteral}{"but LINEAR\_DRAG is not."}, units=\textcolor{stringliteral}{"m"}, fail\_if\_missing=.true.) \textcolor{comment}{! Rescaled later}
2120   \textcolor{keywordflow}{if} (cs%bottomdraglaw) \textcolor{keywordflow}{then}
2121     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CDRAG"}, cs%cdrag, &
2122                  \textcolor{stringliteral}{"CDRAG is the drag coefficient relating the magnitude of "}//&
2123                  \textcolor{stringliteral}{"the velocity field to the bottom stress. CDRAG is only "}//&
2124                  \textcolor{stringliteral}{"used if BOTTOMDRAGLAW is defined."}, units=\textcolor{stringliteral}{"nondim"}, &
2125                  default=0.003)
2126     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"BBL\_USE\_TIDAL\_BG"}, cs%BBL\_use\_tidal\_bg, &
2127                  \textcolor{stringliteral}{"Flag to use the tidal RMS amplitude in place of constant "}//&
2128                  \textcolor{stringliteral}{"background velocity for computing u* in the BBL. "}//&
2129                  \textcolor{stringliteral}{"This flag is only used when BOTTOMDRAGLAW is true and "}//&
2130                  \textcolor{stringliteral}{"LINEAR\_DRAG is false."}, default=.false.)
2131     \textcolor{keywordflow}{if} (cs%BBL\_use\_tidal\_bg) \textcolor{keywordflow}{then}
2132       \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"TIDEAMP\_FILE"}, tideamp\_file, &
2133                    \textcolor{stringliteral}{"The path to the file containing the spatially varying "}//&
2134                    \textcolor{stringliteral}{"tidal amplitudes with INT\_TIDE\_DISSIPATION."}, default=\textcolor{stringliteral}{"tideamp.nc"})
2135     \textcolor{keywordflow}{else}
2136       \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DRAG\_BG\_VEL"}, cs%drag\_bg\_vel, &
2137                    \textcolor{stringliteral}{"DRAG\_BG\_VEL is either the assumed bottom velocity (with "}//&
2138                    \textcolor{stringliteral}{"LINEAR\_DRAG) or an unresolved  velocity that is "}//&
2139                    \textcolor{stringliteral}{"combined with the resolved velocity to estimate the "}//&
2140                    \textcolor{stringliteral}{"velocity magnitude.  DRAG\_BG\_VEL is only used when "}//&
2141                    \textcolor{stringliteral}{"BOTTOMDRAGLAW is defined."}, units=\textcolor{stringliteral}{"m s-1"}, default=0.0, scale=us%m\_s\_to\_L\_T)
2142 \textcolor{keywordflow}{    endif}
2143     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_REGRIDDING"}, use\_regridding, &
2144          do\_not\_log = .true., default = .false. )
2145     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"BBL\_USE\_EOS"}, cs%BBL\_use\_EOS, &
2146                  \textcolor{stringliteral}{"If true, use the equation of state in determining the "}//&
2147                  \textcolor{stringliteral}{"properties of the bottom boundary layer.  Otherwise use "}//&
2148                  \textcolor{stringliteral}{"the layer target potential densities.  The default of "}//&
2149                  \textcolor{stringliteral}{"this is determined by USE\_REGRIDDING."}, default=use\_regridding)
2150     \textcolor{keywordflow}{if} (use\_regridding .and. (.not. cs%BBL\_use\_EOS)) &
2151       \textcolor{keyword}{call }mom\_error(fatal,\textcolor{stringliteral}{"When using MOM6 in ALE mode it is required to "}//&
2152            \textcolor{stringliteral}{"set BBL\_USE\_EOS to True"})
2153 \textcolor{keywordflow}{  endif}
2154   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"BBL\_THICK\_MIN"}, cs%BBL\_thick\_min, &
2155                  \textcolor{stringliteral}{"The minimum bottom boundary layer thickness that can be "}//&
2156                  \textcolor{stringliteral}{"used with BOTTOMDRAGLAW. This might be "}//&
2157                  \textcolor{stringliteral}{"Kv/(cdrag*drag\_bg\_vel) to give Kv as the minimum "}//&
2158                  \textcolor{stringliteral}{"near-bottom viscosity."}, units=\textcolor{stringliteral}{"m"}, default=0.0)  \textcolor{comment}{! Rescaled later}
2159   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HTBL\_SHELF\_MIN"}, cs%Htbl\_shelf\_min, &
2160                  \textcolor{stringliteral}{"The minimum top boundary layer thickness that can be "}//&
2161                  \textcolor{stringliteral}{"used with BOTTOMDRAGLAW. This might be "}//&
2162                  \textcolor{stringliteral}{"Kv/(cdrag*drag\_bg\_vel) to give Kv as the minimum "}//&
2163                  \textcolor{stringliteral}{"near-top viscosity."}, units=\textcolor{stringliteral}{"m"}, default=cs%BBL\_thick\_min, scale=gv%m\_to\_H)
2164   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HTBL\_SHELF"}, cs%Htbl\_shelf, &
2165                  \textcolor{stringliteral}{"The thickness over which near-surface velocities are "}//&
2166                  \textcolor{stringliteral}{"averaged for the drag law under an ice shelf.  By "}//&
2167                  \textcolor{stringliteral}{"default this is the same as HBBL"}, units=\textcolor{stringliteral}{"m"}, default=cs%Hbbl, scale=gv%m\_to\_H)
2168   \textcolor{comment}{! These unit conversions are out outside the get\_param calls because the are also defaults.}
2169   cs%Hbbl = cs%Hbbl * gv%m\_to\_H                   \textcolor{comment}{! Rescale}
2170   cs%BBL\_thick\_min = cs%BBL\_thick\_min * gv%m\_to\_H \textcolor{comment}{! Rescale}
2171 
2172   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KV"}, kv\_background, &
2173                  \textcolor{stringliteral}{"The background kinematic viscosity in the interior. "}//&
2174                  \textcolor{stringliteral}{"The molecular value, ~1e-6 m2 s-1, may be used."}, &
2175                  units=\textcolor{stringliteral}{"m2 s-1"}, fail\_if\_missing=.true.)
2176 
2177   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_KPP"}, use\_kpp, &
2178                  \textcolor{stringliteral}{"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "}//&
2179                  \textcolor{stringliteral}{"to calculate diffusivities and non-local transport in the OBL."}, &
2180                  do\_not\_log=.true., default=.false.)
2181 
2182   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KV\_BBL\_MIN"}, cs%KV\_BBL\_min, &
2183                  \textcolor{stringliteral}{"The minimum viscosities in the bottom boundary layer."}, &
2184                  units=\textcolor{stringliteral}{"m2 s-1"}, default=kv\_background, scale=us%m2\_s\_to\_Z2\_T)
2185   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KV\_TBL\_MIN"}, cs%KV\_TBL\_min, &
2186                  \textcolor{stringliteral}{"The minimum viscosities in the top boundary layer."}, &
2187                  units=\textcolor{stringliteral}{"m2 s-1"}, default=kv\_background, scale=us%m2\_s\_to\_Z2\_T)
2188   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CORRECT\_BBL\_BOUNDS"}, cs%correct\_BBL\_bounds, &
2189                  \textcolor{stringliteral}{"If true, uses the correct bounds on the BBL thickness and "}//&
2190                  \textcolor{stringliteral}{"viscosity so that the bottom layer feels the intended drag."}, &
2191                  default=.false.)
2192 
2193   \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then}
2194     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"SMAG\_LAP\_CONST"}, smag\_const1, default=-1.0)
2195 
2196     csmag\_chan\_dflt = 0.15
2197     \textcolor{keywordflow}{if} (smag\_const1 >= 0.0) csmag\_chan\_dflt = smag\_const1
2198 
2199     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"SMAG\_CONST\_CHANNEL"}, cs%c\_Smag, &
2200                  \textcolor{stringliteral}{"The nondimensional Laplacian Smagorinsky constant used "}//&
2201                  \textcolor{stringliteral}{"in calculating the channel drag if it is enabled.  The "}//&
2202                  \textcolor{stringliteral}{"default is to use the same value as SMAG\_LAP\_CONST if "}//&
2203                  \textcolor{stringliteral}{"it is defined, or 0.15 if it is not. The value used is "}//&
2204                  \textcolor{stringliteral}{"also 0.15 if the specified value is negative."}, &
2205                  units=\textcolor{stringliteral}{"nondim"}, default=csmag\_chan\_dflt)
2206     \textcolor{keywordflow}{if} (cs%c\_Smag < 0.0) cs%c\_Smag = 0.15
2207 \textcolor{keywordflow}{  endif}
2208 
2209   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MLE\_USE\_PBL\_MLD"}, mle\_use\_pbl\_mld, &
2210                  default=.false., do\_not\_log=.true.)
2211 
2212   \textcolor{keywordflow}{if} (cs%RiNo\_mix .and. kappa\_shear\_at\_vertex(param\_file)) \textcolor{keywordflow}{then}
2213     \textcolor{comment}{! This is necessary for reproduciblity across restarts in non-symmetric mode.}
2214     \textcolor{keyword}{call }pass\_var(visc%Kv\_shear\_Bu, g%Domain, position=corner, complete=.true.)
2215 \textcolor{keywordflow}{  endif}
2216 
2217   \textcolor{keywordflow}{if} (cs%bottomdraglaw) \textcolor{keywordflow}{then}
2218     \textcolor{keyword}{allocate}(visc%bbl\_thick\_u(isdb:iedb,jsd:jed)) ; visc%bbl\_thick\_u(:,:) = 0.0
2219     \textcolor{keyword}{allocate}(visc%kv\_bbl\_u(isdb:iedb,jsd:jed)) ; visc%kv\_bbl\_u(:,:) = 0.0
2220     \textcolor{keyword}{allocate}(visc%bbl\_thick\_v(isd:ied,jsdb:jedb)) ; visc%bbl\_thick\_v(:,:) = 0.0
2221     \textcolor{keyword}{allocate}(visc%kv\_bbl\_v(isd:ied,jsdb:jedb)) ; visc%kv\_bbl\_v(:,:) = 0.0
2222     \textcolor{keyword}{allocate}(visc%ustar\_bbl(isd:ied,jsd:jed)) ; visc%ustar\_bbl(:,:) = 0.0
2223     \textcolor{keyword}{allocate}(visc%TKE\_bbl(isd:ied,jsd:jed)) ; visc%TKE\_bbl(:,:) = 0.0
2224 
2225     cs%id\_bbl\_thick\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'bbl\_thick\_u'}, &
2226        diag%axesCu1, time, \textcolor{stringliteral}{'BBL thickness at u points'}, \textcolor{stringliteral}{'m'}, conversion=us%Z\_to\_m)
2227     cs%id\_kv\_bbl\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'kv\_bbl\_u'}, diag%axesCu1, &
2228        time, \textcolor{stringliteral}{'BBL viscosity at u points'}, \textcolor{stringliteral}{'m2 s-1'}, conversion=us%Z2\_T\_to\_m2\_s)
2229     cs%id\_bbl\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'bbl\_u'}, diag%axesCu1, &
2230        time, \textcolor{stringliteral}{'BBL mean u current'}, \textcolor{stringliteral}{'m s-1'}, conversion=us%L\_T\_to\_m\_s)
2231     \textcolor{keywordflow}{if} (cs%id\_bbl\_u>0) \textcolor{keywordflow}{then}
2232       \textcolor{keyword}{allocate}(cs%bbl\_u(isdb:iedb,jsd:jed)) ; cs%bbl\_u(:,:) = 0.0
2233 \textcolor{keywordflow}{    endif}
2234     cs%id\_bbl\_thick\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'bbl\_thick\_v'}, &
2235        diag%axesCv1, time, \textcolor{stringliteral}{'BBL thickness at v points'}, \textcolor{stringliteral}{'m'}, conversion=us%Z\_to\_m)
2236     cs%id\_kv\_bbl\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'kv\_bbl\_v'}, diag%axesCv1, &
2237        time, \textcolor{stringliteral}{'BBL viscosity at v points'}, \textcolor{stringliteral}{'m2 s-1'}, conversion=us%Z2\_T\_to\_m2\_s)
2238     cs%id\_bbl\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'bbl\_v'}, diag%axesCv1, &
2239        time, \textcolor{stringliteral}{'BBL mean v current'}, \textcolor{stringliteral}{'m s-1'}, conversion=us%L\_T\_to\_m\_s)
2240     \textcolor{keywordflow}{if} (cs%id\_bbl\_v>0) \textcolor{keywordflow}{then}
2241       \textcolor{keyword}{allocate}(cs%bbl\_v(isd:ied,jsdb:jedb)) ; cs%bbl\_v(:,:) = 0.0
2242 \textcolor{keywordflow}{    endif}
2243     \textcolor{keywordflow}{if} (cs%BBL\_use\_tidal\_bg) \textcolor{keywordflow}{then}
2244       \textcolor{keyword}{allocate}(cs%tideamp(isd:ied,jsd:jed)) ; cs%tideamp(:,:) = 0.0
2245       filename = trim(cs%inputdir) // trim(tideamp\_file)
2246       \textcolor{keyword}{call }log\_param(param\_file, mdl, \textcolor{stringliteral}{"INPUTDIR/TIDEAMP\_FILE"}, filename)
2247       \textcolor{keyword}{call }mom\_read\_data(filename, \textcolor{stringliteral}{'tideamp'}, cs%tideamp, g%domain, timelevel=1, scale=us%m\_to\_Z*us%T\_to\_s)
2248       \textcolor{keyword}{call }pass\_var(cs%tideamp,g%domain)
2249 \textcolor{keywordflow}{    endif}
2250 \textcolor{keywordflow}{  endif}
2251   \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then}
2252     \textcolor{keyword}{allocate}(visc%Ray\_u(isdb:iedb,jsd:jed,nz)) ; visc%Ray\_u(:,:,:) = 0.0
2253     \textcolor{keyword}{allocate}(visc%Ray\_v(isd:ied,jsdb:jedb,nz)) ; visc%Ray\_v(:,:,:) = 0.0
2254     cs%id\_Ray\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Rayleigh\_u'}, diag%axesCuL, &
2255        time, \textcolor{stringliteral}{'Rayleigh drag velocity at u points'}, \textcolor{stringliteral}{'m s-1'}, conversion=us%Z\_to\_m*us%s\_to\_T)
2256     cs%id\_Ray\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Rayleigh\_v'}, diag%axesCvL, &
2257        time, \textcolor{stringliteral}{'Rayleigh drag velocity at v points'}, \textcolor{stringliteral}{'m s-1'}, conversion=us%Z\_to\_m*us%s\_to\_T)
2258 \textcolor{keywordflow}{  endif}
2259 
2260 
2261   \textcolor{keywordflow}{if} (cs%dynamic\_viscous\_ML) \textcolor{keywordflow}{then}
2262     \textcolor{keyword}{allocate}(visc%nkml\_visc\_u(isdb:iedb,jsd:jed)) ; visc%nkml\_visc\_u(:,:) = 0.0
2263     \textcolor{keyword}{allocate}(visc%nkml\_visc\_v(isd:ied,jsdb:jedb)) ; visc%nkml\_visc\_v(:,:) = 0.0
2264     cs%id\_nkml\_visc\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'nkml\_visc\_u'}, &
2265        diag%axesCu1, time, \textcolor{stringliteral}{'Number of layers in viscous mixed layer at u points'}, \textcolor{stringliteral}{'m'})
2266     cs%id\_nkml\_visc\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'nkml\_visc\_v'}, &
2267        diag%axesCv1, time, \textcolor{stringliteral}{'Number of layers in viscous mixed layer at v points'}, \textcolor{stringliteral}{'m'})
2268 \textcolor{keywordflow}{  endif}
2269 
2270   \textcolor{keyword}{call }register\_restart\_field\_as\_obsolete(\textcolor{stringliteral}{'Kd\_turb'},\textcolor{stringliteral}{'Kd\_shear'}, restart\_cs)
2271   \textcolor{keyword}{call }register\_restart\_field\_as\_obsolete(\textcolor{stringliteral}{'Kv\_turb'},\textcolor{stringliteral}{'Kv\_shear'}, restart\_cs)
2272 
2273   \textcolor{comment}{! Account for possible changes in dimensional scaling for variables that have been}
2274   \textcolor{comment}{! read from a restart file.}
2275   z\_rescale = 1.0
2276   \textcolor{keywordflow}{if} ((us%m\_to\_Z\_restart /= 0.0) .and. (us%m\_to\_Z\_restart /= us%m\_to\_Z)) &
2277     z\_rescale = us%m\_to\_Z / us%m\_to\_Z\_restart
2278   i\_t\_rescale = 1.0
2279   \textcolor{keywordflow}{if} ((us%s\_to\_T\_restart /= 0.0) .and. (us%s\_to\_T\_restart /= us%s\_to\_T)) &
2280     i\_t\_rescale = us%s\_to\_T\_restart / us%s\_to\_T
2281   z2\_t\_rescale = z\_rescale**2*i\_t\_rescale
2282 
2283   \textcolor{keywordflow}{if} (z2\_t\_rescale /= 1.0) \textcolor{keywordflow}{then}
2284     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kd\_shear)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (query\_initialized(visc%Kd\_shear, \textcolor{stringliteral}{"Kd\_shear"}, restart\_cs)) \textcolor{keywordflow}{
      then}
2285       \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
2286         visc%Kd\_shear(i,j,k) = z2\_t\_rescale * visc%Kd\_shear(i,j,k)
2287 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
2288 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
2289 
2290     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_shear)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (query\_initialized(visc%Kv\_shear, \textcolor{stringliteral}{"Kv\_shear"}, restart\_cs)) \textcolor{keywordflow}{
      then}
2291       \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
2292         visc%Kv\_shear(i,j,k) = z2\_t\_rescale * visc%Kv\_shear(i,j,k)
2293 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
2294 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
2295 
2296     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_shear\_Bu)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (query\_initialized(visc%Kv\_shear\_Bu, \textcolor{stringliteral}{"Kv\_shear\_Bu"}, 
      restart\_cs)) \textcolor{keywordflow}{then}
2297       \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is-1,ie
2298         visc%Kv\_shear\_Bu(i,j,k) = z2\_t\_rescale * visc%Kv\_shear\_Bu(i,j,k)
2299 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
2300 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
2301 \textcolor{keywordflow}{  endif}
2302 
2303   \textcolor{keywordflow}{if} (mle\_use\_pbl\_mld .and. (z\_rescale /= 1.0)) \textcolor{keywordflow}{then}
2304     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%MLD)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (query\_initialized(visc%MLD, \textcolor{stringliteral}{"MLD"}, restart\_cs)) \textcolor{keywordflow}{then}
2305       \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=is,ie
2306         visc%MLD(i,j) = z\_rescale * visc%MLD(i,j)
2307 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
2308 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
2309 \textcolor{keywordflow}{  endif}
2310 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__set__visc_ae2d9d9f74c1e9aec56257cfad372b0fd}\label{namespacemom__set__visc_ae2d9d9f74c1e9aec56257cfad372b0fd}} 
\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}!set\+\_\+visc\+\_\+register\+\_\+restarts@{set\+\_\+visc\+\_\+register\+\_\+restarts}}
\index{set\+\_\+visc\+\_\+register\+\_\+restarts@{set\+\_\+visc\+\_\+register\+\_\+restarts}!mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}
\subsubsection{\texorpdfstring{set\+\_\+visc\+\_\+register\+\_\+restarts()}{set\_visc\_register\_restarts()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+set\+\_\+visc\+::set\+\_\+visc\+\_\+register\+\_\+restarts (\begin{DoxyParamCaption}\item[{type(hor\+\_\+index\+\_\+type), intent(in)}]{HI,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(param\+\_\+file\+\_\+type), intent(in)}]{param\+\_\+file,  }\item[{type(vertvisc\+\_\+type), intent(inout)}]{visc,  }\item[{type(mom\+\_\+restart\+\_\+cs), pointer}]{restart\+\_\+\+CS }\end{DoxyParamCaption})}



Register any fields associated with the vertvisc\+\_\+type. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em hi} & A horizontal index type structure.\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure.\\
\hline
\mbox{\tt in}  & {\em param\+\_\+file} & A structure to parse for run-\/time parameters.\\
\hline
\mbox{\tt in,out}  & {\em visc} & A structure containing vertical viscosities and related fields. Allocated here.\\
\hline
 & {\em restart\+\_\+cs} & A pointer to the restart control structure. \\
\hline
\end{DoxyParams}


Definition at line 1887 of file M\+O\+M\+\_\+set\+\_\+viscosity.\+F90.


\begin{DoxyCode}
1887   \textcolor{keywordtype}{type}(hor\_index\_type),    \textcolor{keywordtype}{intent(in)}    :: hi\textcolor{comment}{         !< A horizontal index type structure.}
1888   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{         !< The ocean's vertical grid structure.}
1889   \textcolor{keywordtype}{type}(param\_file\_type),   \textcolor{keywordtype}{intent(in)}    :: param\_file\textcolor{comment}{ !< A structure to parse for run-time}
1890 \textcolor{comment}{                                                       !! parameters.}
1891   \textcolor{keywordtype}{type}(vertvisc\_type),     \textcolor{keywordtype}{intent(inout)} :: visc\textcolor{comment}{       !< A structure containing vertical}
1892 \textcolor{comment}{                                                       !! viscosities and related fields.}
1893 \textcolor{comment}{                                                       !! Allocated here.}
1894   \textcolor{keywordtype}{type}(mom\_restart\_cs),    \textcolor{keywordtype}{pointer}       :: restart\_cs\textcolor{comment}{ !< A pointer to the restart control structure.}
1895   \textcolor{comment}{! Local variables}
1896   \textcolor{keywordtype}{logical} :: use\_kappa\_shear, ks\_at\_vertex
1897   \textcolor{keywordtype}{logical} :: adiabatic, usekpp, useepbl
1898   \textcolor{keywordtype}{logical} :: use\_cvmix\_shear, mle\_use\_pbl\_mld, use\_cvmix\_conv
1899   \textcolor{keywordtype}{integer} :: isd, ied, jsd, jed, nz
1900   \textcolor{keywordtype}{real} :: hfreeze\textcolor{comment}{ !< If hfreeze > 0 [m], melt potential will be computed.}
1901   \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_set\_visc"}  \textcolor{comment}{! This module's name.}
1902   isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1903 
1904   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"ADIABATIC"}, adiabatic, default=.false., &
1905                  do\_not\_log=.true.)
1906 
1907   use\_kappa\_shear = .false. ; ks\_at\_vertex = .false. ; use\_cvmix\_shear = .false.
1908   usekpp = .false. ; useepbl = .false. ; use\_cvmix\_conv = .false.
1909 
1910   \textcolor{keywordflow}{if} (.not.adiabatic) \textcolor{keywordflow}{then}
1911     use\_kappa\_shear = kappa\_shear\_is\_used(param\_file)
1912     ks\_at\_vertex    = kappa\_shear\_at\_vertex(param\_file)
1913     use\_cvmix\_shear = cvmix\_shear\_is\_used(param\_file)
1914     use\_cvmix\_conv  = cvmix\_conv\_is\_used(param\_file)
1915     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"USE\_KPP"}, usekpp, &
1916                  \textcolor{stringliteral}{"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "}//&
1917                  \textcolor{stringliteral}{"to calculate diffusivities and non-local transport in the OBL."}, &
1918                  default=.false., do\_not\_log=.true.)
1919     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"ENERGETICS\_SFC\_PBL"}, useepbl, &
1920                  \textcolor{stringliteral}{"If true, use an implied energetics planetary boundary "}//&
1921                  \textcolor{stringliteral}{"layer scheme to determine the diffusivity and viscosity "}//&
1922                  \textcolor{stringliteral}{"in the surface boundary layer."}, default=.false., do\_not\_log=.true.)
1923 \textcolor{keywordflow}{  endif}
1924 
1925   \textcolor{keywordflow}{if} (use\_kappa\_shear .or. usekpp .or. useepbl .or. use\_cvmix\_shear .or. use\_cvmix\_conv) \textcolor{keywordflow}{then}
1926     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%Kd\_shear, isd, ied, jsd, jed, nz+1)
1927     \textcolor{keyword}{call }register\_restart\_field(visc%Kd\_shear, \textcolor{stringliteral}{"Kd\_shear"}, .false., restart\_cs, &
1928                   \textcolor{stringliteral}{"Shear-driven turbulent diffusivity at interfaces"}, \textcolor{stringliteral}{"m2 s-1"}, z\_grid=\textcolor{stringliteral}{'i'})
1929 \textcolor{keywordflow}{  endif}
1930   \textcolor{keywordflow}{if} (usekpp .or. useepbl .or. use\_cvmix\_shear .or. use\_cvmix\_conv .or. &
1931       (use\_kappa\_shear .and. .not.ks\_at\_vertex )) \textcolor{keywordflow}{then}
1932     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%Kv\_shear, isd, ied, jsd, jed, nz+1)
1933     \textcolor{keyword}{call }register\_restart\_field(visc%Kv\_shear, \textcolor{stringliteral}{"Kv\_shear"}, .false., restart\_cs, &
1934                   \textcolor{stringliteral}{"Shear-driven turbulent viscosity at interfaces"}, \textcolor{stringliteral}{"m2 s-1"}, z\_grid=\textcolor{stringliteral}{'i'})
1935 \textcolor{keywordflow}{  endif}
1936   \textcolor{keywordflow}{if} (use\_kappa\_shear .and. ks\_at\_vertex) \textcolor{keywordflow}{then}
1937     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%TKE\_turb, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1938     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%Kv\_shear\_Bu, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1939     \textcolor{keyword}{call }register\_restart\_field(visc%Kv\_shear\_Bu, \textcolor{stringliteral}{"Kv\_shear\_Bu"}, .false., restart\_cs, &
1940                   \textcolor{stringliteral}{"Shear-driven turbulent viscosity at vertex interfaces"}, \textcolor{stringliteral}{"m2 s-1"}, &
1941                   hor\_grid=\textcolor{stringliteral}{"Bu"}, z\_grid=\textcolor{stringliteral}{'i'})
1942   \textcolor{keywordflow}{elseif} (use\_kappa\_shear) \textcolor{keywordflow}{then}
1943     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%TKE\_turb, isd, ied, jsd, jed, nz+1)
1944 \textcolor{keywordflow}{  endif}
1945 
1946   \textcolor{keywordflow}{if} (usekpp) \textcolor{keywordflow}{then}
1947     \textcolor{comment}{! MOM\_bkgnd\_mixing uses Kv\_slow when KPP is defined.}
1948     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%Kv\_slow, isd, ied, jsd, jed, nz+1)
1949 \textcolor{keywordflow}{  endif}
1950 
1951   \textcolor{comment}{! visc%MLD is used to communicate the state of the (e)PBL or KPP to the rest of the model}
1952   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MLE\_USE\_PBL\_MLD"}, mle\_use\_pbl\_mld, &
1953                  default=.false., do\_not\_log=.true.)
1954   \textcolor{comment}{! visc%MLD needs to be allocated when melt potential is computed (HFREEZE>0)}
1955   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HFREEZE"}, hfreeze, &
1956                  default=-1.0, do\_not\_log=.true.)
1957 
1958   \textcolor{keywordflow}{if} (mle\_use\_pbl\_mld) \textcolor{keywordflow}{then}
1959     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%MLD, isd, ied, jsd, jed)
1960     \textcolor{keyword}{call }register\_restart\_field(visc%MLD, \textcolor{stringliteral}{"MLD"}, .false., restart\_cs, &
1961                   \textcolor{stringliteral}{"Instantaneous active mixing layer depth"}, \textcolor{stringliteral}{"m"})
1962 \textcolor{keywordflow}{  endif}
1963 
1964   \textcolor{keywordflow}{if} (hfreeze >= 0.0 .and. .not.mle\_use\_pbl\_mld) \textcolor{keywordflow}{then}
1965     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%MLD, isd, ied, jsd, jed)
1966 \textcolor{keywordflow}{  endif}
1967 
1968 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__set__visc_a9865113fe07928e7a240c2868ed45e5f}\label{namespacemom__set__visc_a9865113fe07928e7a240c2868ed45e5f}} 
\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}!set\+\_\+viscous\+\_\+bbl@{set\+\_\+viscous\+\_\+bbl}}
\index{set\+\_\+viscous\+\_\+bbl@{set\+\_\+viscous\+\_\+bbl}!mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}
\subsubsection{\texorpdfstring{set\+\_\+viscous\+\_\+bbl()}{set\_viscous\_bbl()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+set\+\_\+visc\+::set\+\_\+viscous\+\_\+bbl (\begin{DoxyParamCaption}\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(g)), intent(in)}]{v,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(g)), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{type(vertvisc\+\_\+type), intent(inout)}]{visc,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(inout)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\hyperlink{structmom__set__visc_1_1set__visc__cs}{set\+\_\+visc\+\_\+cs}), pointer}]{CS,  }\item[{logical, intent(in), optional}]{symmetrize }\end{DoxyParamCaption})}



Calculates the thickness of the bottom boundary layer and the viscosity within that layer. 

A drag law is used, either linearized about an assumed bottom velocity or using the actual near-\/bottom velocities combined with an assumed unresolved velocity. The bottom boundary layer thickness is limited by a combination of stratification and rotation, as in the paper of Killworth and Edwards, J\+PO 1999. It is not necessary to calculate the thickness and viscosity every time step; instead previous values may be used.\hypertarget{namespacemom__set__visc_set_viscous_BBL}{}\subsection{Viscous Bottom Boundary Layer}\label{namespacemom__set__visc_set_viscous_BBL}
If \hyperlink{structmom__set__visc_1_1set__visc__cs_aca87adb9636e724f7ccb78e7c640f399}{set\+\_\+visc\+\_\+cs.\+bottomdraglaw} is True then a bottom boundary layer viscosity and thickness are calculated so that the bottom stress is \[ \mathbf{\tau}_b = C_d | U_{bbl} | \mathbf{u}_{bbl} \] If \hyperlink{structmom__set__visc_1_1set__visc__cs_aca87adb9636e724f7ccb78e7c640f399}{set\+\_\+visc\+\_\+cs.\+bottomdraglaw} is True then the term $|U_{bbl}|$ is set equal to the value in \hyperlink{structmom__set__visc_1_1set__visc__cs_a8e5f72291c7e5c5e602c9b4765d63ae0}{set\+\_\+visc\+\_\+cs.\+drag\+\_\+bg\+\_\+vel} so that $C_d |U_{bbl}|$ becomes a Rayleigh bottom drag. Otherwise $|U_{bbl}|$ is found by averaging the flow over the bottom \hyperlink{structmom__set__visc_1_1set__visc__cs_a20768c8872b3eab25e59dfe7c5ff62a0}{set\+\_\+visc\+\_\+cs.\+hbbl} of the model, adding the amplitude of tides \hyperlink{structmom__set__visc_1_1set__visc__cs_a765c0eec288ee86e83afc501b757a821}{set\+\_\+visc\+\_\+cs.\+tideamp} and a constant \hyperlink{structmom__set__visc_1_1set__visc__cs_a8e5f72291c7e5c5e602c9b4765d63ae0}{set\+\_\+visc\+\_\+cs.\+drag\+\_\+bg\+\_\+vel}. For these calculations the vertical grid at the velocity component locations is found by \[ \begin{array}{ll} \frac{2 h^- h^+}{h^- + h^+} & u \left( h^+ - h^-\right) >= 0 \\ \frac{1}{2} \left( h^- + h^+ \right) & u \left( h^+ - h^-\right) < 0 \end{array} \] which biases towards the thin cell if the thin cell is upwind. Biasing the grid toward thin upwind cells helps increase the effect of viscosity and inhibits flow out of these thin cells.

After diagnosing $|U_{bbl}|$ over a fixed depth an active viscous boundary layer thickness is found using the ideas of Killworth and Edwards, 1999 (hereafter K\+W99). K\+W99 solve the equation \[ \left( \frac{h_{bbl}}{h_f} \right)^2 + \frac{h_{bbl}}{h_N} = 1 \] for the boundary layer depth $h_{bbl}$. Here \[ h_f = \frac{C_n u_*}{f} \] is the rotation controlled boundary layer depth in the absence of stratification. $u_*$ is the surface friction speed given by \[ u_*^2 = C_d |U_{bbl}|^2 \] and is a function of near bottom model flow. \[ h_N = \frac{C_i u_*}{N} = \frac{ (C_i u_* )^2 }{g^\prime} \] is the stratification controlled boundary layer depth. The non-\/dimensional parameters $C_n=0.5$ and $C_i=20$ are suggested by Zilitinkevich and Mironov, 1996.

If a Richardson number dependent mixing scheme is being used, as indicated by \hyperlink{structmom__set__visc_1_1set__visc__cs_abf6d0c52b4c7697848e11790fed715dc}{set\+\_\+visc\+\_\+cs.\+rino\+\_\+mix}, then the boundary layer thickness is bounded to be no larger than a half of \hyperlink{structmom__set__visc_1_1set__visc__cs_a20768c8872b3eab25e59dfe7c5ff62a0}{set\+\_\+visc\+\_\+cs.\+hbbl} .

\begin{DoxyRefDesc}{Todo}
\item[\hyperlink{todo__todo000007}{Todo}]Channel drag needs to be explained\end{DoxyRefDesc}


A B\+BL viscosity is calculated so that the no-\/slip boundary condition in the vertical viscosity solver implies the stress $\mathbf{\tau}_b$.\hypertarget{namespacemom__set__visc_set_viscous_BBL_ref}{}\subsubsection{References}\label{namespacemom__set__visc_set_viscous_BBL_ref}
\begin{DoxyItemize}
\item Killworth, P. D., and N. R. Edwards, 1999\+: A Turbulent Bottom Boundary Layer Code for Use in Numerical Ocean Models. J. Phys. Oceanogr., 29, 1221-\/1238, \href{https://doi.org/10.1175/1520-0485(1999)029<1221:ATBBLC>2.0.CO;2}{\tt doi\+:10.\+1175/1520-\/0485(1999)029$<$1221\+:A\+T\+B\+B\+LC$>$2.\+0.\+CO;2} \item Zilitinkevich, S., Mironov, D.\+V., 1996\+: A multi-\/limit formulation for the equilibrium depth of a stably stratified boundary layer. Boundary-\/\+Layer Meteorology 81, 325-\/351. \href{https://doi.org/10.1007/BF02430334}{\tt doi\+:10.\+1007/\+B\+F02430334}\end{DoxyItemize}

\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in,out}  & {\em g} & The ocean\textquotesingle{}s grid structure.\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure.\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em u} & The zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em v} & The meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em tv} & A structure containing pointers to any available thermodynamic fields. Absent fields have N\+U\+LL ptrs..\\
\hline
\mbox{\tt in,out}  & {\em visc} & A structure containing vertical viscosities and related fields.\\
\hline
 & {\em cs} & The control structure returned by a previous call to set\+\_\+visc\+\_\+init.\\
\hline
\mbox{\tt in}  & {\em symmetrize} & If present and true, do extra calculations of those values in visc that would be calculated with symmetric memory. \\
\hline
\end{DoxyParams}


Definition at line 193 of file M\+O\+M\+\_\+set\+\_\+viscosity.\+F90.


\begin{DoxyCode}
193   \textcolor{keywordtype}{type}(ocean\_grid\_type),    \textcolor{keywordtype}{intent(inout)} :: g\textcolor{comment}{    !< The ocean's grid structure.}
194   \textcolor{keywordtype}{type}(verticalgrid\_type),  \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure.}
195   \textcolor{keywordtype}{type}(unit\_scale\_type),    \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{   !< A dimensional unit scaling type}
196   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))}, &
197                             \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{    !< The zonal velocity [L T-1 ~> m s-1].}
198   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))}, &
199                             \textcolor{keywordtype}{intent(in)}    :: v\textcolor{comment}{    !< The meridional velocity [L T-1 ~> m s-1].}
200   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))},  &
201                             \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2].}
202   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),    \textcolor{keywordtype}{intent(in)}    :: tv\textcolor{comment}{   !< A structure containing pointers to any}
203 \textcolor{comment}{                                                  !! available thermodynamic fields. Absent fields}
204 \textcolor{comment}{                                                  !! have NULL ptrs..}
205   \textcolor{keywordtype}{type}(vertvisc\_type),      \textcolor{keywordtype}{intent(inout)} :: visc\textcolor{comment}{ !< A structure containing vertical viscosities and}
206 \textcolor{comment}{                                                  !! related fields.}
207   \textcolor{keywordtype}{type}(set\_visc\_cs),        \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{   !< The control structure returned by a previous}
208 \textcolor{comment}{                                                  !! call to set\_visc\_init.}
209   \textcolor{keywordtype}{logical},        \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: symmetrize\textcolor{comment}{ !< If present and true, do extra calculations}
210 \textcolor{comment}{                                                  !! of those values in visc that would be}
211 \textcolor{comment}{                                                  !! calculated with symmetric memory.}
212 
213   \textcolor{comment}{! Local variables}
214   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
215     ustar, &    \textcolor{comment}{!   The bottom friction velocity [Z T-1 ~> m s-1].}
216     t\_eos, &    \textcolor{comment}{!   The temperature used to calculate the partial derivatives}
217                 \textcolor{comment}{! of density with T and S [degC].}
218     s\_eos, &    \textcolor{comment}{!   The salinity used to calculate the partial derivatives}
219                 \textcolor{comment}{! of density with T and S [ppt].}
220     dr\_dt, &    \textcolor{comment}{!   Partial derivative of the density in the bottom boundary}
221                 \textcolor{comment}{! layer with temperature [R degC-1 ~> kg m-3 degC-1].}
222     dr\_ds, &    \textcolor{comment}{!   Partial derivative of the density in the bottom boundary}
223                 \textcolor{comment}{! layer with salinity [R ppt-1 ~> kg m-3 ppt-1].}
224     press       \textcolor{comment}{!   The pressure at which dR\_dT and dR\_dS are evaluated [R L2 T-2 ~> Pa].}
225   \textcolor{keywordtype}{real} :: htot      \textcolor{comment}{! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].}
226   \textcolor{keywordtype}{real} :: htot\_vel  \textcolor{comment}{! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].}
227 
228   \textcolor{keywordtype}{real} :: rhtot \textcolor{comment}{! Running sum of thicknesses times the layer potential}
229                 \textcolor{comment}{! densities [H R ~> kg m-2 or kg2 m-5].}
230   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G))} :: &
231     d\_u, &      \textcolor{comment}{! Bottom depth interpolated to u points [Z ~> m].}
232     mask\_u      \textcolor{comment}{! A mask that disables any contributions from u points that}
233                 \textcolor{comment}{! are land or past open boundary conditions [nondim], 0 or 1.}
234   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G))} :: &
235     d\_v, &      \textcolor{comment}{! Bottom depth interpolated to v points [Z ~> m].}
236     mask\_v      \textcolor{comment}{! A mask that disables any contributions from v points that}
237                 \textcolor{comment}{! are land or past open boundary conditions [nondim], 0 or 1.}
238   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(G))} :: &
239     h\_at\_vel, & \textcolor{comment}{! Layer thickness at a velocity point, using an upwind-biased}
240                 \textcolor{comment}{! second order accurate estimate based on the previous velocity}
241                 \textcolor{comment}{! direction [H ~> m or kg m-2].}
242     h\_vel, &    \textcolor{comment}{! Arithmetic mean of the layer thicknesses adjacent to a}
243                 \textcolor{comment}{! velocity point [H ~> m or kg m-2].}
244     t\_vel, &    \textcolor{comment}{! Arithmetic mean of the layer temperatures adjacent to a}
245                 \textcolor{comment}{! velocity point [degC].}
246     s\_vel, &    \textcolor{comment}{! Arithmetic mean of the layer salinities adjacent to a}
247                 \textcolor{comment}{! velocity point [ppt].}
248     rml\_vel     \textcolor{comment}{! Arithmetic mean of the layer coordinate densities adjacent}
249                 \textcolor{comment}{! to a velocity point [R ~> kg m-3].}
250 
251   \textcolor{keywordtype}{real} :: h\_vel\_pos        \textcolor{comment}{! The arithmetic mean thickness at a velocity point}
252                            \textcolor{comment}{! plus H\_neglect to avoid 0 values [H ~> m or kg m-2].}
253   \textcolor{keywordtype}{real} :: ustarsq          \textcolor{comment}{! 400 times the square of ustar, times}
254                            \textcolor{comment}{! Rho0 divided by G\_Earth and the conversion}
255                            \textcolor{comment}{! from m to thickness units [H R ~> kg m-2 or kg2 m-5].}
256   \textcolor{keywordtype}{real} :: cdrag\_sqrt\_z     \textcolor{comment}{! Square root of the drag coefficient, times a unit conversion}
257                            \textcolor{comment}{! factor from lateral lengths to vertical depths [Z L-1 ~> 1].}
258   \textcolor{keywordtype}{real} :: cdrag\_sqrt       \textcolor{comment}{! Square root of the drag coefficient [nondim].}
259   \textcolor{keywordtype}{real} :: oldfn            \textcolor{comment}{! The integrated energy required to}
260                            \textcolor{comment}{! entrain up to the bottom of the layer,}
261                            \textcolor{comment}{! divided by G\_Earth [H R ~> kg m-2 or kg2 m-5].}
262   \textcolor{keywordtype}{real} :: dfn              \textcolor{comment}{! The increment in oldfn for entraining}
263                            \textcolor{comment}{! the layer [H R ~> kg m-2 or kg2 m-5].}
264   \textcolor{keywordtype}{real} :: dh               \textcolor{comment}{! The increment in layer thickness from}
265                            \textcolor{comment}{! the present layer [H ~> m or kg m-2].}
266   \textcolor{keywordtype}{real} :: bbl\_thick        \textcolor{comment}{! The thickness of the bottom boundary layer [H ~> m or kg m-2].}
267   \textcolor{keywordtype}{real} :: bbl\_thick\_z      \textcolor{comment}{! The thickness of the bottom boundary layer [Z ~> m].}
268   \textcolor{keywordtype}{real} :: kv\_bbl           \textcolor{comment}{! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].}
269   \textcolor{keywordtype}{real} :: c2f              \textcolor{comment}{! C2f = 2*f at velocity points [T-1 ~> s-1].}
270 
271   \textcolor{keywordtype}{real} :: u\_bg\_sq          \textcolor{comment}{! The square of an assumed background}
272                            \textcolor{comment}{! velocity, for calculating the mean}
273                            \textcolor{comment}{! magnitude near the bottom for use in the}
274                            \textcolor{comment}{! quadratic bottom drag [L2 T-2 ~> m2 s-2].}
275   \textcolor{keywordtype}{real} :: hwtot            \textcolor{comment}{! Sum of the thicknesses used to calculate}
276                            \textcolor{comment}{! the near-bottom velocity magnitude [H ~> m or kg m-2].}
277   \textcolor{keywordtype}{real} :: hutot            \textcolor{comment}{! Running sum of thicknesses times the}
278                            \textcolor{comment}{! velocity magnitudes [H T T-1 ~> m2 s-1 or kg m-1 s-1].}
279   \textcolor{keywordtype}{real} :: thtot            \textcolor{comment}{! Running sum of thickness times temperature [degC H ~> degC m or degC kg m-2].}
280   \textcolor{keywordtype}{real} :: shtot            \textcolor{comment}{! Running sum of thickness times salinity [ppt H ~> ppt m or ppt kg m-2].}
281   \textcolor{keywordtype}{real} :: hweight          \textcolor{comment}{! The thickness of a layer that is within Hbbl}
282                            \textcolor{comment}{! of the bottom [H ~> m or kg m-2].}
283   \textcolor{keywordtype}{real} :: v\_at\_u, u\_at\_v   \textcolor{comment}{! v at a u point or vice versa [L T-1 ~> m s-1].}
284   \textcolor{keywordtype}{real} :: rho0x400\_g       \textcolor{comment}{! 400*Rho0/G\_Earth, times unit conversion factors}
285                            \textcolor{comment}{! [R T2 H Z-2 ~> kg s2 m-4 or kg2 s2 m-7].}
286                            \textcolor{comment}{! The 400 is a constant proposed by Killworth and Edwards, 1999.}
287   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),max(GV%nk\_rho\_varies,1))} :: &
288     rml                    \textcolor{comment}{! The mixed layer coordinate density [R ~> kg m-3].}
289   \textcolor{keywordtype}{real} :: p\_ref(szi\_(g))   \textcolor{comment}{!   The pressure used to calculate the coordinate}
290                            \textcolor{comment}{! density [R L2 T-2 ~> Pa] (usually set to 2e7 Pa = 2000 dbar).}
291 
292   \textcolor{keywordtype}{real} :: d\_vel            \textcolor{comment}{! The bottom depth at a velocity point [H ~> m or kg m-2].}
293   \textcolor{keywordtype}{real} :: dp, dm           \textcolor{comment}{! The depths at the edges of a velocity cell [H ~> m or kg m-2].}
294   \textcolor{keywordtype}{real} :: a                \textcolor{comment}{! a is the curvature of the bottom depth across a}
295                            \textcolor{comment}{! cell, times the cell width squared [H ~> m or kg m-2].}
296   \textcolor{keywordtype}{real} :: a\_3, a\_12        \textcolor{comment}{! a/3 and a/12 [H ~> m or kg m-2].}
297   \textcolor{keywordtype}{real} :: c24\_a            \textcolor{comment}{! 24/a [H-1 ~> m-1 or m2 kg-1].}
298   \textcolor{keywordtype}{real} :: slope            \textcolor{comment}{! The absolute value of the bottom depth slope across}
299                            \textcolor{comment}{! a cell times the cell width [H ~> m or kg m-2].}
300   \textcolor{keywordtype}{real} :: apb\_4a, ax2\_3apb \textcolor{comment}{! Various nondimensional ratios of a and slope.}
301   \textcolor{keywordtype}{real} :: a2x48\_apb3, iapb, ibma\_2 \textcolor{comment}{! Combinations of a and slope [H-1 ~> m-1 or m2 kg-1].}
302   \textcolor{comment}{! All of the following "volumes" have units of thickness because they are normalized}
303   \textcolor{comment}{! by the full horizontal area of a velocity cell.}
304   \textcolor{keywordtype}{real} :: vol\_open         \textcolor{comment}{! The cell volume above which it is open [H ~> m or kg m-2].}
305   \textcolor{keywordtype}{real} :: vol\_direct       \textcolor{comment}{! With less than Vol\_direct [H ~> m or kg m-2], there is a direct}
306                            \textcolor{comment}{! solution of a cubic equation for L.}
307   \textcolor{keywordtype}{real} :: vol\_2\_reg        \textcolor{comment}{! The cell volume above which there are two separate}
308                            \textcolor{comment}{! open areas that must be integrated [H ~> m or kg m-2].}
309   \textcolor{keywordtype}{real} :: vol              \textcolor{comment}{! The volume below the interface whose normalized}
310                            \textcolor{comment}{! width is being sought [H ~> m or kg m-2].}
311   \textcolor{keywordtype}{real} :: vol\_below        \textcolor{comment}{! The volume below the interface below the one that}
312                            \textcolor{comment}{! is currently under consideration [H ~> m or kg m-2].}
313   \textcolor{keywordtype}{real} :: vol\_err          \textcolor{comment}{! The error in the volume with the latest estimate of}
314                            \textcolor{comment}{! L, or the error for the interface below [H ~> m or kg m-2].}
315   \textcolor{keywordtype}{real} :: vol\_quit         \textcolor{comment}{! The volume error below which to quit iterating [H ~> m or kg m-2].}
316   \textcolor{keywordtype}{real} :: vol\_tol          \textcolor{comment}{! A volume error tolerance [H ~> m or kg m-2].}
317   \textcolor{keywordtype}{real} :: l(szk\_(g)+1)     \textcolor{comment}{! The fraction of the full cell width that is open at}
318                            \textcolor{comment}{! the depth of each interface [nondim].}
319   \textcolor{keywordtype}{real} :: l\_direct         \textcolor{comment}{! The value of L above volume Vol\_direct [nondim].}
320   \textcolor{keywordtype}{real} :: l\_max, l\_min     \textcolor{comment}{! Upper and lower bounds on the correct value for L  [nondim].}
321   \textcolor{keywordtype}{real} :: vol\_err\_max      \textcolor{comment}{! The volume errors for the upper and lower bounds on}
322   \textcolor{keywordtype}{real} :: vol\_err\_min      \textcolor{comment}{! the correct value for L [H ~> m or kg m-2].}
323   \textcolor{keywordtype}{real} :: vol\_0            \textcolor{comment}{! A deeper volume with known width L0 [H ~> m or kg m-2].}
324   \textcolor{keywordtype}{real} :: l0               \textcolor{comment}{! The value of L above volume Vol\_0 [nondim].}
325   \textcolor{keywordtype}{real} :: dvol             \textcolor{comment}{! vol - Vol\_0 [H ~> m or kg m-2].}
326   \textcolor{keywordtype}{real} :: dv\_dl2           \textcolor{comment}{! The partial derivative of volume with L squared}
327                            \textcolor{comment}{! evaluated at L=L0 [H ~> m or kg m-2].}
328   \textcolor{keywordtype}{real} :: h\_neglect        \textcolor{comment}{! A thickness that is so small it is usually lost}
329                            \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
330   \textcolor{keywordtype}{real} :: usth             \textcolor{comment}{! ustar converted to units of H T-1 [H T-1 ~> m s-1 or kg m-2 s-1].}
331   \textcolor{keywordtype}{real} :: root             \textcolor{comment}{! A temporary variable [H T-1 ~> m s-1 or kg m-2 s-1].}
332 
333   \textcolor{keywordtype}{real} :: cell\_width       \textcolor{comment}{! The transverse width of the velocity cell [L ~> m].}
334   \textcolor{keywordtype}{real} :: rayleigh         \textcolor{comment}{! A nondimensional value that is multiplied by the layer's}
335                            \textcolor{comment}{! velocity magnitude to give the Rayleigh drag velocity, times}
336                            \textcolor{comment}{! a lateral to vertical distance conversion factor [Z L-1 ~> 1].}
337   \textcolor{keywordtype}{real} :: gam              \textcolor{comment}{! The ratio of the change in the open interface width}
338                            \textcolor{comment}{! to the open interface width atop a cell [nondim].}
339   \textcolor{keywordtype}{real} :: bbl\_frac         \textcolor{comment}{! The fraction of a layer's drag that goes into the}
340                            \textcolor{comment}{! viscous bottom boundary layer [nondim].}
341   \textcolor{keywordtype}{real} :: bbl\_visc\_frac    \textcolor{comment}{! The fraction of all the drag that is expressed as}
342                            \textcolor{comment}{! a viscous bottom boundary layer [nondim].}
343   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: c1\_3 = 1.0/3.0, c1\_6 = 1.0/6.0, c1\_12 = 1.0/12.0
344   \textcolor{keywordtype}{real} :: c2pi\_3           \textcolor{comment}{! An irrational constant, 2/3 pi.}
345   \textcolor{keywordtype}{real} :: tmp              \textcolor{comment}{! A temporary variable.}
346   \textcolor{keywordtype}{real} :: tmp\_val\_m1\_to\_p1
347   \textcolor{keywordtype}{real} :: curv\_tol         \textcolor{comment}{! Numerator of curvature cubed, used to estimate}
348                            \textcolor{comment}{! accuracy of a single L(:) Newton iteration}
349   \textcolor{keywordtype}{logical} :: use\_l0, do\_one\_l\_iter    \textcolor{comment}{! Control flags for L(:) Newton iteration}
350   \textcolor{keywordtype}{logical} :: use\_bbl\_eos, do\_i(szib\_(g))
351   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(2)} :: eosdom \textcolor{comment}{! The computational domain for the equation of state}
352   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, m, n, k2, nkmb, nkml
353   \textcolor{keywordtype}{integer} :: itt, maxitt=20
354   \textcolor{keywordtype}{type}(ocean\_obc\_type), \textcolor{keywordtype}{pointer} :: obc => null()
355 
356   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
357   isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
358   nkmb = gv%nk\_rho\_varies ; nkml = gv%nkml
359   h\_neglect = gv%H\_subroundoff
360   rho0x400\_g = 400.0*(gv%Rho0 / (us%L\_to\_Z**2 * gv%g\_Earth)) * gv%Z\_to\_H
361   vol\_quit = 0.9*gv%Angstrom\_H + h\_neglect
362   c2pi\_3 = 8.0*atan(1.0)/3.0
363 
364   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal,\textcolor{stringliteral}{"MOM\_set\_viscosity(BBL): "}//&
365          \textcolor{stringliteral}{"Module must be initialized before it is used."})
366   \textcolor{keywordflow}{if} (.not.cs%bottomdraglaw) \textcolor{keywordflow}{return}
367 
368   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(symmetrize)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (symmetrize) \textcolor{keywordflow}{then}
369     jsq = js-1 ; isq = is-1
370 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
371 
372   \textcolor{keywordflow}{if} (cs%debug) \textcolor{keywordflow}{then}
373     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"Start set\_viscous\_BBL [uv]"}, u, v, g%HI, haloshift=1, scale=us%L\_T\_to\_m\_s)
374     \textcolor{keyword}{call }hchksum(h,\textcolor{stringliteral}{"Start set\_viscous\_BBL h"}, g%HI, haloshift=1, scale=gv%H\_to\_m)
375     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%T)) \textcolor{keyword}{call }hchksum(tv%T, \textcolor{stringliteral}{"Start set\_viscous\_BBL T"}, g%HI, haloshift=1)
376     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%S)) \textcolor{keyword}{call }hchksum(tv%S, \textcolor{stringliteral}{"Start set\_viscous\_BBL S"}, g%HI, haloshift=1)
377 \textcolor{keywordflow}{  endif}
378 
379   use\_bbl\_eos = \textcolor{keyword}{associated}(tv%eqn\_of\_state) .and. cs%BBL\_use\_EOS
380   obc => cs%OBC
381 
382   u\_bg\_sq = cs%drag\_bg\_vel * cs%drag\_bg\_vel
383   cdrag\_sqrt = sqrt(cs%cdrag)
384   cdrag\_sqrt\_z = us%L\_to\_Z * sqrt(cs%cdrag)
385   k2 = max(nkmb+1, 2)
386 
387 \textcolor{comment}{!  With a linear drag law, the friction velocity is already known.}
388 \textcolor{comment}{!  if (CS%linear\_drag) ustar(:) = cdrag\_sqrt\_Z*CS%drag\_bg\_vel}
389 
390   \textcolor{keywordflow}{if} ((nkml>0) .and. .not.use\_bbl\_eos) \textcolor{keywordflow}{then}
391     eosdom(1) = isq - (g%isd-1) ;  eosdom(2) = g%iec+1 - (g%isd-1)
392     \textcolor{keywordflow}{do} i=isq,ieq+1 ; p\_ref(i) = tv%P\_Ref ;\textcolor{keywordflow}{ enddo}
393     \textcolor{comment}{!$OMP parallel do default(shared)}
394     \textcolor{keywordflow}{do} k=1,nkmb ; \textcolor{keywordflow}{do} j=jsq,jeq+1
395       \textcolor{keyword}{call }calculate\_density(tv%T(:,j,k), tv%S(:,j,k), p\_ref, rml(:,j,k), tv%eqn\_of\_state, &
396                              eosdom)
397 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo}
398 \textcolor{keywordflow}{  endif}
399 
400   \textcolor{comment}{!$OMP parallel do default(shared)}
401   \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is-1,ie+1
402     d\_v(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))
403     mask\_v(i,j) = g%mask2dCv(i,j)
404 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
405   \textcolor{comment}{!$OMP parallel do default(shared)}
406   \textcolor{keywordflow}{do} j=js-1,je+1 ; \textcolor{keywordflow}{do} i=is-1,ie
407     d\_u(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))
408     mask\_u(i,j) = g%mask2dCu(i,j)
409 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
410 
411   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} n=1,obc%number\_of\_segments
412     \textcolor{keywordflow}{if} (.not. obc%segment(n)%on\_pe) cycle
413     \textcolor{comment}{! Use a one-sided projection of bottom depths at OBC points.}
414     i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
415     \textcolor{keywordflow}{if} (obc%segment(n)%is\_N\_or\_S .and. (j >= js-1) .and. (j <= je)) \textcolor{keywordflow}{then}
416       \textcolor{keywordflow}{do} i = max(is-1,obc%segment(n)%HI%isd), min(ie+1,obc%segment(n)%HI%ied)
417         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_n) d\_v(i,j) = g%bathyT(i,j)
418         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_s) d\_v(i,j) = g%bathyT(i,j+1)
419 \textcolor{keywordflow}{      enddo}
420     \textcolor{keywordflow}{elseif} (obc%segment(n)%is\_E\_or\_W .and. (i >= is-1) .and. (i <= ie)) \textcolor{keywordflow}{then}
421       \textcolor{keywordflow}{do} j = max(js-1,obc%segment(n)%HI%jsd), min(je+1,obc%segment(n)%HI%jed)
422         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_e) d\_u(i,j) = g%bathyT(i,j)
423         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_w) d\_u(i,j) = g%bathyT(i+1,j)
424 \textcolor{keywordflow}{      enddo}
425 \textcolor{keywordflow}{    endif}
426 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ endif}
427   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} n=1,obc%number\_of\_segments
428     \textcolor{comment}{! Now project bottom depths across cell-corner points in the OBCs.  The two}
429     \textcolor{comment}{! projections have to occur in sequence and can not be combined easily.}
430     \textcolor{keywordflow}{if} (.not. obc%segment(n)%on\_pe) cycle
431     \textcolor{comment}{! Use a one-sided projection of bottom depths at OBC points.}
432     i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
433     \textcolor{keywordflow}{if} (obc%segment(n)%is\_N\_or\_S .and. (j >= js-1) .and. (j <= je)) \textcolor{keywordflow}{then}
434       \textcolor{keywordflow}{do} i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
435         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_n) \textcolor{keywordflow}{then}
436           d\_u(i,j+1) = d\_u(i,j) ;  mask\_u(i,j+1) = 0.0
437         \textcolor{keywordflow}{elseif} (obc%segment(n)%direction == obc\_direction\_s) \textcolor{keywordflow}{then}
438           d\_u(i,j) = d\_u(i,j+1) ; mask\_u(i,j) = 0.0
439 \textcolor{keywordflow}{        endif}
440 \textcolor{keywordflow}{      enddo}
441     \textcolor{keywordflow}{elseif} (obc%segment(n)%is\_E\_or\_W .and. (i >= is-1) .and. (i <= ie)) \textcolor{keywordflow}{then}
442       \textcolor{keywordflow}{do} j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
443         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_e) \textcolor{keywordflow}{then}
444           d\_v(i+1,j) = d\_v(i,j) ; mask\_v(i+1,j) = 0.0
445         \textcolor{keywordflow}{elseif} (obc%segment(n)%direction == obc\_direction\_w) \textcolor{keywordflow}{then}
446           d\_v(i,j) = d\_v(i+1,j) ;  mask\_v(i,j) = 0.0
447 \textcolor{keywordflow}{        endif}
448 \textcolor{keywordflow}{      enddo}
449 \textcolor{keywordflow}{    endif}
450 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ endif}
451 
452   \textcolor{keywordflow}{if} (.not.use\_bbl\_eos) rml\_vel(:,:) = 0.0
453 
454   \textcolor{comment}{!$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,nz,nkmb, &}
455   \textcolor{comment}{!$OMP                                     nkml,Isq,Ieq,Jsq,Jeq,h\_neglect,Rho0x400\_G,C2pi\_3, &}
456   \textcolor{comment}{!$OMP                                     U\_bg\_sq,cdrag\_sqrt\_Z,cdrag\_sqrt,K2,use\_BBL\_EOS,   &}
457   \textcolor{comment}{!$OMP                                     OBC,maxitt,D\_u,D\_v,mask\_u,mask\_v) &}
458   \textcolor{comment}{!$OMP                              firstprivate(Vol\_quit)}
459   \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} m=1,2
460 
461     \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then}
462       \textcolor{comment}{! m=1 refers to u-points}
463       \textcolor{keywordflow}{if} (j<g%Jsc) cycle
464       is = isq ; ie = ieq
465       \textcolor{keywordflow}{do} i=is,ie
466         do\_i(i) = .false.
467         \textcolor{keywordflow}{if} (g%mask2dCu(i,j) > 0) do\_i(i) = .true.
468 \textcolor{keywordflow}{      enddo}
469     \textcolor{keywordflow}{else}
470       \textcolor{comment}{! m=2 refers to v-points}
471       is = g%isc ; ie = g%iec
472       \textcolor{keywordflow}{do} i=is,ie
473         do\_i(i) = .false.
474         \textcolor{keywordflow}{if} (g%mask2dCv(i,j) > 0) do\_i(i) = .true.
475 \textcolor{keywordflow}{      enddo}
476 \textcolor{keywordflow}{    endif}
477 
478     \textcolor{comment}{! Calculate thickness at velocity points (u or v depending on value of m).}
479     \textcolor{comment}{! Also interpolate the ML density or T/S properties.}
480     \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then} \textcolor{comment}{! u-points}
481       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
482         \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
483           \textcolor{keywordflow}{if} (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) \textcolor{keywordflow}{then}
484             \textcolor{comment}{! If the flow is from thin to thick then bias towards the thinner thickness}
485             h\_at\_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
486                             (h(i,j,k) + h(i+1,j,k) + h\_neglect)
487           \textcolor{keywordflow}{else}
488             \textcolor{comment}{! If the flow is from thick to thin then use the simple average thickness}
489             h\_at\_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
490 \textcolor{keywordflow}{          endif}
491 \textcolor{keywordflow}{        endif}
492         h\_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
493 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
494       \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
495         \textcolor{comment}{! Perhaps these should be thickness weighted.}
496         t\_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
497         s\_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
498 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ; \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{do} k=1,nkmb ; \textcolor{keywordflow}{do} i=is,ie
499         rml\_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i+1,j,k))
500 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
501     \textcolor{keywordflow}{else} \textcolor{comment}{! v-points}
502       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
503         \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
504           \textcolor{keywordflow}{if} (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) \textcolor{keywordflow}{then}
505             \textcolor{comment}{! If the flow is from thin to thick then bias towards the thinner thickness}
506             h\_at\_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
507                             (h(i,j,k) + h(i,j+1,k) + h\_neglect)
508           \textcolor{keywordflow}{else}
509             \textcolor{comment}{! If the flow is from thick to thin then use the simple average thickness}
510             h\_at\_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
511 \textcolor{keywordflow}{          endif}
512 \textcolor{keywordflow}{        endif}
513         h\_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
514 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
515       \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
516         t\_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
517         s\_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
518 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ; \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{do} k=1,nkmb ; \textcolor{keywordflow}{do} i=is,ie
519         rml\_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i,j+1,k))
520 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
521 \textcolor{keywordflow}{    endif}
522 
523     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (obc%number\_of\_segments > 0) \textcolor{keywordflow}{then}
524       \textcolor{comment}{! Apply a zero gradient projection of thickness across OBC points.}
525       \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then}
526         \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_u(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
527           \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_e) \textcolor{keywordflow}{then}
528             \textcolor{keywordflow}{do} k=1,nz
529               h\_at\_vel(i,k) = h(i,j,k) ; h\_vel(i,k) = h(i,j,k)
530 \textcolor{keywordflow}{            enddo}
531             \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then}
532               \textcolor{keywordflow}{do} k=1,nz
533                 t\_vel(i,k) = tv%T(i,j,k) ; s\_vel(i,k) = tv%S(i,j,k)
534 \textcolor{keywordflow}{              enddo}
535             \textcolor{keywordflow}{else}
536               \textcolor{keywordflow}{do} k=1,nkmb
537                 rml\_vel(i,k) = rml(i,j,k)
538 \textcolor{keywordflow}{              enddo}
539 \textcolor{keywordflow}{            endif}
540           \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_w) \textcolor{keywordflow}{then}
541             \textcolor{keywordflow}{do} k=1,nz
542               h\_at\_vel(i,k) = h(i+1,j,k) ; h\_vel(i,k) = h(i+1,j,k)
543 \textcolor{keywordflow}{            enddo}
544             \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then}
545               \textcolor{keywordflow}{do} k=1,nz
546                 t\_vel(i,k) = tv%T(i+1,j,k) ; s\_vel(i,k) = tv%S(i+1,j,k)
547 \textcolor{keywordflow}{              enddo}
548             \textcolor{keywordflow}{else}
549               \textcolor{keywordflow}{do} k=1,nkmb
550                 rml\_vel(i,k) = rml(i+1,j,k)
551 \textcolor{keywordflow}{              enddo}
552 \textcolor{keywordflow}{            endif}
553 \textcolor{keywordflow}{          endif}
554 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo}
555       \textcolor{keywordflow}{else}
556         \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_v(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
557           \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_n) \textcolor{keywordflow}{then}
558             \textcolor{keywordflow}{do} k=1,nz
559               h\_at\_vel(i,k) = h(i,j,k) ; h\_vel(i,k) = h(i,j,k)
560 \textcolor{keywordflow}{            enddo}
561             \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then}
562               \textcolor{keywordflow}{do} k=1,nz
563                 t\_vel(i,k) = tv%T(i,j,k) ; s\_vel(i,k) = tv%S(i,j,k)
564 \textcolor{keywordflow}{              enddo}
565             \textcolor{keywordflow}{else}
566               \textcolor{keywordflow}{do} k=1,nkmb
567                 rml\_vel(i,k) = rml(i,j,k)
568 \textcolor{keywordflow}{              enddo}
569 \textcolor{keywordflow}{            endif}
570           \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_s) \textcolor{keywordflow}{then}
571             \textcolor{keywordflow}{do} k=1,nz
572               h\_at\_vel(i,k) = h(i,j+1,k) ; h\_vel(i,k) = h(i,j+1,k)
573 \textcolor{keywordflow}{            enddo}
574             \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then}
575               \textcolor{keywordflow}{do} k=1,nz
576                 t\_vel(i,k) = tv%T(i,j+1,k) ; s\_vel(i,k) = tv%S(i,j+1,k)
577 \textcolor{keywordflow}{              enddo}
578             \textcolor{keywordflow}{else}
579               \textcolor{keywordflow}{do} k=1,nkmb
580                 rml\_vel(i,k) = rml(i,j+1,k)
581 \textcolor{keywordflow}{              enddo}
582 \textcolor{keywordflow}{            endif}
583 \textcolor{keywordflow}{          endif}
584 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo}
585 \textcolor{keywordflow}{      endif}
586 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
587 
588     \textcolor{keywordflow}{if} (use\_bbl\_eos .or. .not.cs%linear\_drag) \textcolor{keywordflow}{then}
589       \textcolor{comment}{! Calculate the mean velocity magnitude over the bottommost CS%Hbbl of}
590       \textcolor{comment}{! the water column for determining the quadratic bottom drag.}
591       \textcolor{comment}{! Used in ustar(i)}
592       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
593         htot\_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
594         thtot = 0.0 ; shtot = 0.0
595         \textcolor{keywordflow}{do} k=nz,1,-1
596 
597           \textcolor{keywordflow}{if} (htot\_vel>=cs%Hbbl) \textcolor{keywordflow}{exit} \textcolor{comment}{! terminate the k loop}
598 
599           hweight = min(cs%Hbbl - htot\_vel, h\_at\_vel(i,k))
600           \textcolor{keywordflow}{if} (hweight < 1.5*gv%Angstrom\_H + h\_neglect) cycle
601 
602           htot\_vel  = htot\_vel + h\_at\_vel(i,k)
603           hwtot = hwtot + hweight
604 
605           \textcolor{keywordflow}{if} ((.not.cs%linear\_drag) .and. (hweight >= 0.0)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then}
606             v\_at\_u = set\_v\_at\_u(v, h, g, i, j, k, mask\_v, obc)
607             \textcolor{keywordflow}{if} (cs%BBL\_use\_tidal\_bg) \textcolor{keywordflow}{then}
608               u\_bg\_sq = 0.5*( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
609                               g%mask2dT(i+1,j)*(cs%tideamp(i+1,j)*cs%tideamp(i+1,j)) )
610 \textcolor{keywordflow}{            endif}
611             hutot = hutot + hweight * sqrt(u(i,j,k)*u(i,j,k) + &
612                                            v\_at\_u*v\_at\_u + u\_bg\_sq)
613           \textcolor{keywordflow}{else}
614             u\_at\_v = set\_u\_at\_v(u, h, g, i, j, k, mask\_u, obc)
615             \textcolor{keywordflow}{if} (cs%BBL\_use\_tidal\_bg) \textcolor{keywordflow}{then}
616               u\_bg\_sq = 0.5*( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
617                               g%mask2dT(i,j+1)*(cs%tideamp(i,j+1)*cs%tideamp(i,j+1)) )
618 \textcolor{keywordflow}{            endif}
619             hutot = hutot + hweight * sqrt(v(i,j,k)*v(i,j,k) + &
620                                            u\_at\_v*u\_at\_v + u\_bg\_sq)
621 \textcolor{keywordflow}{          endif} ;\textcolor{keywordflow}{ endif}
622 
623           \textcolor{keywordflow}{if} (use\_bbl\_eos .and. (hweight >= 0.0)) \textcolor{keywordflow}{then}
624             thtot = thtot + hweight * t\_vel(i,k)
625             shtot = shtot + hweight * s\_vel(i,k)
626 \textcolor{keywordflow}{          endif}
627 \textcolor{keywordflow}{        enddo} \textcolor{comment}{! end of k loop}
628 
629         \textcolor{comment}{! Set u* based on u*^2 = Cdrag u\_bbl^2}
630         \textcolor{keywordflow}{if} (.not.cs%linear\_drag .and. (hwtot > 0.0)) \textcolor{keywordflow}{then}
631           ustar(i) = cdrag\_sqrt\_z*hutot/hwtot
632         \textcolor{keywordflow}{else}
633           ustar(i) = cdrag\_sqrt\_z*cs%drag\_bg\_vel
634 \textcolor{keywordflow}{        endif}
635 
636         \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (hwtot > 0.0) \textcolor{keywordflow}{then}
637           t\_eos(i) = thtot/hwtot ; s\_eos(i) = shtot/hwtot
638         \textcolor{keywordflow}{else}
639           t\_eos(i) = 0.0 ; s\_eos(i) = 0.0
640 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ endif}
641 
642         \textcolor{comment}{! Diagnostic BBL flow speed at u- and v-points.}
643         \textcolor{keywordflow}{if} (cs%id\_bbl\_u>0 .and. m==1) \textcolor{keywordflow}{then}
644           \textcolor{keywordflow}{if} (hwtot > 0.0) cs%bbl\_u(i,j) = hutot/hwtot
645         \textcolor{keywordflow}{elseif} (cs%id\_bbl\_v>0 .and. m==2) \textcolor{keywordflow}{then}
646           \textcolor{keywordflow}{if} (hwtot > 0.0) cs%bbl\_v(i,j) = hutot/hwtot
647 \textcolor{keywordflow}{        endif}
648 
649 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
650     \textcolor{keywordflow}{else}
651       \textcolor{keywordflow}{do} i=is,ie ; ustar(i) = cdrag\_sqrt\_z*cs%drag\_bg\_vel ;\textcolor{keywordflow}{ enddo}
652 \textcolor{keywordflow}{    endif} \textcolor{comment}{! Not linear\_drag}
653 
654     \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then}
655       \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%p\_surf)) \textcolor{keywordflow}{then}
656         \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie ; press(i) = 0.5*(tv%p\_surf(i,j) + tv%p\_surf(i+1,j)) ;\textcolor{keywordflow}{ enddo}
657         \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{do} i=is,ie ; press(i) = 0.5*(tv%p\_surf(i,j) + tv%p\_surf(i,j+1)) ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
658       \textcolor{keywordflow}{else}
659         \textcolor{keywordflow}{do} i=is,ie ; press(i) = 0.0 ;\textcolor{keywordflow}{ enddo}
660 \textcolor{keywordflow}{      endif}
661       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (.not.do\_i(i)) \textcolor{keywordflow}{then} ; t\_eos(i) = 0.0 ; s\_eos(i) = 0.0 ;\textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo}
662       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
663         press(i) = press(i) + (gv%H\_to\_RZ*gv%g\_Earth) * h\_vel(i,k)
664 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
665       \textcolor{keyword}{call }calculate\_density\_derivs(t\_eos, s\_eos, press, dr\_dt, dr\_ds, tv%eqn\_of\_state, &
666                                     (/is-g%IsdB+1,ie-g%IsdB+1/) )
667 \textcolor{keywordflow}{    endif}
668 
669     \textcolor{comment}{! Find a BBL thickness given by equation 2.20 of Killworth and Edwards, 1999:}
670     \textcolor{comment}{!    ( f h / Cn u* )^2 + ( N h / Ci u* ) = 1}
671     \textcolor{comment}{! where Cn=0.5 and Ci=20 (constants suggested by Zilitinkevich and Mironov, 1996).}
672     \textcolor{comment}{! Eq. 2.20 can be expressed in terms of boundary layer thicknesses limited by}
673     \textcolor{comment}{! rotation (h\_f) and stratification (h\_N):}
674     \textcolor{comment}{!    ( h / h\_f )^2 + ( h / h\_N ) = 1}
675     \textcolor{comment}{! When stratification dominates h\_N<<h\_f, and vice versa.}
676     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
677       \textcolor{comment}{! The 400.0 in this expression is the square of a Ci introduced in KW99, eq. 2.22.}
678       ustarsq = rho0x400\_g * ustar(i)**2 \textcolor{comment}{! Note not in units of u*^2 but [H R]}
679       htot = 0.0
680 
681       \textcolor{comment}{! Calculate the thickness of a stratification limited BBL ignoring rotation:}
682       \textcolor{comment}{!   h\_N = Ci u* / N          (limit of KW99 eq. 2.20 for |f|->0)}
683       \textcolor{comment}{! For layer mode, N^2 = g'/h. Since (Ci u*)^2 = (h\_N N)^2 = h\_N g' then}
684       \textcolor{comment}{!   h\_N = (Ci u*)^2 / g'     (KW99, eq, 2.22)}
685       \textcolor{comment}{! Starting from the bottom, integrate the stratification upward until h\_N N balances Ci u*}
686       \textcolor{comment}{! or in layer mode}
687       \textcolor{comment}{!   h\_N Delta rho ~ (Ci u*)^2 rho0 / g}
688       \textcolor{comment}{! where the rhs is stored in variable ustarsq.}
689       \textcolor{comment}{! The method was described in Stephens and Hallberg 2000 (unpublished and lost manuscript).}
690       \textcolor{keywordflow}{if} (use\_bbl\_eos) \textcolor{keywordflow}{then}
691         thtot = 0.0 ; shtot = 0.0 ; oldfn = 0.0
692         \textcolor{keywordflow}{do} k=nz,2,-1
693           \textcolor{keywordflow}{if} (h\_at\_vel(i,k) <= 0.0) cycle
694 
695           \textcolor{comment}{! Delta rho * h\_bbl assuming everything below is homogenized}
696           oldfn = dr\_dt(i)*(thtot - t\_vel(i,k)*htot) + &
697                   dr\_ds(i)*(shtot - s\_vel(i,k)*htot)
698           \textcolor{keywordflow}{if} (oldfn >= ustarsq) \textcolor{keywordflow}{exit}
699 
700           \textcolor{comment}{! Local Delta rho * h\_bbl  at interface}
701           dfn = (dr\_dt(i)*(t\_vel(i,k) - t\_vel(i,k-1)) + &
702                  dr\_ds(i)*(s\_vel(i,k) - s\_vel(i,k-1))) * &
703                 (h\_at\_vel(i,k) + htot)
704 
705           \textcolor{keywordflow}{if} ((oldfn + dfn) <= ustarsq) \textcolor{keywordflow}{then}
706             \textcolor{comment}{! Use whole layer}
707             dh = h\_at\_vel(i,k)
708           \textcolor{keywordflow}{else}
709             \textcolor{comment}{! Use only part of the layer}
710             dh = h\_at\_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
711 \textcolor{keywordflow}{          endif}
712 
713           \textcolor{comment}{! Increment total BBL thickness and cumulative T and S}
714           htot = htot + dh
715           thtot = thtot + t\_vel(i,k)*dh ; shtot = shtot + s\_vel(i,k)*dh
716 \textcolor{keywordflow}{        enddo}
717         \textcolor{keywordflow}{if} ((oldfn < ustarsq) .and. h\_at\_vel(i,1) > 0.0) \textcolor{keywordflow}{then}
718           \textcolor{comment}{! Layer 1 might be part of the BBL.}
719           \textcolor{keywordflow}{if} (dr\_dt(i) * (thtot - t\_vel(i,1)*htot) + &
720               dr\_ds(i) * (shtot - s\_vel(i,1)*htot) < ustarsq) &
721             htot = htot + h\_at\_vel(i,1)
722 \textcolor{keywordflow}{        endif} \textcolor{comment}{! Examination of layer 1.}
723       \textcolor{keywordflow}{else}  \textcolor{comment}{! Use Rlay and/or the coordinate density as density variables.}
724         rhtot = 0.0
725         \textcolor{keywordflow}{do} k=nz,k2,-1
726           oldfn = rhtot - gv%Rlay(k)*htot
727           dfn = (gv%Rlay(k) - gv%Rlay(k-1))*(h\_at\_vel(i,k)+htot)
728 
729           \textcolor{keywordflow}{if} (oldfn >= ustarsq) \textcolor{keywordflow}{then}
730             cycle
731           \textcolor{keywordflow}{elseif} ((oldfn + dfn) <= ustarsq) \textcolor{keywordflow}{then}
732             dh = h\_at\_vel(i,k)
733           \textcolor{keywordflow}{else}
734             dh = h\_at\_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
735 \textcolor{keywordflow}{          endif}
736 
737           htot = htot + dh
738           rhtot = rhtot + gv%Rlay(k)*dh
739 \textcolor{keywordflow}{        enddo}
740         \textcolor{keywordflow}{if} (nkml>0) \textcolor{keywordflow}{then}
741           \textcolor{keywordflow}{do} k=nkmb,2,-1
742             oldfn = rhtot - rml\_vel(i,k)*htot
743             dfn = (rml\_vel(i,k) - rml\_vel(i,k-1)) * (h\_at\_vel(i,k)+htot)
744 
745             \textcolor{keywordflow}{if} (oldfn >= ustarsq) \textcolor{keywordflow}{then}
746               cycle
747             \textcolor{keywordflow}{elseif} ((oldfn + dfn) <= ustarsq) \textcolor{keywordflow}{then}
748               dh = h\_at\_vel(i,k)
749             \textcolor{keywordflow}{else}
750               dh = h\_at\_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
751 \textcolor{keywordflow}{            endif}
752 
753             htot = htot + dh
754             rhtot = rhtot + rml\_vel(i,k)*dh
755 \textcolor{keywordflow}{          enddo}
756           \textcolor{keywordflow}{if} (rhtot - rml\_vel(i,1)*htot < ustarsq) htot = htot + h\_at\_vel(i,1)
757         \textcolor{keywordflow}{else}
758           \textcolor{keywordflow}{if} (rhtot - gv%Rlay(1)*htot < ustarsq) htot = htot + h\_at\_vel(i,1)
759 \textcolor{keywordflow}{        endif}
760 \textcolor{keywordflow}{      endif} \textcolor{comment}{! use\_BBL\_EOS}
761 
762       \textcolor{comment}{! Value of 2*f at u- or v-points.}
763       \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then} ; c2f = g%CoriolisBu(i,j-1) + g%CoriolisBu(i,j)
764       \textcolor{keywordflow}{else} ; c2f = g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j) ;\textcolor{keywordflow}{ endif}
765 
766       \textcolor{comment}{! The thickness of a rotation limited BBL ignoring stratification is}
767       \textcolor{comment}{!   h\_f ~ Cn u* / f        (limit of KW99 eq. 2.20 for N->0).}
768       \textcolor{comment}{! The buoyancy limit of BBL thickness (h\_N) is already in the variable htot from above.}
769       \textcolor{comment}{! Substituting x = h\_N/h into KW99 eq. 2.20 yields the quadratic}
770       \textcolor{comment}{!   x^2 - x = (h\_N / h\_f)^2}
771       \textcolor{comment}{! for which the positive root is}
772       \textcolor{comment}{!   xp = 1/2 + sqrt( 1/4 + (h\_N/h\_f)^2 )}
773       \textcolor{comment}{! and thus h\_bbl = h\_N / xp . Since h\_f = Cn u*/f and Cn=0.5}
774       \textcolor{comment}{!   xp = 1/2 + sqrt( 1/4 + (2 f h\_N/u*)^2 )}
775       \textcolor{comment}{! To avoid dividing by zero if u*=0 then}
776       \textcolor{comment}{!   xp u* = 1/2 u* + sqrt( 1/4 u*^2 + (2 f h\_N)^2 )}
777       \textcolor{keywordflow}{if} (cs%cdrag * u\_bg\_sq <= 0.0) \textcolor{keywordflow}{then}
778         \textcolor{comment}{! This avoids NaNs and overflows, and could be used in all cases,}
779         \textcolor{comment}{! but is not bitwise identical to the current code.}
780         usth = ustar(i)*gv%Z\_to\_H ; root = sqrt(0.25*usth**2 + (htot*c2f)**2)
781         \textcolor{keywordflow}{if} (htot*usth <= (cs%BBL\_thick\_min+h\_neglect) * (0.5*usth + root)) \textcolor{keywordflow}{then}
782           bbl\_thick = cs%BBL\_thick\_min
783         \textcolor{keywordflow}{else}
784           \textcolor{comment}{! The following expression reads}
785           \textcolor{comment}{!   h\_bbl = h\_N u* / ( 1/2 u* + sqrt( 1/4 u*^2 + ( 2 f h\_N )^2 ) )}
786           \textcolor{comment}{! which is h\_bbl = h\_N u*/(xp u*) as described above.}
787           bbl\_thick = (htot * usth) / (0.5*usth + root)
788 \textcolor{keywordflow}{        endif}
789       \textcolor{keywordflow}{else}
790         \textcolor{comment}{! The following expression reads}
791         \textcolor{comment}{!   h\_bbl = h\_N / ( 1/2 + sqrt( 1/4 + ( 2 f h\_N / u* )^2 ) )}
792         \textcolor{comment}{! which is h\_bbl = h\_N/xp as described above.}
793         bbl\_thick = htot / (0.5 + sqrt(0.25 + htot*htot*c2f*c2f/ &
794                                        ((ustar(i)*ustar(i)) * (gv%Z\_to\_H**2)) ) )
795 
796         \textcolor{keywordflow}{if} (bbl\_thick < cs%BBL\_thick\_min) bbl\_thick = cs%BBL\_thick\_min
797 \textcolor{keywordflow}{      endif}
798 
799       \textcolor{comment}{! If there is Richardson number dependent mixing, that determines}
800       \textcolor{comment}{! the vertical extent of the bottom boundary layer, and there is no}
801       \textcolor{comment}{! need to set that scale here.  In fact, viscously reducing the}
802       \textcolor{comment}{! shears over an excessively large region reduces the efficacy of}
803       \textcolor{comment}{! the Richardson number dependent mixing.}
804       \textcolor{comment}{! In other words, if using RiNo\_mix then CS%Hbbl acts as an upper bound on}
805       \textcolor{comment}{! bbl\_thick.}
806       \textcolor{keywordflow}{if} ((bbl\_thick > 0.5*cs%Hbbl) .and. (cs%RiNo\_mix)) bbl\_thick = 0.5*cs%Hbbl
807 
808       \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then}
809         \textcolor{comment}{! The drag within the bottommost bbl\_thick is applied as a part of}
810         \textcolor{comment}{! an enhanced bottom viscosity, while above this the drag is applied}
811         \textcolor{comment}{! directly to the layers in question as a Rayleigh drag term.}
812         \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then}
813           d\_vel = d\_u(i,j)
814           tmp = g%mask2dCu(i,j+1) * d\_u(i,j+1)
815           dp = 2.0 * d\_vel * tmp / (d\_vel + tmp)
816           tmp = g%mask2dCu(i,j-1) * d\_u(i,j-1)
817           dm = 2.0 * d\_vel * tmp / (d\_vel + tmp)
818         \textcolor{keywordflow}{else}
819           d\_vel = d\_v(i,j)
820           tmp = g%mask2dCv(i+1,j) * d\_v(i+1,j)
821           dp = 2.0 * d\_vel * tmp / (d\_vel + tmp)
822           tmp = g%mask2dCv(i-1,j) * d\_v(i-1,j)
823           dm = 2.0 * d\_vel * tmp / (d\_vel + tmp)
824 \textcolor{keywordflow}{        endif}
825         \textcolor{keywordflow}{if} (dm > dp) \textcolor{keywordflow}{then} ; tmp = dp ; dp = dm ; dm = tmp ;\textcolor{keywordflow}{ endif}
826 
827         \textcolor{comment}{! Convert the D's to the units of thickness.}
828         dp = gv%Z\_to\_H*dp ; dm = gv%Z\_to\_H*dm ; d\_vel = gv%Z\_to\_H*d\_vel
829 
830         a\_3 = (dp + dm - 2.0*d\_vel) ; a = 3.0*a\_3 ; a\_12 = 0.25*a\_3
831         slope = dp - dm
832         \textcolor{comment}{! If the curvature is small enough, there is no reason not to assume}
833         \textcolor{comment}{! a uniformly sloping or flat bottom.}
834         \textcolor{keywordflow}{if} (abs(a) < 1e-2*(slope + cs%BBL\_thick\_min)) a = 0.0
835         \textcolor{comment}{! Each cell extends from x=-1/2 to 1/2, and has a topography}
836         \textcolor{comment}{! given by D(x) = a*x^2 + b*x + D - a/12.}
837 
838         \textcolor{comment}{! Calculate the volume above which the entire cell is open and the}
839         \textcolor{comment}{! other volumes at which the equation that is solved for L changes.}
840         \textcolor{keywordflow}{if} (a > 0.0) \textcolor{keywordflow}{then}
841           \textcolor{keywordflow}{if} (slope >= a) \textcolor{keywordflow}{then}
842             vol\_open = d\_vel - dm ; vol\_2\_reg = vol\_open
843           \textcolor{keywordflow}{else}
844             tmp = slope/a
845             vol\_open = 0.25*slope*tmp + c1\_12*a
846             vol\_2\_reg = 0.5*tmp**2 * (a - c1\_3*slope)
847 \textcolor{keywordflow}{          endif}
848           \textcolor{comment}{! Define some combinations of a & b for later use.}
849           c24\_a = 24.0/a ; iapb = 1.0/(a+slope)
850           apb\_4a = (slope+a)/(4.0*a) ; a2x48\_apb3 = (48.0*(a*a))*(iapb**3)
851           ax2\_3apb = 2.0*c1\_3*a*iapb
852         \textcolor{keywordflow}{elseif} (a == 0.0) \textcolor{keywordflow}{then}
853           vol\_open = 0.5*slope
854           \textcolor{keywordflow}{if} (slope > 0) iapb = 1.0/slope
855         \textcolor{keywordflow}{else} \textcolor{comment}{! a < 0.0}
856           vol\_open = d\_vel - dm
857           \textcolor{keywordflow}{if} (slope >= -a) \textcolor{keywordflow}{then}
858             iapb = 1.0e30 ; \textcolor{keywordflow}{if} (slope+a /= 0.0) iapb = 1.0/(a+slope)
859             vol\_direct = 0.0 ; l\_direct = 0.0 ; c24\_a = 0.0
860           \textcolor{keywordflow}{else}
861             c24\_a = 24.0/a ; iapb = 1.0/(a+slope)
862             l\_direct = 1.0 + slope/a \textcolor{comment}{! L\_direct < 1 because a < 0}
863             vol\_direct = -c1\_6*a*l\_direct**3
864 \textcolor{keywordflow}{          endif}
865           ibma\_2 = 2.0 / (slope - a)
866 \textcolor{keywordflow}{        endif}
867 
868         l(nz+1) = 0.0 ; vol = 0.0 ; vol\_err = 0.0 ; bbl\_visc\_frac = 0.0
869         \textcolor{comment}{! Determine the normalized open length at each interface.}
870         \textcolor{keywordflow}{do} k=nz,1,-1
871           vol\_below = vol
872 
873           vol = vol + h\_vel(i,k)
874           h\_vel\_pos = h\_vel(i,k) + h\_neglect
875 
876           \textcolor{keywordflow}{if} (vol >= vol\_open) \textcolor{keywordflow}{then} ; l(k) = 1.0
877           \textcolor{keywordflow}{elseif} (a == 0) \textcolor{keywordflow}{then} \textcolor{comment}{! The bottom has no curvature.}
878             l(k) = sqrt(2.0*vol*iapb)
879           \textcolor{keywordflow}{elseif} (a > 0) \textcolor{keywordflow}{then}
880             \textcolor{comment}{! There may be a minimum depth, and there are}
881             \textcolor{comment}{! analytic expressions for L for all cases.}
882             \textcolor{keywordflow}{if} (vol < vol\_2\_reg) \textcolor{keywordflow}{then}
883               \textcolor{comment}{! In this case, there is a contiguous open region and}
884               \textcolor{comment}{!   vol = 0.5*L^2*(slope + a/3*(3-4L)).}
885               \textcolor{keywordflow}{if} (a2x48\_apb3*vol < 1e-8) \textcolor{keywordflow}{then} \textcolor{comment}{! Could be 1e-7?}
886                 \textcolor{comment}{! There is a very good approximation here for massless layers.}
887                 l0 = sqrt(2.0*vol*iapb) ; l(k) = l0*(1.0 + ax2\_3apb*l0)
888               \textcolor{keywordflow}{else}
889                 l(k) = apb\_4a * (1.0 - &
890                          2.0 * cos(c1\_3*acos(a2x48\_apb3*vol - 1.0) - c2pi\_3))
891 \textcolor{keywordflow}{              endif}
892               \textcolor{comment}{! To check the answers.}
893               \textcolor{comment}{! Vol\_err = 0.5*(L(K)*L(K))*(slope + a\_3*(3.0-4.0*L(K))) - vol}
894             \textcolor{keywordflow}{else} \textcolor{comment}{! There are two separate open regions.}
895               \textcolor{comment}{!   vol = slope^2/4a + a/12 - (a/12)*(1-L)^2*(1+2L)}
896               \textcolor{comment}{! At the deepest volume, L = slope/a, at the top L = 1.}
897               \textcolor{comment}{!L(K) = 0.5 - cos(C1\_3*acos(1.0 - C24\_a*(Vol\_open - vol)) - C2pi\_3)}
898               tmp\_val\_m1\_to\_p1 = 1.0 - c24\_a*(vol\_open - vol)
899               tmp\_val\_m1\_to\_p1 = max(-1., min(1., tmp\_val\_m1\_to\_p1))
900               l(k) = 0.5 - cos(c1\_3*acos(tmp\_val\_m1\_to\_p1) - c2pi\_3)
901               \textcolor{comment}{! To check the answers.}
902               \textcolor{comment}{! Vol\_err = Vol\_open - a\_12*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol}
903 \textcolor{keywordflow}{            endif}
904           \textcolor{keywordflow}{else} \textcolor{comment}{! a < 0.}
905             \textcolor{keywordflow}{if} (vol <= vol\_direct) \textcolor{keywordflow}{then}
906               \textcolor{comment}{! Both edges of the cell are bounded by walls.}
907               l(k) = (-0.25*c24\_a*vol)**c1\_3
908             \textcolor{keywordflow}{else}
909               \textcolor{comment}{! x\_R is at 1/2 but x\_L is in the interior & L is found by solving}
910               \textcolor{comment}{!   vol = 0.5*L^2*(slope + a/3*(3-4L))}
911 
912               \textcolor{comment}{!  Vol\_err = 0.5*(L(K+1)*L(K+1))*(slope + a\_3*(3.0-4.0*L(K+1))) - vol\_below}
913               \textcolor{comment}{! Change to ...}
914               \textcolor{comment}{!   if (min(Vol\_below + Vol\_err, vol) <= Vol\_direct) then ?}
915               \textcolor{keywordflow}{if} (vol\_below + vol\_err <= vol\_direct) \textcolor{keywordflow}{then}
916                 l0 = l\_direct ; vol\_0 = vol\_direct
917               \textcolor{keywordflow}{else}
918                 l0 = l(k+1) ; vol\_0 = vol\_below + vol\_err
919                 \textcolor{comment}{! Change to   Vol\_0 = min(Vol\_below + Vol\_err, vol) ?}
920 \textcolor{keywordflow}{              endif}
921 
922               \textcolor{comment}{!   Try a relatively simple solution that usually works well}
923               \textcolor{comment}{! for massless layers.}
924               dv\_dl2 = 0.5*(slope+a) - a*l0 ; dvol = (vol-vol\_0)
925            \textcolor{comment}{!  dV\_dL2 = 0.5*(slope+a) - a*L0 ; dVol = max(vol-Vol\_0, 0.0)}
926 
927               use\_l0 = .false.
928               do\_one\_l\_iter = .false.
929               \textcolor{keywordflow}{if} (cs%answers\_2018) \textcolor{keywordflow}{then}
930                 curv\_tol = gv%Angstrom\_H*dv\_dl2**2 &
931                            * (0.25 * dv\_dl2 * gv%Angstrom\_H - a * l0 * dvol)
932                 do\_one\_l\_iter = (a * a * dvol**3) < curv\_tol
933               \textcolor{keywordflow}{else}
934                 \textcolor{comment}{! The following code is more robust when GV%Angstrom\_H=0, but}
935                 \textcolor{comment}{! it changes answers.}
936                 use\_l0 = (dvol <= 0.)
937 
938                 vol\_tol = max(0.5 * gv%Angstrom\_H + gv%H\_subroundoff, 1e-14 * vol)
939                 vol\_quit = max(0.9 * gv%Angstrom\_H + gv%H\_subroundoff, 1e-14 * vol)
940 
941                 curv\_tol = vol\_tol * dv\_dl2**2 &
942                            * (dv\_dl2 * vol\_tol - 2.0 * a * l0 * dvol)
943                 do\_one\_l\_iter = (a * a * dvol**3) < curv\_tol
944 \textcolor{keywordflow}{              endif}
945 
946               \textcolor{keywordflow}{if} (use\_l0) \textcolor{keywordflow}{then}
947                 l(k) = l0
948                 vol\_err = 0.5*(l(k)*l(k))*(slope + a\_3*(3.0-4.0*l(k))) - vol
949               \textcolor{keywordflow}{elseif} (do\_one\_l\_iter) \textcolor{keywordflow}{then}
950                 \textcolor{comment}{! One iteration of Newton's method should give an estimate}
951                 \textcolor{comment}{! that is accurate to within Vol\_tol.}
952                 l(k) = sqrt(l0*l0 + dvol / dv\_dl2)
953                 vol\_err = 0.5*(l(k)*l(k))*(slope + a\_3*(3.0-4.0*l(k))) - vol
954               \textcolor{keywordflow}{else}
955                 \textcolor{keywordflow}{if} (dv\_dl2*(1.0-l0*l0) < dvol + &
956                     dv\_dl2 * (vol\_open - vol)*ibma\_2) \textcolor{keywordflow}{then}
957                   l\_max = sqrt(1.0 - (vol\_open - vol)*ibma\_2)
958                 \textcolor{keywordflow}{else}
959                   l\_max = sqrt(l0*l0 + dvol / dv\_dl2)
960 \textcolor{keywordflow}{                endif}
961                 l\_min = sqrt(l0*l0 + dvol / (0.5*(slope+a) - a*l\_max))
962 
963                 vol\_err\_min = 0.5*(l\_min**2)*(slope + a\_3*(3.0-4.0*l\_min)) - vol
964                 vol\_err\_max = 0.5*(l\_max**2)*(slope + a\_3*(3.0-4.0*l\_max)) - vol
965            \textcolor{comment}{!    if ((abs(Vol\_err\_min) <= Vol\_quit) .or. (Vol\_err\_min >= Vol\_err\_max)) then}
966                 \textcolor{keywordflow}{if} (abs(vol\_err\_min) <= vol\_quit) \textcolor{keywordflow}{then}
967                   l(k) = l\_min ; vol\_err = vol\_err\_min
968                 \textcolor{keywordflow}{else}
969                   l(k) = sqrt((l\_min**2*vol\_err\_max - l\_max**2*vol\_err\_min) / &
970                               (vol\_err\_max - vol\_err\_min))
971                   \textcolor{keywordflow}{do} itt=1,maxitt
972                     vol\_err = 0.5*(l(k)*l(k))*(slope + a\_3*(3.0-4.0*l(k))) - vol
973                     \textcolor{keywordflow}{if} (abs(vol\_err) <= vol\_quit) \textcolor{keywordflow}{exit}
974                     \textcolor{comment}{! Take a Newton's method iteration. This equation has proven}
975                     \textcolor{comment}{! robust enough not to need bracketing.}
976                     l(k) = l(k) - vol\_err / (l(k)* (slope + a - 2.0*a*l(k)))
977                     \textcolor{comment}{! This would be a Newton's method iteration for L^2:}
978                     \textcolor{comment}{!   L(K) = sqrt(L(K)*L(K) - Vol\_err / (0.5*(slope+a) - a*L(K)))}
979 \textcolor{keywordflow}{                  enddo}
980 \textcolor{keywordflow}{                endif} \textcolor{comment}{! end of iterative solver}
981 \textcolor{keywordflow}{              endif} \textcolor{comment}{! end of 1-boundary alternatives.}
982 \textcolor{keywordflow}{            endif} \textcolor{comment}{! end of a<0 cases.}
983 \textcolor{keywordflow}{          endif}
984 
985           \textcolor{comment}{! Determine the drag contributing to the bottom boundary layer}
986           \textcolor{comment}{! and the Raleigh drag that acts on each layer.}
987           \textcolor{keywordflow}{if} (l(k) > l(k+1)) \textcolor{keywordflow}{then}
988             \textcolor{keywordflow}{if} (vol\_below < bbl\_thick) \textcolor{keywordflow}{then}
989               bbl\_frac = (1.0-vol\_below/bbl\_thick)**2
990               bbl\_visc\_frac = bbl\_visc\_frac + bbl\_frac*(l(k) - l(k+1))
991             \textcolor{keywordflow}{else}
992               bbl\_frac = 0.0
993 \textcolor{keywordflow}{            endif}
994 
995             \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then} ; cell\_width = g%dy\_Cu(i,j)
996             \textcolor{keywordflow}{else} ; cell\_width = g%dx\_Cv(i,j) ;\textcolor{keywordflow}{ endif}
997             gam = 1.0 - l(k+1)/l(k)
998             rayleigh = us%L\_to\_Z * cs%cdrag * (l(k)-l(k+1)) * (1.0-bbl\_frac) * &
999                 (12.0*cs%c\_Smag*h\_vel\_pos) /  (12.0*cs%c\_Smag*h\_vel\_pos + &
1000                  us%L\_to\_Z*gv%Z\_to\_H * cs%cdrag * gam*(1.0-gam)*(1.0-1.5*gam) * l(k)**2 * cell\_width)
1001           \textcolor{keywordflow}{else} \textcolor{comment}{! This layer feels no drag.}
1002             rayleigh = 0.0
1003 \textcolor{keywordflow}{          endif}
1004 
1005           \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then}
1006             \textcolor{keywordflow}{if} (rayleigh > 0.0) \textcolor{keywordflow}{then}
1007               v\_at\_u = set\_v\_at\_u(v, h, g, i, j, k, mask\_v, obc)
1008               visc%Ray\_u(i,j,k) = rayleigh*sqrt(u(i,j,k)*u(i,j,k) + &
1009                                                 v\_at\_u*v\_at\_u + u\_bg\_sq)
1010             \textcolor{keywordflow}{else} ; visc%Ray\_u(i,j,k) = 0.0 ;\textcolor{keywordflow}{ endif}
1011           \textcolor{keywordflow}{else}
1012             \textcolor{keywordflow}{if} (rayleigh > 0.0) \textcolor{keywordflow}{then}
1013               u\_at\_v = set\_u\_at\_v(u, h, g, i, j, k, mask\_u, obc)
1014               visc%Ray\_v(i,j,k) = rayleigh*sqrt(v(i,j,k)*v(i,j,k) + &
1015                                                 u\_at\_v*u\_at\_v + u\_bg\_sq)
1016             \textcolor{keywordflow}{else} ; visc%Ray\_v(i,j,k) = 0.0 ;\textcolor{keywordflow}{ endif}
1017 \textcolor{keywordflow}{          endif}
1018 
1019 \textcolor{keywordflow}{        enddo} \textcolor{comment}{! k loop to determine L(K).}
1020 
1021         \textcolor{comment}{! Set the near-bottom viscosity to a value which will give}
1022         \textcolor{comment}{! the correct stress when the shear occurs over bbl\_thick.}
1023         \textcolor{comment}{! See next block for explanation.}
1024         bbl\_thick\_z = bbl\_thick * gv%H\_to\_Z
1025         \textcolor{keywordflow}{if} (cs%correct\_BBL\_bounds .and. &
1026             cdrag\_sqrt*ustar(i)*bbl\_thick\_z*bbl\_visc\_frac <= cs%Kv\_BBL\_min) \textcolor{keywordflow}{then}
1027           \textcolor{comment}{! If the bottom stress implies less viscosity than Kv\_BBL\_min then}
1028           \textcolor{comment}{! set kv\_bbl to the bound and recompute bbl\_thick to be consistent}
1029           \textcolor{comment}{! but with a ridiculously large upper bound on thickness (for Cd u*=0)}
1030           kv\_bbl = cs%Kv\_BBL\_min
1031           \textcolor{keywordflow}{if} (cdrag\_sqrt*ustar(i)*bbl\_visc\_frac*g%Rad\_Earth*us%m\_to\_Z > kv\_bbl) \textcolor{keywordflow}{then}
1032             bbl\_thick\_z = kv\_bbl / ( cdrag\_sqrt*ustar(i)*bbl\_visc\_frac )
1033           \textcolor{keywordflow}{else}
1034             bbl\_thick\_z = g%Rad\_Earth * us%m\_to\_Z
1035 \textcolor{keywordflow}{          endif}
1036         \textcolor{keywordflow}{else}
1037           kv\_bbl = cdrag\_sqrt*ustar(i)*bbl\_thick\_z*bbl\_visc\_frac
1038 \textcolor{keywordflow}{        endif}
1039 
1040       \textcolor{keywordflow}{else} \textcolor{comment}{! Not Channel\_drag.}
1041         \textcolor{comment}{! Set the near-bottom viscosity to a value which will give}
1042         \textcolor{comment}{! the correct stress when the shear occurs over bbl\_thick.}
1043         \textcolor{comment}{! - The bottom stress is tau\_b = Cdrag * u\_bbl^2}
1044         \textcolor{comment}{! - u\_bbl was calculated by averaging flow over CS%Hbbl}
1045         \textcolor{comment}{!   (and includes unresolved tidal components)}
1046         \textcolor{comment}{! - u\_bbl is embedded in u* since u*^2 = Cdrag u\_bbl^2}
1047         \textcolor{comment}{! - The average shear in the BBL is du/dz = 2 * u\_bbl / h\_bbl}
1048         \textcolor{comment}{!   (which assumes a linear profile, hence the "2")}
1049         \textcolor{comment}{! - bbl\_thick was bounded to <= 0.5 * CS%Hbbl}
1050         \textcolor{comment}{! - The viscous stress kv\_bbl du/dz should balance tau\_b}
1051         \textcolor{comment}{!      Cdrag u\_bbl^2 = kv\_bbl du/dz}
1052         \textcolor{comment}{!                    = 2 kv\_bbl u\_bbl}
1053         \textcolor{comment}{! so}
1054         \textcolor{comment}{!      kv\_bbl = 0.5 h\_bbl Cdrag u\_bbl}
1055         \textcolor{comment}{!             = 0.5 h\_bbl sqrt(Cdrag) u*}
1056         bbl\_thick\_z = bbl\_thick * gv%H\_to\_Z
1057         \textcolor{keywordflow}{if} (cs%correct\_BBL\_bounds .and. &
1058             cdrag\_sqrt*ustar(i)*bbl\_thick\_z <= cs%Kv\_BBL\_min) \textcolor{keywordflow}{then}
1059           \textcolor{comment}{! If the bottom stress implies less viscosity than Kv\_BBL\_min then}
1060           \textcolor{comment}{! set kv\_bbl to the bound and recompute bbl\_thick to be consistent}
1061           \textcolor{comment}{! but with a ridiculously large upper bound on thickness (for Cd u*=0)}
1062           kv\_bbl = cs%Kv\_BBL\_min
1063           \textcolor{keywordflow}{if} (cdrag\_sqrt*ustar(i)*g%Rad\_Earth*us%m\_to\_Z > kv\_bbl) \textcolor{keywordflow}{then}
1064             bbl\_thick\_z = kv\_bbl / ( cdrag\_sqrt*ustar(i) )
1065           \textcolor{keywordflow}{else}
1066             bbl\_thick\_z = g%Rad\_Earth * us%m\_to\_Z
1067 \textcolor{keywordflow}{          endif}
1068         \textcolor{keywordflow}{else}
1069           kv\_bbl = cdrag\_sqrt*ustar(i)*bbl\_thick\_z
1070 \textcolor{keywordflow}{        endif}
1071 \textcolor{keywordflow}{      endif}
1072       kv\_bbl = max(cs%Kv\_BBL\_min, kv\_bbl)
1073       \textcolor{keywordflow}{if} (m==1) \textcolor{keywordflow}{then}
1074         visc%Kv\_bbl\_u(i,j) = kv\_bbl
1075         visc%bbl\_thick\_u(i,j) = bbl\_thick\_z
1076       \textcolor{keywordflow}{else}
1077         visc%Kv\_bbl\_v(i,j) = kv\_bbl
1078         visc%bbl\_thick\_v(i,j) = bbl\_thick\_z
1079 \textcolor{keywordflow}{      endif}
1080 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end of i loop}
1081 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end of m & j loops}
1082 
1083 \textcolor{comment}{! Offer diagnostics for averaging}
1084   \textcolor{keywordflow}{if} (cs%id\_bbl\_thick\_u > 0) &
1085     \textcolor{keyword}{call }post\_data(cs%id\_bbl\_thick\_u, visc%bbl\_thick\_u, cs%diag)
1086   \textcolor{keywordflow}{if} (cs%id\_kv\_bbl\_u > 0) &
1087     \textcolor{keyword}{call }post\_data(cs%id\_kv\_bbl\_u, visc%kv\_bbl\_u, cs%diag)
1088   \textcolor{keywordflow}{if} (cs%id\_bbl\_u > 0) &
1089     \textcolor{keyword}{call }post\_data(cs%id\_bbl\_u, cs%bbl\_u, cs%diag)
1090   \textcolor{keywordflow}{if} (cs%id\_bbl\_thick\_v > 0) &
1091     \textcolor{keyword}{call }post\_data(cs%id\_bbl\_thick\_v, visc%bbl\_thick\_v, cs%diag)
1092   \textcolor{keywordflow}{if} (cs%id\_kv\_bbl\_v > 0) &
1093     \textcolor{keyword}{call }post\_data(cs%id\_kv\_bbl\_v, visc%kv\_bbl\_v, cs%diag)
1094   \textcolor{keywordflow}{if} (cs%id\_bbl\_v > 0) &
1095     \textcolor{keyword}{call }post\_data(cs%id\_bbl\_v, cs%bbl\_v, cs%diag)
1096   \textcolor{keywordflow}{if} (cs%id\_Ray\_u > 0) &
1097     \textcolor{keyword}{call }post\_data(cs%id\_Ray\_u, visc%Ray\_u, cs%diag)
1098   \textcolor{keywordflow}{if} (cs%id\_Ray\_v > 0) &
1099     \textcolor{keyword}{call }post\_data(cs%id\_Ray\_v, visc%Ray\_v, cs%diag)
1100 
1101   \textcolor{keywordflow}{if} (cs%debug) \textcolor{keywordflow}{then}
1102     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Ray\_u) .and. \textcolor{keyword}{associated}(visc%Ray\_v)) &
1103         \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"Ray [uv]"}, visc%Ray\_u, visc%Ray\_v, g%HI, haloshift=0, scale=us%Z\_to\_m*us%s\_to\_T)
1104     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%kv\_bbl\_u) .and. \textcolor{keyword}{associated}(visc%kv\_bbl\_v)) &
1105         \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"kv\_bbl\_[uv]"}, visc%kv\_bbl\_u, visc%kv\_bbl\_v, g%HI, &
1106                       haloshift=0, scale=us%Z2\_T\_to\_m2\_s, scalar\_pair=.true.)
1107     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%bbl\_thick\_u) .and. \textcolor{keyword}{associated}(visc%bbl\_thick\_v)) &
1108         \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"bbl\_thick\_[uv]"}, visc%bbl\_thick\_u, visc%bbl\_thick\_v, &
1109                       g%HI, haloshift=0, scale=us%Z\_to\_m, scalar\_pair=.true.)
1110 \textcolor{keywordflow}{  endif}
1111 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__set__visc_aba41cd4f8baa1cda9036d97087ce8a22}\label{namespacemom__set__visc_aba41cd4f8baa1cda9036d97087ce8a22}} 
\index{mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}!set\+\_\+viscous\+\_\+ml@{set\+\_\+viscous\+\_\+ml}}
\index{set\+\_\+viscous\+\_\+ml@{set\+\_\+viscous\+\_\+ml}!mom\+\_\+set\+\_\+visc@{mom\+\_\+set\+\_\+visc}}
\subsubsection{\texorpdfstring{set\+\_\+viscous\+\_\+ml()}{set\_viscous\_ml()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+set\+\_\+visc\+::set\+\_\+viscous\+\_\+ml (\begin{DoxyParamCaption}\item[{real, dimension( g \%isdb\+: g \%iedb, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{u,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsdb\+: g \%jedb, g \%ke), intent(in)}]{v,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, g \%ke), intent(in)}]{h,  }\item[{type(thermo\+\_\+var\+\_\+ptrs), intent(in)}]{tv,  }\item[{type(mech\+\_\+forcing), intent(in)}]{forces,  }\item[{type(vertvisc\+\_\+type), intent(inout)}]{visc,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(inout)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\hyperlink{structmom__set__visc_1_1set__visc__cs}{set\+\_\+visc\+\_\+cs}), pointer}]{CS,  }\item[{logical, intent(in), optional}]{symmetrize }\end{DoxyParamCaption})}



Calculates the thickness of the surface boundary layer for applying an elevated viscosity. 

A bulk Richardson criterion or the thickness of the topmost N\+K\+ML layers (with a bulk mixed layer) are currently used. The thicknesses are given in terms of fractional layers, so that this thickness will move as the thickness of the topmost layers change.


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in,out}  & {\em g} & The ocean\textquotesingle{}s grid structure.\\
\hline
\mbox{\tt in}  & {\em gv} & The ocean\textquotesingle{}s vertical grid structure.\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em u} & The zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em v} & The meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em h} & Layer thicknesses \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em tv} & A structure containing pointers to any available thermodynamic fields. Absent fields have N\+U\+LL ptrs.\\
\hline
\mbox{\tt in}  & {\em forces} & A structure with the driving mechanical forces\\
\hline
\mbox{\tt in,out}  & {\em visc} & A structure containing vertical viscosities and related fields.\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}.\\
\hline
 & {\em cs} & The control structure returned by a previous call to set\+\_\+visc\+\_\+init.\\
\hline
\mbox{\tt in}  & {\em symmetrize} & If present and true, do extra calculations of those values in visc that would be calculated with symmetric memory. \\
\hline
\end{DoxyParams}


Definition at line 1208 of file M\+O\+M\+\_\+set\+\_\+viscosity.\+F90.


\begin{DoxyCode}
1208   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(inout)} :: g\textcolor{comment}{    !< The ocean's grid structure.}
1209   \textcolor{keywordtype}{type}(verticalgrid\_type), \textcolor{keywordtype}{intent(in)}    :: gv\textcolor{comment}{   !< The ocean's vertical grid structure.}
1210   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: us\textcolor{comment}{   !< A dimensional unit scaling type}
1211   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(G))}, &
1212                            \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{    !< The zonal velocity [L T-1 ~> m s-1].}
1213   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(G))}, &
1214                            \textcolor{keywordtype}{intent(in)}    :: v\textcolor{comment}{    !< The meridional velocity [L T-1 ~> m s-1].}
1215   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(G))},  &
1216                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{    !< Layer thicknesses [H ~> m or kg m-2].}
1217   \textcolor{keywordtype}{type}(thermo\_var\_ptrs),   \textcolor{keywordtype}{intent(in)}    :: tv\textcolor{comment}{   !< A structure containing pointers to any available}
1218 \textcolor{comment}{                                                 !! thermodynamic fields. Absent fields have}
1219 \textcolor{comment}{                                                 !! NULL ptrs.}
1220   \textcolor{keywordtype}{type}(mech\_forcing),      \textcolor{keywordtype}{intent(in)}    :: forces\textcolor{comment}{ !< A structure with the driving mechanical forces}
1221   \textcolor{keywordtype}{type}(vertvisc\_type),     \textcolor{keywordtype}{intent(inout)} :: visc\textcolor{comment}{ !< A structure containing vertical viscosities and}
1222 \textcolor{comment}{                                                 !! related fields.}
1223   \textcolor{keywordtype}{real},                    \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{   !< Time increment [T ~> s].}
1224   \textcolor{keywordtype}{type}(set\_visc\_cs),       \textcolor{keywordtype}{pointer}       :: cs\textcolor{comment}{   !< The control structure returned by a previous}
1225 \textcolor{comment}{                                                 !! call to set\_visc\_init.}
1226   \textcolor{keywordtype}{logical},        \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: symmetrize\textcolor{comment}{ !< If present and true, do extra calculations}
1227 \textcolor{comment}{                                                 !! of those values in visc that would be}
1228 \textcolor{comment}{                                                 !! calculated with symmetric memory.}
1229 
1230   \textcolor{comment}{! Local variables}
1231   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
1232     htot, &     \textcolor{comment}{!   The total depth of the layers being that are within the}
1233                 \textcolor{comment}{! surface mixed layer [H ~> m or kg m-2].}
1234     thtot, &    \textcolor{comment}{!   The integrated temperature of layers that are within the}
1235                 \textcolor{comment}{! surface mixed layer [H degC ~> m degC or kg degC m-2].}
1236     shtot, &    \textcolor{comment}{!   The integrated salt of layers that are within the}
1237                 \textcolor{comment}{! surface mixed layer [H ppt ~> m ppt or kg ppt m-2].}
1238     rhtot, &    \textcolor{comment}{!   The integrated density of layers that are within the surface mixed layer}
1239                 \textcolor{comment}{! [H R ~> kg m-2 or kg2 m-5].  Rhtot is only used if no}
1240                 \textcolor{comment}{! equation of state is used.}
1241     uhtot, &    \textcolor{comment}{!   The depth integrated zonal and meridional velocities within}
1242     vhtot, &    \textcolor{comment}{! the surface mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].}
1243     idecay\_len\_tke, & \textcolor{comment}{! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].}
1244     dr\_dt, &    \textcolor{comment}{!   Partial derivative of the density at the base of layer nkml}
1245                 \textcolor{comment}{! (roughly the base of the mixed layer) with temperature [R degC-1 ~> kg m-3 degC-1].}
1246     dr\_ds, &    \textcolor{comment}{!   Partial derivative of the density at the base of layer nkml}
1247                 \textcolor{comment}{! (roughly the base of the mixed layer) with salinity [R ppt-1 ~> kg m-3 ppt-1].}
1248     ustar, &    \textcolor{comment}{!   The surface friction velocity under ice shelves [Z T-1 ~> m s-1].}
1249     press, &    \textcolor{comment}{! The pressure at which dR\_dT and dR\_dS are evaluated [R L2 T-2 ~> Pa].}
1250     t\_eos, &    \textcolor{comment}{! The potential temperature at which dR\_dT and dR\_dS are evaluated [degC]}
1251     s\_eos       \textcolor{comment}{! The salinity at which dR\_dT and dR\_dS are evaluated [ppt].}
1252   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G))} :: &
1253     mask\_u      \textcolor{comment}{! A mask that disables any contributions from u points that}
1254                 \textcolor{comment}{! are land or past open boundary conditions [nondim], 0 or 1.}
1255   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G))} :: &
1256     mask\_v      \textcolor{comment}{! A mask that disables any contributions from v points that}
1257                 \textcolor{comment}{! are land or past open boundary conditions [nondim], 0 or 1.}
1258   \textcolor{keywordtype}{real} :: h\_at\_vel(szib\_(g),szk\_(g))\textcolor{comment}{! Layer thickness at velocity points,}
1259                 \textcolor{comment}{! using an upwind-biased second order accurate estimate based}
1260                 \textcolor{comment}{! on the previous velocity direction [H ~> m or kg m-2].}
1261   \textcolor{keywordtype}{integer} :: k\_massive(szib\_(g)) \textcolor{comment}{! The k-index of the deepest layer yet found}
1262                 \textcolor{comment}{! that has more than h\_tiny thickness and will be in the}
1263                 \textcolor{comment}{! viscous mixed layer.}
1264   \textcolor{keywordtype}{real} :: uh2   \textcolor{comment}{! The squared magnitude of the difference between the velocity}
1265                 \textcolor{comment}{! integrated through the mixed layer and the velocity of the}
1266                 \textcolor{comment}{! interior layer layer times the depth of the the mixed layer}
1267                 \textcolor{comment}{! [H2 Z2 T-2 ~> m4 s-2 or kg2 m-2 s-2].}
1268   \textcolor{keywordtype}{real} :: htot\_vel  \textcolor{comment}{! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].}
1269   \textcolor{keywordtype}{real} :: hwtot     \textcolor{comment}{! Sum of the thicknesses used to calculate}
1270                     \textcolor{comment}{! the near-bottom velocity magnitude [H ~> m or kg m-2].}
1271   \textcolor{keywordtype}{real} :: hutot     \textcolor{comment}{! Running sum of thicknesses times the}
1272                     \textcolor{comment}{! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].}
1273   \textcolor{keywordtype}{real} :: hweight   \textcolor{comment}{! The thickness of a layer that is within Hbbl}
1274                     \textcolor{comment}{! of the bottom [H ~> m or kg m-2].}
1275   \textcolor{keywordtype}{real} :: tbl\_thick\_z  \textcolor{comment}{! The thickness of the top boundary layer [Z ~> m].}
1276 
1277   \textcolor{keywordtype}{real} :: hlay      \textcolor{comment}{! The layer thickness at velocity points [H ~> m or kg m-2].}
1278   \textcolor{keywordtype}{real} :: i\_2hlay   \textcolor{comment}{! 1 / 2*hlay [H-1 ~> m-1 or m2 kg-1].}
1279   \textcolor{keywordtype}{real} :: t\_lay     \textcolor{comment}{! The layer temperature at velocity points [degC].}
1280   \textcolor{keywordtype}{real} :: s\_lay     \textcolor{comment}{! The layer salinity at velocity points [ppt].}
1281   \textcolor{keywordtype}{real} :: rlay      \textcolor{comment}{! The layer potential density at velocity points [R ~> kg m-3].}
1282   \textcolor{keywordtype}{real} :: rlb       \textcolor{comment}{! The potential density of the layer below [R ~> kg m-3].}
1283   \textcolor{keywordtype}{real} :: v\_at\_u    \textcolor{comment}{! The meridonal velocity at a zonal velocity point [L T-1 ~> m s-1].}
1284   \textcolor{keywordtype}{real} :: u\_at\_v    \textcolor{comment}{! The zonal velocity at a meridonal velocity point [L T-1 ~> m s-1].}
1285   \textcolor{keywordtype}{real} :: ghprime   \textcolor{comment}{! The mixed-layer internal gravity wave speed squared, based}
1286                     \textcolor{comment}{! on the mixed layer thickness and density difference across}
1287                     \textcolor{comment}{! the base of the mixed layer [L2 T-2 ~> m2 s-2].}
1288   \textcolor{keywordtype}{real} :: ribulk    \textcolor{comment}{! The bulk Richardson number below which water is in the}
1289                     \textcolor{comment}{! viscous mixed layer, including reduction for turbulent}
1290                     \textcolor{comment}{! decay. Nondimensional.}
1291   \textcolor{keywordtype}{real} :: dt\_rho0   \textcolor{comment}{! The time step divided by the conversion from the layer}
1292                     \textcolor{comment}{! thickness to layer mass [T H Z-1 R-1 ~> s m3 kg-1 or s].}
1293   \textcolor{keywordtype}{real} :: g\_h\_rho0  \textcolor{comment}{!   The gravitational acceleration times the conversion from H to m divided}
1294                     \textcolor{comment}{! by the mean density [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].}
1295   \textcolor{keywordtype}{real} :: ustarsq     \textcolor{comment}{! 400 times the square of ustar, times}
1296                       \textcolor{comment}{! Rho0 divided by G\_Earth and the conversion}
1297                       \textcolor{comment}{! from m to thickness units [H R ~> kg m-2 or kg2 m-5].}
1298   \textcolor{keywordtype}{real} :: cdrag\_sqrt\_z  \textcolor{comment}{! Square root of the drag coefficient, times a unit conversion}
1299                       \textcolor{comment}{! factor from lateral lengths to vertical depths [Z L-1 ~> 1].}
1300   \textcolor{keywordtype}{real} :: cdrag\_sqrt  \textcolor{comment}{! Square root of the drag coefficient [nondim].}
1301   \textcolor{keywordtype}{real} :: oldfn       \textcolor{comment}{! The integrated energy required to}
1302                       \textcolor{comment}{! entrain up to the bottom of the layer,}
1303                       \textcolor{comment}{! divided by G\_Earth [H R ~> kg m-2 or kg2 m-5].}
1304   \textcolor{keywordtype}{real} :: dfn         \textcolor{comment}{! The increment in oldfn for entraining}
1305                       \textcolor{comment}{! the layer [H R ~> kg m-2 or kg2 m-5].}
1306   \textcolor{keywordtype}{real} :: dh          \textcolor{comment}{! The increment in layer thickness from}
1307                       \textcolor{comment}{! the present layer [H ~> m or kg m-2].}
1308   \textcolor{keywordtype}{real} :: u\_bg\_sq   \textcolor{comment}{! The square of an assumed background velocity, for}
1309                     \textcolor{comment}{! calculating the mean magnitude near the top for use in}
1310                     \textcolor{comment}{! the quadratic surface drag [L2 T-2 ~> m2 s-2].}
1311   \textcolor{keywordtype}{real} :: h\_tiny    \textcolor{comment}{! A very small thickness [H ~> m or kg m-2]. Layers that are less than}
1312                     \textcolor{comment}{! h\_tiny can not be the deepest in the viscous mixed layer.}
1313   \textcolor{keywordtype}{real} :: absf      \textcolor{comment}{! The absolute value of f averaged to velocity points [T-1 ~> s-1].}
1314   \textcolor{keywordtype}{real} :: u\_star    \textcolor{comment}{! The friction velocity at velocity points [Z T-1 ~> m s-1].}
1315   \textcolor{keywordtype}{real} :: h\_neglect \textcolor{comment}{! A thickness that is so small it is usually lost}
1316                     \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
1317   \textcolor{keywordtype}{real} :: rho0x400\_g \textcolor{comment}{! 400*Rho0/G\_Earth, times unit conversion factors}
1318                      \textcolor{comment}{! [R T2 H Z-2 ~> kg s2 m-4 or kg2 s2 m-7].}
1319                      \textcolor{comment}{! The 400 is a constant proposed by Killworth and Edwards, 1999.}
1320   \textcolor{keywordtype}{real} :: ustar1    \textcolor{comment}{! ustar [H T-1 ~> m s-1 or kg m-2 s-1]}
1321   \textcolor{keywordtype}{real} :: h2f2      \textcolor{comment}{! (h*2*f)^2 [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]}
1322   \textcolor{keywordtype}{logical} :: use\_eos, do\_any, do\_any\_shelf, do\_i(szib\_(g))
1323   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, k2, nkmb, nkml, n
1324   \textcolor{keywordtype}{type}(ocean\_obc\_type), \textcolor{keywordtype}{pointer} :: obc => null()
1325 
1326   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1327   isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1328   nkmb = gv%nk\_rho\_varies ; nkml = gv%nkml
1329 
1330   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal,\textcolor{stringliteral}{"MOM\_set\_viscosity(visc\_ML): "}//&
1331          \textcolor{stringliteral}{"Module must be initialized before it is used."})
1332   \textcolor{keywordflow}{if} (.not.(cs%dynamic\_viscous\_ML .or. \textcolor{keyword}{associated}(forces%frac\_shelf\_u) .or. &
1333             \textcolor{keyword}{associated}(forces%frac\_shelf\_v)) ) \textcolor{keywordflow}{return}
1334 
1335   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(symmetrize)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (symmetrize) \textcolor{keywordflow}{then}
1336     jsq = js-1 ; isq = is-1
1337 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ endif}
1338 
1339   rho0x400\_g = 400.0*(gv%Rho0/(us%L\_to\_Z**2 * gv%g\_Earth)) * gv%Z\_to\_H
1340   u\_bg\_sq = cs%drag\_bg\_vel * cs%drag\_bg\_vel
1341   cdrag\_sqrt = sqrt(cs%cdrag)
1342   cdrag\_sqrt\_z = us%L\_to\_Z * sqrt(cs%cdrag)
1343 
1344   obc => cs%OBC
1345   use\_eos = \textcolor{keyword}{associated}(tv%eqn\_of\_state)
1346   dt\_rho0 = dt / gv%H\_to\_RZ
1347   h\_neglect = gv%H\_subroundoff
1348   h\_tiny = 2.0*gv%Angstrom\_H + h\_neglect
1349   g\_h\_rho0 = (gv%g\_Earth*gv%H\_to\_Z) / (gv%Rho0)
1350 
1351   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(forces%frac\_shelf\_u) .neqv. \textcolor{keyword}{associated}(forces%frac\_shelf\_v)) &
1352     \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"set\_viscous\_ML: one of forces%frac\_shelf\_u and "}//&
1353                    \textcolor{stringliteral}{"forces%frac\_shelf\_v is associated, but the other is not."})
1354 
1355   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(forces%frac\_shelf\_u)) \textcolor{keywordflow}{then}
1356     \textcolor{comment}{! This configuration has ice shelves, and the appropriate variables need to be}
1357     \textcolor{comment}{! allocated.  If the arrays have already been allocated, these calls do nothing.}
1358     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%tauy\_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1359     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%tbl\_thick\_shelf\_u, g%IsdB, g%IedB, g%jsd, g%jed)
1360     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%tbl\_thick\_shelf\_v, g%isd, g%ied, g%JsdB, g%JedB)
1361     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%kv\_tbl\_shelf\_u, g%IsdB, g%IedB, g%jsd, g%jed)
1362     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%kv\_tbl\_shelf\_v, g%isd, g%ied, g%JsdB, g%JedB)
1363     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%taux\_shelf, g%IsdB, g%IedB, g%jsd, g%jed)
1364     \textcolor{keyword}{call }safe\_alloc\_ptr(visc%tauy\_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1365 
1366     \textcolor{comment}{!  With a linear drag law under shelves, the friction velocity is already known.}
1367 \textcolor{comment}{!    if (CS%linear\_drag) ustar(:) = cdrag\_sqrt\_Z*CS%drag\_bg\_vel}
1368 \textcolor{keywordflow}{  endif}
1369 
1370   \textcolor{comment}{!$OMP parallel do default(shared)}
1371   \textcolor{keywordflow}{do} j=js-1,je ; \textcolor{keywordflow}{do} i=is-1,ie+1
1372     mask\_v(i,j) = g%mask2dCv(i,j)
1373 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
1374   \textcolor{comment}{!$OMP parallel do default(shared)}
1375   \textcolor{keywordflow}{do} j=js-1,je+1 ; \textcolor{keywordflow}{do} i=is-1,ie
1376     mask\_u(i,j) = g%mask2dCu(i,j)
1377 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo}
1378 
1379   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} n=1,obc%number\_of\_segments
1380     \textcolor{comment}{! Now project bottom depths across cell-corner points in the OBCs.  The two}
1381     \textcolor{comment}{! projections have to occur in sequence and can not be combined easily.}
1382     \textcolor{keywordflow}{if} (.not. obc%segment(n)%on\_pe) cycle
1383     \textcolor{comment}{! Use a one-sided projection of bottom depths at OBC points.}
1384     i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
1385     \textcolor{keywordflow}{if} (obc%segment(n)%is\_N\_or\_S .and. (j >= js-1) .and. (j <= je)) \textcolor{keywordflow}{then}
1386       \textcolor{keywordflow}{do} i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
1387         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_n) mask\_u(i,j+1) = 0.0
1388         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_s) mask\_u(i,j) = 0.0
1389 \textcolor{keywordflow}{      enddo}
1390     \textcolor{keywordflow}{elseif} (obc%segment(n)%is\_E\_or\_W .and. (i >= is-1) .and. (i <= je)) \textcolor{keywordflow}{then}
1391       \textcolor{keywordflow}{do} j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
1392         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_e) mask\_v(i+1,j) = 0.0
1393         \textcolor{keywordflow}{if} (obc%segment(n)%direction == obc\_direction\_w) mask\_v(i,j) = 0.0
1394 \textcolor{keywordflow}{      enddo}
1395 \textcolor{keywordflow}{    endif}
1396 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ endif}
1397 
1398   \textcolor{comment}{!$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use\_EOS,dt\_Rho0, &}
1399   \textcolor{comment}{!$OMP                                     h\_neglect,h\_tiny,g\_H\_Rho0,js,je,OBC,Isq,Ieq,nz,  &}
1400   \textcolor{comment}{!$OMP                                     U\_bg\_sq,mask\_v,cdrag\_sqrt,cdrag\_sqrt\_Z,Rho0x400\_G,nkml)}
1401   \textcolor{keywordflow}{do} j=js,je  \textcolor{comment}{! u-point loop}
1402     \textcolor{keywordflow}{if} (cs%dynamic\_viscous\_ML) \textcolor{keywordflow}{then}
1403       do\_any = .false.
1404       \textcolor{keywordflow}{do} i=isq,ieq
1405         htot(i) = 0.0
1406         \textcolor{keywordflow}{if} (g%mask2dCu(i,j) < 0.5) \textcolor{keywordflow}{then}
1407           do\_i(i) = .false. ; visc%nkml\_visc\_u(i,j) = nkml
1408         \textcolor{keywordflow}{else}
1409           do\_i(i) = .true. ; do\_any = .true.
1410           k\_massive(i) = nkml
1411           thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1412           uhtot(i) = dt\_rho0 * forces%taux(i,j)
1413           vhtot(i) = 0.25 * dt\_rho0 * ((forces%tauy(i,j) + forces%tauy(i+1,j-1)) + &
1414                                        (forces%tauy(i,j-1) + forces%tauy(i+1,j)))
1415 
1416           \textcolor{keywordflow}{if} (cs%omega\_frac >= 1.0) \textcolor{keywordflow}{then} ; absf = 2.0*cs%omega ; \textcolor{keywordflow}{else}
1417             absf = 0.5*(abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i,j-1)))
1418             \textcolor{keywordflow}{if} (cs%omega\_frac > 0.0) &
1419               absf = sqrt(cs%omega\_frac*4.0*cs%omega**2 + (1.0-cs%omega\_frac)*absf**2)
1420 \textcolor{keywordflow}{          endif}
1421           u\_star = max(cs%ustar\_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i+1,j)))
1422           idecay\_len\_tke(i) = ((absf / u\_star) * cs%TKE\_decay) * gv%H\_to\_Z
1423 \textcolor{keywordflow}{        endif}
1424 \textcolor{keywordflow}{      enddo}
1425 
1426       \textcolor{keywordflow}{if} (do\_any) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz
1427         \textcolor{keywordflow}{if} (k > nkml) \textcolor{keywordflow}{then}
1428           do\_any = .false.
1429           \textcolor{keywordflow}{if} (use\_eos .and. (k==nkml+1)) \textcolor{keywordflow}{then}
1430             \textcolor{comment}{! Find dRho/dT and dRho\_dS.}
1431             \textcolor{keywordflow}{do} i=isq,ieq
1432               press(i) = (gv%H\_to\_RZ*gv%g\_Earth) * htot(i)
1433               \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%p\_surf)) press(i) = press(i) + 0.5*(tv%p\_surf(i,j)+tv%p\_surf(i+1,j))
1434               k2 = max(1,nkml)
1435               i\_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h\_neglect)
1436               t\_eos(i) = (h(i,j,k2)*tv%T(i,j,k2) + h(i+1,j,k2)*tv%T(i+1,j,k2)) * i\_2hlay
1437               s\_eos(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i+1,j,k2)*tv%S(i+1,j,k2)) * i\_2hlay
1438 \textcolor{keywordflow}{            enddo}
1439             \textcolor{keyword}{call }calculate\_density\_derivs(t\_eos, s\_eos, press, dr\_dt, dr\_ds, tv%eqn\_of\_state, &
1440                                           (/isq-g%IsdB+1,ieq-g%IsdB+1/) )
1441 \textcolor{keywordflow}{          endif}
1442 
1443           \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1444 
1445             hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
1446             \textcolor{keywordflow}{if} (hlay > h\_tiny) \textcolor{keywordflow}{then} \textcolor{comment}{! Only consider non-vanished layers.}
1447               i\_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
1448               v\_at\_u = 0.5 * (h(i,j,k)   * (v(i,j,k) + v(i,j-1,k)) + &
1449                               h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k))) * i\_2hlay
1450               uh2 = ((uhtot(i) - htot(i)*u(i,j,k))**2 + (vhtot(i) - htot(i)*v\_at\_u)**2)
1451 
1452               \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1453                 t\_lay = (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k)) * i\_2hlay
1454                 s\_lay = (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k)) * i\_2hlay
1455                 ghprime = g\_h\_rho0 * (dr\_dt(i) * (t\_lay*htot(i) - thtot(i)) + &
1456                                       dr\_ds(i) * (s\_lay*htot(i) - shtot(i)))
1457               \textcolor{keywordflow}{else}
1458                 ghprime = g\_h\_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1459 \textcolor{keywordflow}{              endif}
1460 
1461               \textcolor{keywordflow}{if} (ghprime > 0.0) \textcolor{keywordflow}{then}
1462                 ribulk = cs%bulk\_Ri\_ML * exp(-htot(i) * idecay\_len\_tke(i))
1463                 \textcolor{keywordflow}{if} (ribulk * uh2 <= (htot(i)**2) * ghprime) \textcolor{keywordflow}{then}
1464                   visc%nkml\_visc\_u(i,j) = \textcolor{keywordtype}{real}(k\_massive(i))
1465                   do\_i(i) = .false.
1466                 \textcolor{keywordflow}{elseif} (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) \textcolor{keywordflow}{then}
1467                   visc%nkml\_visc\_u(i,j) = \textcolor{keywordtype}{real(k-1)} + &
1468                     ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1469                   do\_i(i) = .false.
1470 \textcolor{keywordflow}{                endif}
1471 \textcolor{keywordflow}{              endif}
1472               k\_massive(i) = k
1473 \textcolor{keywordflow}{            endif} \textcolor{comment}{! hlay > h\_tiny}
1474 
1475             \textcolor{keywordflow}{if} (do\_i(i)) do\_any = .true.
1476 \textcolor{keywordflow}{          endif} ;\textcolor{keywordflow}{ enddo}
1477 
1478           \textcolor{keywordflow}{if} (.not.do\_any) \textcolor{keywordflow}{exit} \textcolor{comment}{! All columns are done.}
1479 \textcolor{keywordflow}{        endif}
1480 
1481         \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1482           htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k))
1483           uhtot(i) = uhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(i,j,k)
1484           vhtot(i) = vhtot(i) + 0.25 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1485                                         h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k)))
1486           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1487             thtot(i) = thtot(i) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k))
1488             shtot(i) = shtot(i) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k))
1489           \textcolor{keywordflow}{else}
1490             rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%Rlay(k)
1491 \textcolor{keywordflow}{          endif}
1492 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo}
1493 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ endif}
1494 
1495       \textcolor{keywordflow}{if} (do\_any) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1496         visc%nkml\_visc\_u(i,j) = k\_massive(i)
1497 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
1498 \textcolor{keywordflow}{    endif} \textcolor{comment}{! dynamic\_viscous\_ML}
1499 
1500     do\_any\_shelf = .false.
1501     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(forces%frac\_shelf\_u)) \textcolor{keywordflow}{then}
1502       \textcolor{keywordflow}{do} i=isq,ieq
1503         \textcolor{keywordflow}{if} (forces%frac\_shelf\_u(i,j)*g%mask2dCu(i,j) == 0.0) \textcolor{keywordflow}{then}
1504           do\_i(i) = .false.
1505           visc%tbl\_thick\_shelf\_u(i,j) = 0.0 ; visc%kv\_tbl\_shelf\_u(i,j) = 0.0
1506         \textcolor{keywordflow}{else}
1507           do\_i(i) = .true. ; do\_any\_shelf = .true.
1508 \textcolor{keywordflow}{        endif}
1509 \textcolor{keywordflow}{      enddo}
1510 \textcolor{keywordflow}{    endif}
1511 
1512     \textcolor{keywordflow}{if} (do\_any\_shelf) \textcolor{keywordflow}{then}
1513       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1514         \textcolor{keywordflow}{if} (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) \textcolor{keywordflow}{then}
1515           h\_at\_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
1516                           (h(i,j,k) + h(i+1,j,k) + h\_neglect)
1517         \textcolor{keywordflow}{else}
1518           h\_at\_vel(i,k) =  0.5 * (h(i,j,k) + h(i+1,j,k))
1519 \textcolor{keywordflow}{        endif}
1520       \textcolor{keywordflow}{else}
1521         h\_at\_vel(i,k) = 0.0 ; ustar(i) = 0.0
1522 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1523 
1524       \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1525         htot\_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1526         thtot(i) = 0.0 ; shtot(i) = 0.0
1527         \textcolor{keywordflow}{if} (use\_eos .or. .not.cs%linear\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz
1528           \textcolor{keywordflow}{if} (htot\_vel>=cs%Htbl\_shelf) \textcolor{keywordflow}{exit} \textcolor{comment}{! terminate the k loop}
1529           hweight = min(cs%Htbl\_shelf - htot\_vel, h\_at\_vel(i,k))
1530           \textcolor{keywordflow}{if} (hweight <= 1.5*gv%Angstrom\_H + h\_neglect) cycle
1531 
1532           htot\_vel  = htot\_vel + h\_at\_vel(i,k)
1533           hwtot = hwtot + hweight
1534 
1535           \textcolor{keywordflow}{if} (.not.cs%linear\_drag) \textcolor{keywordflow}{then}
1536             v\_at\_u = set\_v\_at\_u(v, h, g, i, j, k, mask\_v, obc)
1537             hutot = hutot + hweight * sqrt(u(i,j,k)**2 + &
1538                                            v\_at\_u**2 + u\_bg\_sq)
1539 \textcolor{keywordflow}{          endif}
1540           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1541             thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1542             shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1543 \textcolor{keywordflow}{          endif}
1544 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ endif}
1545 
1546         \textcolor{keywordflow}{if} ((.not.cs%linear\_drag) .and. (hwtot > 0.0)) \textcolor{keywordflow}{then}
1547           ustar(i) = cdrag\_sqrt\_z * hutot/hwtot
1548         \textcolor{keywordflow}{else}
1549           ustar(i) = cdrag\_sqrt\_z * cs%drag\_bg\_vel
1550 \textcolor{keywordflow}{        endif}
1551 
1552         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (hwtot > 0.0) \textcolor{keywordflow}{then}
1553           t\_eos(i) = thtot(i)/hwtot ; s\_eos(i) = shtot(i)/hwtot
1554         \textcolor{keywordflow}{else}
1555           t\_eos(i) = 0.0 ; s\_eos(i) = 0.0
1556 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ endif}
1557 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! I-loop}
1558 
1559       \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1560         \textcolor{keyword}{call }calculate\_density\_derivs(t\_eos, s\_eos, forces%p\_surf(:,j), dr\_dt, dr\_ds, &
1561                                       tv%eqn\_of\_state, (/isq-g%IsdB+1,ieq-g%IsdB+1/) )
1562 \textcolor{keywordflow}{      endif}
1563 
1564       \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1565   \textcolor{comment}{!  The 400.0 in this expression is the square of a constant proposed}
1566   \textcolor{comment}{!  by Killworth and Edwards, 1999, in equation (2.20).}
1567         ustarsq = rho0x400\_g * ustar(i)**2
1568         htot(i) = 0.0
1569         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1570           thtot(i) = 0.0 ; shtot(i) = 0.0
1571           \textcolor{keywordflow}{do} k=1,nz-1
1572             \textcolor{keywordflow}{if} (h\_at\_vel(i,k) <= 0.0) cycle
1573             t\_lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1574             s\_lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1575             oldfn = dr\_dt(i)*(t\_lay*htot(i) - thtot(i)) + dr\_ds(i)*(s\_lay*htot(i) - shtot(i))
1576             \textcolor{keywordflow}{if} (oldfn >= ustarsq) \textcolor{keywordflow}{exit}
1577 
1578             dfn = (dr\_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - t\_lay) + &
1579                    dr\_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - s\_lay)) * &
1580                   (h\_at\_vel(i,k)+htot(i))
1581             \textcolor{keywordflow}{if} ((oldfn + dfn) <= ustarsq) \textcolor{keywordflow}{then}
1582               dh = h\_at\_vel(i,k)
1583             \textcolor{keywordflow}{else}
1584               dh = h\_at\_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1585 \textcolor{keywordflow}{            endif}
1586 
1587             htot(i) = htot(i) + dh
1588             thtot(i) = thtot(i) + t\_lay*dh ; shtot(i) = shtot(i) + s\_lay*dh
1589 \textcolor{keywordflow}{          enddo}
1590           \textcolor{keywordflow}{if} ((oldfn < ustarsq) .and. (h\_at\_vel(i,nz) > 0.0)) \textcolor{keywordflow}{then}
1591             t\_lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
1592             s\_lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
1593             \textcolor{keywordflow}{if} (dr\_dt(i)*(t\_lay*htot(i) - thtot(i)) + &
1594                 dr\_ds(i)*(s\_lay*htot(i) - shtot(i)) < ustarsq) &
1595               htot(i) = htot(i) + h\_at\_vel(i,nz)
1596 \textcolor{keywordflow}{          endif} \textcolor{comment}{! Examination of layer nz.}
1597         \textcolor{keywordflow}{else}  \textcolor{comment}{! Use Rlay as the density variable.}
1598           rhtot = 0.0
1599           \textcolor{keywordflow}{do} k=1,nz-1
1600             rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1601 
1602             oldfn = rlay*htot(i) - rhtot(i)
1603             \textcolor{keywordflow}{if} (oldfn >= ustarsq) \textcolor{keywordflow}{exit}
1604 
1605             dfn = (rlb - rlay)*(h\_at\_vel(i,k)+htot(i))
1606             \textcolor{keywordflow}{if} ((oldfn + dfn) <= ustarsq) \textcolor{keywordflow}{then}
1607               dh = h\_at\_vel(i,k)
1608             \textcolor{keywordflow}{else}
1609               dh = h\_at\_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1610 \textcolor{keywordflow}{            endif}
1611 
1612             htot(i) = htot(i) + dh
1613             rhtot(i) = rhtot(i) + rlay*dh
1614 \textcolor{keywordflow}{          enddo}
1615           \textcolor{keywordflow}{if} (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1616             htot(i) = htot(i) + h\_at\_vel(i,nz)
1617 \textcolor{keywordflow}{        endif} \textcolor{comment}{! use\_EOS}
1618 
1619        \textcolor{comment}{!visc%tbl\_thick\_shelf\_u(I,j) = GV%H\_to\_Z * max(CS%Htbl\_shelf\_min, &}
1620        \textcolor{comment}{!    htot(I) / (0.5 + sqrt(0.25 + &}
1621        \textcolor{comment}{!                 (htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2 / &}
1622        \textcolor{comment}{!                 (ustar(i)*GV%Z\_to\_H)**2 )) )}
1623         ustar1 = ustar(i)*gv%Z\_to\_H
1624         h2f2 = (htot(i)*(g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j)) + h\_neglect*cs%omega)**2
1625         tbl\_thick\_z = gv%H\_to\_Z * max(cs%Htbl\_shelf\_min, &
1626             ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1627         visc%tbl\_thick\_shelf\_u(i,j) = tbl\_thick\_z
1628         visc%Kv\_tbl\_shelf\_u(i,j) = max(cs%Kv\_TBL\_min, cdrag\_sqrt*ustar(i)*tbl\_thick\_z)
1629 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! I-loop}
1630 \textcolor{keywordflow}{    endif} \textcolor{comment}{! do\_any\_shelf}
1631 
1632 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! j-loop at u-points}
1633 
1634   \textcolor{comment}{!$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use\_EOS,dt\_Rho0, &}
1635   \textcolor{comment}{!$OMP                                     h\_neglect,h\_tiny,g\_H\_Rho0,is,ie,OBC,Jsq,Jeq,nz, &}
1636   \textcolor{comment}{!$OMP                                     U\_bg\_sq,cdrag\_sqrt,cdrag\_sqrt\_Z,Rho0x400\_G,nkml,mask\_u)}
1637   \textcolor{keywordflow}{do} j=jsq,jeq  \textcolor{comment}{! v-point loop}
1638     \textcolor{keywordflow}{if} (cs%dynamic\_viscous\_ML) \textcolor{keywordflow}{then}
1639       do\_any = .false.
1640       \textcolor{keywordflow}{do} i=is,ie
1641         htot(i) = 0.0
1642         \textcolor{keywordflow}{if} (g%mask2dCv(i,j) < 0.5) \textcolor{keywordflow}{then}
1643           do\_i(i) = .false. ; visc%nkml\_visc\_v(i,j) = nkml
1644         \textcolor{keywordflow}{else}
1645           do\_i(i) = .true. ; do\_any = .true.
1646           k\_massive(i) = nkml
1647           thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1648           vhtot(i) = dt\_rho0 * forces%tauy(i,j)
1649           uhtot(i) = 0.25 * dt\_rho0 * ((forces%taux(i,j) + forces%taux(i-1,j+1)) + &
1650                                        (forces%taux(i-1,j) + forces%taux(i,j+1)))
1651 
1652          \textcolor{keywordflow}{if} (cs%omega\_frac >= 1.0) \textcolor{keywordflow}{then} ; absf = 2.0*cs%omega ; \textcolor{keywordflow}{else}
1653            absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1654            \textcolor{keywordflow}{if} (cs%omega\_frac > 0.0) &
1655              absf = sqrt(cs%omega\_frac*4.0*cs%omega**2 + (1.0-cs%omega\_frac)*absf**2)
1656 \textcolor{keywordflow}{         endif}
1657 
1658          u\_star = max(cs%ustar\_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i,j+1)))
1659          idecay\_len\_tke(i) = ((absf / u\_star) * cs%TKE\_decay) * gv%H\_to\_Z
1660 
1661 \textcolor{keywordflow}{        endif}
1662 \textcolor{keywordflow}{      enddo}
1663 
1664       \textcolor{keywordflow}{if} (do\_any) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz
1665         \textcolor{keywordflow}{if} (k > nkml) \textcolor{keywordflow}{then}
1666           do\_any = .false.
1667           \textcolor{keywordflow}{if} (use\_eos .and. (k==nkml+1)) \textcolor{keywordflow}{then}
1668             \textcolor{comment}{! Find dRho/dT and dRho\_dS.}
1669             \textcolor{keywordflow}{do} i=is,ie
1670               press(i) = (gv%H\_to\_RZ * gv%g\_Earth) * htot(i)
1671               \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(tv%p\_surf)) press(i) = press(i) + 0.5*(tv%p\_surf(i,j)+tv%p\_surf(i,j+1))
1672               k2 = max(1,nkml)
1673               i\_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h\_neglect)
1674               t\_eos(i) = (h(i,j,k2)*tv%T(i,j,k2) + h(i,j+1,k2)*tv%T(i,j+1,k2)) * i\_2hlay
1675               s\_eos(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i,j+1,k2)*tv%S(i,j+1,k2)) * i\_2hlay
1676 \textcolor{keywordflow}{            enddo}
1677             \textcolor{keyword}{call }calculate\_density\_derivs(t\_eos, s\_eos, press, dr\_dt, dr\_ds, &
1678                                           tv%eqn\_of\_state, (/is-g%IsdB+1,ie-g%IsdB+1/) )
1679 \textcolor{keywordflow}{          endif}
1680 
1681           \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1682 
1683             hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
1684             \textcolor{keywordflow}{if} (hlay > h\_tiny) \textcolor{keywordflow}{then} \textcolor{comment}{! Only consider non-vanished layers.}
1685               i\_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
1686               u\_at\_v = 0.5 * (h(i,j,k)   * (u(i-1,j,k)   + u(i,j,k)) + &
1687                               h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k))) * i\_2hlay
1688               uh2 = ((uhtot(i) - htot(i)*u\_at\_v)**2 + (vhtot(i) - htot(i)*v(i,j,k))**2)
1689 
1690               \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1691                 t\_lay = (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k)) * i\_2hlay
1692                 s\_lay = (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k)) * i\_2hlay
1693                 ghprime = g\_h\_rho0 * (dr\_dt(i) * (t\_lay*htot(i) - thtot(i)) + &
1694                                       dr\_ds(i) * (s\_lay*htot(i) - shtot(i)))
1695               \textcolor{keywordflow}{else}
1696                 ghprime = g\_h\_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1697 \textcolor{keywordflow}{              endif}
1698 
1699               \textcolor{keywordflow}{if} (ghprime > 0.0) \textcolor{keywordflow}{then}
1700                 ribulk = cs%bulk\_Ri\_ML * exp(-htot(i) * idecay\_len\_tke(i))
1701                 \textcolor{keywordflow}{if} (ribulk * uh2 <= htot(i)**2 * ghprime) \textcolor{keywordflow}{then}
1702                   visc%nkml\_visc\_v(i,j) = \textcolor{keywordtype}{real}(k\_massive(i))
1703                   do\_i(i) = .false.
1704                 \textcolor{keywordflow}{elseif} (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) \textcolor{keywordflow}{then}
1705                   visc%nkml\_visc\_v(i,j) = \textcolor{keywordtype}{real(k-1)} + &
1706                     ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1707                   do\_i(i) = .false.
1708 \textcolor{keywordflow}{                endif}
1709 \textcolor{keywordflow}{              endif}
1710               k\_massive(i) = k
1711 \textcolor{keywordflow}{            endif} \textcolor{comment}{! hlay > h\_tiny}
1712 
1713             \textcolor{keywordflow}{if} (do\_i(i)) do\_any = .true.
1714 \textcolor{keywordflow}{          endif} ;\textcolor{keywordflow}{ enddo}
1715 
1716           \textcolor{keywordflow}{if} (.not.do\_any) \textcolor{keywordflow}{exit} \textcolor{comment}{! All columns are done.}
1717 \textcolor{keywordflow}{        endif}
1718 
1719         \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1720           htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k))
1721           vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,j,k)
1722           uhtot(i) = uhtot(i) + 0.25 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1723                                         h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k)))
1724           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1725             thtot(i) = thtot(i) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k))
1726             shtot(i) = shtot(i) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k))
1727           \textcolor{keywordflow}{else}
1728             rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%Rlay(k)
1729 \textcolor{keywordflow}{          endif}
1730 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo}
1731 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ endif}
1732 
1733       \textcolor{keywordflow}{if} (do\_any) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1734         visc%nkml\_visc\_v(i,j) = k\_massive(i)
1735 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
1736 \textcolor{keywordflow}{    endif} \textcolor{comment}{! dynamic\_viscous\_ML}
1737 
1738     do\_any\_shelf = .false.
1739     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(forces%frac\_shelf\_v)) \textcolor{keywordflow}{then}
1740       \textcolor{keywordflow}{do} i=is,ie
1741         \textcolor{keywordflow}{if} (forces%frac\_shelf\_v(i,j)*g%mask2dCv(i,j) == 0.0) \textcolor{keywordflow}{then}
1742           do\_i(i) = .false.
1743           visc%tbl\_thick\_shelf\_v(i,j) = 0.0 ; visc%kv\_tbl\_shelf\_v(i,j) = 0.0
1744         \textcolor{keywordflow}{else}
1745           do\_i(i) = .true. ; do\_any\_shelf = .true.
1746 \textcolor{keywordflow}{        endif}
1747 \textcolor{keywordflow}{      enddo}
1748 \textcolor{keywordflow}{    endif}
1749 
1750     \textcolor{keywordflow}{if} (do\_any\_shelf) \textcolor{keywordflow}{then}
1751       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1752         \textcolor{keywordflow}{if} (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) \textcolor{keywordflow}{then}
1753           h\_at\_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
1754                           (h(i,j,k) + h(i,j+1,k) + h\_neglect)
1755         \textcolor{keywordflow}{else}
1756           h\_at\_vel(i,k) =  0.5 * (h(i,j,k) + h(i,j+1,k))
1757 \textcolor{keywordflow}{        endif}
1758       \textcolor{keywordflow}{else}
1759         h\_at\_vel(i,k) = 0.0 ; ustar(i) = 0.0
1760 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1761 
1762       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1763         htot\_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1764         thtot(i) = 0.0 ; shtot(i) = 0.0
1765         \textcolor{keywordflow}{if} (use\_eos .or. .not.cs%linear\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz
1766           \textcolor{keywordflow}{if} (htot\_vel>=cs%Htbl\_shelf) \textcolor{keywordflow}{exit} \textcolor{comment}{! terminate the k loop}
1767           hweight = min(cs%Htbl\_shelf - htot\_vel, h\_at\_vel(i,k))
1768           \textcolor{keywordflow}{if} (hweight <= 1.5*gv%Angstrom\_H + h\_neglect) cycle
1769 
1770           htot\_vel  = htot\_vel + h\_at\_vel(i,k)
1771           hwtot = hwtot + hweight
1772 
1773           \textcolor{keywordflow}{if} (.not.cs%linear\_drag) \textcolor{keywordflow}{then}
1774             u\_at\_v = set\_u\_at\_v(u, h, g, i, j, k, mask\_u, obc)
1775             hutot = hutot + hweight * sqrt(v(i,j,k)**2 + &
1776                                            u\_at\_v**2 + u\_bg\_sq)
1777 \textcolor{keywordflow}{          endif}
1778           \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1779             thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1780             shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1781 \textcolor{keywordflow}{          endif}
1782 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ endif}
1783 
1784         \textcolor{keywordflow}{if} (.not.cs%linear\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (hwtot > 0.0) \textcolor{keywordflow}{then}
1785           ustar(i) = cdrag\_sqrt\_z * hutot/hwtot
1786         \textcolor{keywordflow}{else}
1787           ustar(i) = cdrag\_sqrt\_z * cs%drag\_bg\_vel
1788 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ endif}
1789 
1790         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (hwtot > 0.0) \textcolor{keywordflow}{then}
1791           t\_eos(i) = thtot(i)/hwtot ; s\_eos(i) = shtot(i)/hwtot
1792         \textcolor{keywordflow}{else}
1793           t\_eos(i) = 0.0 ; s\_eos(i) = 0.0
1794 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ endif}
1795 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! I-loop}
1796 
1797       \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1798         \textcolor{keyword}{call }calculate\_density\_derivs(t\_eos, s\_eos, forces%p\_surf(:,j), dr\_dt, dr\_ds, &
1799                                       tv%eqn\_of\_state, (/is-g%IsdB+1,ie-g%IsdB+1/) )
1800 \textcolor{keywordflow}{      endif}
1801 
1802       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1803   \textcolor{comment}{!  The 400.0 in this expression is the square of a constant proposed}
1804   \textcolor{comment}{!  by Killworth and Edwards, 1999, in equation (2.20).}
1805         ustarsq = rho0x400\_g * ustar(i)**2
1806         htot(i) = 0.0
1807         \textcolor{keywordflow}{if} (use\_eos) \textcolor{keywordflow}{then}
1808           thtot(i) = 0.0 ; shtot(i) = 0.0
1809           \textcolor{keywordflow}{do} k=1,nz-1
1810             \textcolor{keywordflow}{if} (h\_at\_vel(i,k) <= 0.0) cycle
1811             t\_lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1812             s\_lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1813             oldfn = dr\_dt(i)*(t\_lay*htot(i) - thtot(i)) + dr\_ds(i)*(s\_lay*htot(i) - shtot(i))
1814             \textcolor{keywordflow}{if} (oldfn >= ustarsq) \textcolor{keywordflow}{exit}
1815 
1816             dfn = (dr\_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - t\_lay) + &
1817                    dr\_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - s\_lay)) * &
1818                   (h\_at\_vel(i,k)+htot(i))
1819             \textcolor{keywordflow}{if} ((oldfn + dfn) <= ustarsq) \textcolor{keywordflow}{then}
1820               dh = h\_at\_vel(i,k)
1821             \textcolor{keywordflow}{else}
1822               dh = h\_at\_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1823 \textcolor{keywordflow}{            endif}
1824 
1825             htot(i) = htot(i) + dh
1826             thtot(i) = thtot(i) + t\_lay*dh ; shtot(i) = shtot(i) + s\_lay*dh
1827 \textcolor{keywordflow}{          enddo}
1828           \textcolor{keywordflow}{if} ((oldfn < ustarsq) .and. (h\_at\_vel(i,nz) > 0.0)) \textcolor{keywordflow}{then}
1829             t\_lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
1830             s\_lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
1831             \textcolor{keywordflow}{if} (dr\_dt(i)*(t\_lay*htot(i) - thtot(i)) + &
1832                 dr\_ds(i)*(s\_lay*htot(i) - shtot(i)) < ustarsq) &
1833               htot(i) = htot(i) + h\_at\_vel(i,nz)
1834 \textcolor{keywordflow}{          endif} \textcolor{comment}{! Examination of layer nz.}
1835         \textcolor{keywordflow}{else}  \textcolor{comment}{! Use Rlay as the density variable.}
1836           rhtot = 0.0
1837           \textcolor{keywordflow}{do} k=1,nz-1
1838             rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1839 
1840             oldfn = rlay*htot(i) - rhtot(i)
1841             \textcolor{keywordflow}{if} (oldfn >= ustarsq) \textcolor{keywordflow}{exit}
1842 
1843             dfn = (rlb - rlay)*(h\_at\_vel(i,k)+htot(i))
1844             \textcolor{keywordflow}{if} ((oldfn + dfn) <= ustarsq) \textcolor{keywordflow}{then}
1845               dh = h\_at\_vel(i,k)
1846             \textcolor{keywordflow}{else}
1847               dh = h\_at\_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1848 \textcolor{keywordflow}{            endif}
1849 
1850             htot(i) = htot(i) + dh
1851             rhtot = rhtot + rlay*dh
1852 \textcolor{keywordflow}{          enddo}
1853           \textcolor{keywordflow}{if} (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1854             htot(i) = htot(i) + h\_at\_vel(i,nz)
1855 \textcolor{keywordflow}{        endif} \textcolor{comment}{! use\_EOS}
1856 
1857        \textcolor{comment}{!visc%tbl\_thick\_shelf\_v(i,J) = GV%H\_to\_Z * max(CS%Htbl\_shelf\_min, &}
1858        \textcolor{comment}{!    htot(i) / (0.5 + sqrt(0.25 + &}
1859        \textcolor{comment}{!        (htot(i)*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))**2 / &}
1860        \textcolor{comment}{!        (ustar(i)*GV%Z\_to\_H)**2 )) )}
1861         ustar1 = ustar(i)*gv%Z\_to\_H
1862         h2f2 = (htot(i)*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) + h\_neglect*cs%omega)**2
1863         tbl\_thick\_z = gv%H\_to\_Z * max(cs%Htbl\_shelf\_min, &
1864             ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1865         visc%tbl\_thick\_shelf\_v(i,j) = tbl\_thick\_z
1866         visc%Kv\_tbl\_shelf\_v(i,j) = max(cs%Kv\_TBL\_min, cdrag\_sqrt*ustar(i)*tbl\_thick\_z)
1867 
1868 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i-loop}
1869 \textcolor{keywordflow}{    endif} \textcolor{comment}{! do\_any\_shelf}
1870 
1871 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! J-loop at v-points}
1872 
1873   \textcolor{keywordflow}{if} (cs%debug) \textcolor{keywordflow}{then}
1874     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%nkml\_visc\_u) .and. \textcolor{keyword}{associated}(visc%nkml\_visc\_v)) &
1875       \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"nkml\_visc\_[uv]"}, visc%nkml\_visc\_u, visc%nkml\_visc\_v, &
1876                     g%HI, haloshift=0, scalar\_pair=.true.)
1877 \textcolor{keywordflow}{  endif}
1878   \textcolor{keywordflow}{if} (cs%id\_nkml\_visc\_u > 0) &
1879     \textcolor{keyword}{call }post\_data(cs%id\_nkml\_visc\_u, visc%nkml\_visc\_u, cs%diag)
1880   \textcolor{keywordflow}{if} (cs%id\_nkml\_visc\_v > 0) &
1881     \textcolor{keyword}{call }post\_data(cs%id\_nkml\_visc\_v, visc%nkml\_visc\_v, cs%diag)
1882 
\end{DoxyCode}
