\hypertarget{namespacemom__vert__friction}{}\section{mom\+\_\+vert\+\_\+friction Module Reference}
\label{namespacemom__vert__friction}\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}


\subsection{Detailed Description}
Implements vertical viscosity (vertvisc) 

\begin{DoxyAuthor}{Author}
Robert Hallberg 
\end{DoxyAuthor}
\begin{DoxyDate}{Date}
April 1994 -\/ October 2006
\end{DoxyDate}
The vertical diffusion of momentum is fully implicit. This is necessary to allow for vanishingly small layers. The coupling is based on the distance between the centers of adjacent layers, except where a layer is close to the bottom compared with a bottom boundary layer thickness when a bottom drag law is used. A stress top b.\+c. and a no slip bottom b.\+c. are used. There is no limit on the time step for vertvisc.

Near the bottom, the horizontal thickness interpolation scheme changes to an upwind biased estimate to control the effect of spurious Montgomery potential gradients at the bottom where nearly massless layers layers ride over the topography. Within a few boundary layer depths of the bottom, the harmonic mean thickness (i.\+e. (2 h+ h-\/) / (h+ + h-\/) ) is used if the velocity is from the thinner side and the arithmetic mean thickness (i.\+e. (h+ + h-\/)/2) is used if the velocity is from the thicker side. Both of these thickness estimates are second order accurate. Above this the arithmetic mean thickness is used.

In addition, vertvisc truncates any velocity component that exceeds maxvel to truncvel. This basically keeps instabilities spatially localized. The number of times the velocity is truncated is reported each time the energies are saved, and if exceeds CSMaxtrunc the model will stop itself and change the time to a large value. This has proven very useful in (1) diagnosing model failures and (2) letting the model settle down to a meaningful integration from a poorly specified initial condition.

The same code is used for the two velocity components, by indirectly referencing the velocities and defining a handful of direction-\/specific defined variables.

Macros written all in capital letters are defined in \mbox{\hyperlink{MOM__memory_8h}{M\+O\+M\+\_\+memory.\+h}}.

A small fragment of the grid is shown below\+: \begin{DoxyVerb}    j+1  x ^ x ^ x   At x:  q
    j+1  > o > o >   At ^:  v, frhatv, tauy
    j    x ^ x ^ x   At >:  u, frhatu, taux
    j    > o > o >   At o:  h
    j-1  x ^ x ^ x
        i-1  i  i+1  At x & ^:
           i  i+1    At > & o:\end{DoxyVerb}


The boundaries always run through q grid points (x). \subsection*{Data Types}
\begin{DoxyCompactItemize}
\item 
type \mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}
\begin{DoxyCompactList}\small\item\em The control structure with parameters and memory for the M\+O\+M\+\_\+vert\+\_\+friction module. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \mbox{\hyperlink{namespacemom__vert__friction_a8f1a390fa24fbe985068fed9ac26873c}{vertvisc}} (u, v, h, forces, visc, dt, O\+BC, A\+Dp, C\+Dp, G, GV, US, CS, taux\+\_\+bot, tauy\+\_\+bot, Waves)
\begin{DoxyCompactList}\small\item\em Perform a fully implicit vertical diffusion of momentum. Stress top and bottom boundary conditions are used. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__vert__friction_ad0b7bdc6695ee0c4207797bc94672863}{vertvisc\+\_\+remnant}} (visc, visc\+\_\+rem\+\_\+u, visc\+\_\+rem\+\_\+v, dt, G, GV, US, CS)
\begin{DoxyCompactList}\small\item\em Calculate the fraction of momentum originally in a layer that remains after a time-\/step of viscosity, and the fraction of a time-\/step\textquotesingle{}s worth of barotropic acceleration that a layer experiences after viscosity is applied. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__vert__friction_ac281f6595593b33436594112785e982b}{vertvisc\+\_\+coef}} (u, v, h, forces, visc, dt, G, GV, US, CS, O\+BC)
\begin{DoxyCompactList}\small\item\em Calculate the coupling coefficients (CSa\+\_\+u and CSa\+\_\+v) and effective layer thicknesses (CSh\+\_\+u and CSh\+\_\+v) for later use in the applying the implicit vertical viscosity via \mbox{\hyperlink{namespacemom__vert__friction_a8f1a390fa24fbe985068fed9ac26873c}{vertvisc()}}. \end{DoxyCompactList}\item 
subroutine \mbox{\hyperlink{namespacemom__vert__friction_aa9e6f1f0d75a54d85b0d0cdad874b41f}{find\+\_\+coupling\+\_\+coef}} (a\+\_\+cpl, hvel, do\+\_\+i, h\+\_\+harm, bbl\+\_\+thick, kv\+\_\+bbl, z\+\_\+i, h\+\_\+ml, dt, j, G, GV, US, CS, visc, forces, work\+\_\+on\+\_\+u, O\+BC, shelf)
\begin{DoxyCompactList}\small\item\em Calculate the \textquotesingle{}coupling coefficient\textquotesingle{} (a\+\_\+cpl) at the interfaces. If B\+O\+T\+T\+O\+M\+D\+R\+A\+G\+L\+AW is defined, the minimum of Hbbl and half the adjacent layer thicknesses are used to calculate a\+\_\+cpl near the bottom. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__vert__friction_ac14bc439f3b7ae3000fa9883c95ebfe2}{vertvisc\+\_\+limit\+\_\+vel}} (u, v, h, A\+Dp, C\+Dp, forces, visc, dt, G, GV, US, CS)
\begin{DoxyCompactList}\small\item\em Velocity components which exceed a threshold for physically reasonable values are truncated. Optionally, any column with excessive velocities may be sent to a diagnostic reporting subroutine. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__vert__friction_a4bd5c8772584e41890bff55ccd52507c}{vertvisc\+\_\+init}} (M\+IS, Time, G, GV, US, param\+\_\+file, diag, A\+Dp, dirs, ntrunc, CS)
\begin{DoxyCompactList}\small\item\em Initialize the vertical friction module. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__vert__friction_a62e586a80ed4bdd3fd27ab62ca4c054f}{updatecfltruncationvalue}} (Time, CS, activate)
\begin{DoxyCompactList}\small\item\em Update the C\+FL truncation value as a function of time. If called with the optional argument activate=.true., record the value of Time as the beginning of the ramp period. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacemom__vert__friction_a158d29cf5c1d79ca7c5f327706855972}{vertvisc\+\_\+end}} (CS)
\begin{DoxyCompactList}\small\item\em Clean up and deallocate the vertical friction module. \end{DoxyCompactList}\end{DoxyCompactItemize}


\subsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespacemom__vert__friction_aa9e6f1f0d75a54d85b0d0cdad874b41f}\label{namespacemom__vert__friction_aa9e6f1f0d75a54d85b0d0cdad874b41f}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!find\+\_\+coupling\+\_\+coef@{find\+\_\+coupling\+\_\+coef}}
\index{find\+\_\+coupling\+\_\+coef@{find\+\_\+coupling\+\_\+coef}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{find\+\_\+coupling\+\_\+coef()}{find\_coupling\_coef()}}
{\footnotesize\ttfamily subroutine mom\+\_\+vert\+\_\+friction\+::find\+\_\+coupling\+\_\+coef (\begin{DoxyParamCaption}\item[{real, dimension(szib\+\_\+(g),szk\+\_\+(gv)+1), intent(out)}]{a\+\_\+cpl,  }\item[{real, dimension(szib\+\_\+(g),szk\+\_\+(gv)), intent(in)}]{hvel,  }\item[{logical, dimension(szib\+\_\+(g)), intent(in)}]{do\+\_\+i,  }\item[{real, dimension(szib\+\_\+(g),szk\+\_\+(gv)), intent(in)}]{h\+\_\+harm,  }\item[{real, dimension(szib\+\_\+(g)), intent(in)}]{bbl\+\_\+thick,  }\item[{real, dimension(szib\+\_\+(g)), intent(in)}]{kv\+\_\+bbl,  }\item[{real, dimension(szib\+\_\+(g),szk\+\_\+(gv)+1), intent(in)}]{z\+\_\+i,  }\item[{real, dimension(szib\+\_\+(g)), intent(out)}]{h\+\_\+ml,  }\item[{real, intent(in)}]{dt,  }\item[{integer, intent(in)}]{j,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS,  }\item[{type(vertvisc\+\_\+type), intent(in)}]{visc,  }\item[{type(mech\+\_\+forcing), intent(in)}]{forces,  }\item[{logical, intent(in)}]{work\+\_\+on\+\_\+u,  }\item[{type(ocean\+\_\+obc\+\_\+type), pointer}]{O\+BC,  }\item[{logical, intent(in), optional}]{shelf }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Calculate the \textquotesingle{}coupling coefficient\textquotesingle{} (a\+\_\+cpl) at the interfaces. If B\+O\+T\+T\+O\+M\+D\+R\+A\+G\+L\+AW is defined, the minimum of Hbbl and half the adjacent layer thicknesses are used to calculate a\+\_\+cpl near the bottom. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt out}  & {\em a\+\_\+cpl} & Coupling coefficient across interfaces \mbox{[}Z T-\/1 $\sim$$>$ m s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em hvel} & Thickness at velocity points \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em do\+\_\+i} & If true, determine coupling coefficient for a column\\
\hline
\mbox{\tt in}  & {\em h\+\_\+harm} & Harmonic mean of thicknesses around a velocity\\
\hline
\mbox{\tt in}  & {\em bbl\+\_\+thick} & Bottom boundary layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em kv\+\_\+bbl} & Bottom boundary layer viscosity \mbox{[}Z2 T-\/1 $\sim$$>$ m2 s-\/1\mbox{]}.\\
\hline
\mbox{\tt in}  & {\em z\+\_\+i} & Estimate of interface heights above the bottom,\\
\hline
\mbox{\tt out}  & {\em h\+\_\+ml} & Mixed layer depth \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em j} & j-\/index to find coupling coefficient for\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
 & {\em cs} & Vertical viscosity control structure\\
\hline
\mbox{\tt in}  & {\em visc} & Structure containing viscosities and bottom drag\\
\hline
\mbox{\tt in}  & {\em forces} & A structure with the driving mechanical forces\\
\hline
\mbox{\tt in}  & {\em work\+\_\+on\+\_\+u} & If true, u-\/points are being calculated, otherwise they are v-\/points\\
\hline
 & {\em obc} & Open boundary condition structure\\
\hline
\mbox{\tt in}  & {\em shelf} & If present and true, use a surface boundary condition appropriate for an ice shelf. \\
\hline
\end{DoxyParams}


Definition at line 1092 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
1092   \textcolor{keywordtype}{type}(ocean\_grid\_type),     \textcolor{keywordtype}{intent(in)}  :: G\textcolor{comment}{  !< Ocean grid structure}
1093   \textcolor{keywordtype}{type}(verticalGrid\_type),   \textcolor{keywordtype}{intent(in)}  :: GV\textcolor{comment}{ !< Ocean vertical grid structure}
1094   \textcolor{keywordtype}{type}(unit\_scale\_type),     \textcolor{keywordtype}{intent(in)}  :: US\textcolor{comment}{ !< A dimensional unit scaling type}
1095   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(GV)+1)}, &
1096                              \textcolor{keywordtype}{intent(out)} :: a\_cpl\textcolor{comment}{ !< Coupling coefficient across interfaces [Z T-1 ~> m
       s-1].}
1097   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(GV))}, &
1098                              \textcolor{keywordtype}{intent(in)}  :: hvel\textcolor{comment}{ !< Thickness at velocity points [H ~> m or kg m-2]}
1099   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{dimension(SZIB\_(G))}, &
1100                              \textcolor{keywordtype}{intent(in)}  :: do\_i\textcolor{comment}{ !< If true, determine coupling coefficient for a column}
1101   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(GV))}, &
1102                              \textcolor{keywordtype}{intent(in)}  :: h\_harm\textcolor{comment}{ !< Harmonic mean of thicknesses around a velocity}
1103 \textcolor{comment}{                                                   !! grid point [H ~> m or kg m-2]}
1104   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))}, \textcolor{keywordtype}{intent(in)}  :: bbl\_thick\textcolor{comment}{ !< Bottom boundary layer thickness [H ~> m or kg m-2]}
1105   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))}, \textcolor{keywordtype}{intent(in)}  :: kv\_bbl\textcolor{comment}{ !< Bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].}
1106   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(GV)+1)}, &
1107                              \textcolor{keywordtype}{intent(in)}  :: z\_i\textcolor{comment}{  !< Estimate of interface heights above the bottom,}
1108 \textcolor{comment}{                                                 !! normalized by the bottom boundary layer thickness}
1109   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))}, \textcolor{keywordtype}{intent(out)} :: h\_ml\textcolor{comment}{ !< Mixed layer depth [H ~> m or kg m-2]}
1110   \textcolor{keywordtype}{integer},                   \textcolor{keywordtype}{intent(in)}  :: j\textcolor{comment}{    !< j-index to find coupling coefficient for}
1111   \textcolor{keywordtype}{real},                      \textcolor{keywordtype}{intent(in)}  :: dt\textcolor{comment}{   !< Time increment [T ~> s]}
1112   \textcolor{keywordtype}{type}(vertvisc\_CS),         \textcolor{keywordtype}{pointer}     :: CS\textcolor{comment}{   !< Vertical viscosity control structure}
1113   \textcolor{keywordtype}{type}(vertvisc\_type),       \textcolor{keywordtype}{intent(in)}  :: visc\textcolor{comment}{ !< Structure containing viscosities and bottom drag}
1114   \textcolor{keywordtype}{type}(mech\_forcing),        \textcolor{keywordtype}{intent(in)}  :: forces\textcolor{comment}{ !< A structure with the driving mechanical forces}
1115   \textcolor{keywordtype}{logical},                   \textcolor{keywordtype}{intent(in)}  :: work\_on\_u\textcolor{comment}{ !< If true, u-points are being calculated,}
1116 \textcolor{comment}{                                                  !! otherwise they are v-points}
1117   \textcolor{keywordtype}{type}(ocean\_OBC\_type),      \textcolor{keywordtype}{pointer}     :: OBC\textcolor{comment}{   !< Open boundary condition structure}
1118   \textcolor{keywordtype}{logical},         \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: shelf\textcolor{comment}{ !< If present and true, use a surface boundary}
1119 \textcolor{comment}{                                                  !! condition appropriate for an ice shelf.}
1120 
1121   \textcolor{comment}{! Local variables}
1122 
1123   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
1124     u\_star, &   \textcolor{comment}{! ustar at a velocity point [Z T-1 ~> m s-1].}
1125     absf, &     \textcolor{comment}{! The average of the neighboring absolute values of f [T-1 ~> s-1].}
1126 \textcolor{comment}{!      h\_ml, &  ! The mixed layer depth [H ~> m or kg m-2].}
1127     nk\_visc, &  \textcolor{comment}{! The (real) interface index of the base of mixed layer.}
1128     z\_t, &      \textcolor{comment}{! The distance from the top, sometimes normalized}
1129                 \textcolor{comment}{! by Hmix, [H ~> m or kg m-2] or [nondim].}
1130     kv\_tbl, &   \textcolor{comment}{! The viscosity in a top boundary layer under ice [Z2 T-1 ~> m2 s-1].}
1131     tbl\_thick
1132   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(GV)+1)} :: &
1133     Kv\_tot, &   \textcolor{comment}{! The total viscosity at an interface [Z2 T-1 ~> m2 s-1].}
1134     Kv\_add      \textcolor{comment}{! A viscosity to add [Z2 T-1 ~> m2 s-1].}
1135   \textcolor{keywordtype}{real} :: h\_shear \textcolor{comment}{! The distance over which shears occur [H ~> m or kg m-2].}
1136   \textcolor{keywordtype}{real} :: r       \textcolor{comment}{! A thickness to compare with Hbbl [H ~> m or kg m-2].}
1137   \textcolor{keywordtype}{real} :: visc\_ml \textcolor{comment}{! The mixed layer viscosity [Z2 T-1 ~> m2 s-1].}
1138   \textcolor{keywordtype}{real} :: I\_Hmix  \textcolor{comment}{! The inverse of the mixed layer thickness [H-1 ~> m-1 or m2 kg-1].}
1139   \textcolor{keywordtype}{real} :: a\_ml    \textcolor{comment}{! The layer coupling coefficient across an interface in}
1140                   \textcolor{comment}{! the mixed layer [Z T-1 ~> m s-1].}
1141   \textcolor{keywordtype}{real} :: I\_amax  \textcolor{comment}{! The inverse of the maximum coupling coefficient [T Z-1 ~> s m-1].}
1142   \textcolor{keywordtype}{real} :: temp1   \textcolor{comment}{! A temporary variable [H Z ~> m2 or kg m-1]}
1143   \textcolor{keywordtype}{real} :: h\_neglect   \textcolor{comment}{! A thickness that is so small it is usually lost}
1144                       \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
1145   \textcolor{keywordtype}{real} :: z2      \textcolor{comment}{! A copy of z\_i [nondim]}
1146   \textcolor{keywordtype}{real} :: botfn   \textcolor{comment}{! A function that is 1 at the bottom and small far from it [nondim]}
1147   \textcolor{keywordtype}{real} :: topfn   \textcolor{comment}{! A function that is 1 at the top and small far from it [nondim]}
1148   \textcolor{keywordtype}{real} :: kv\_top  \textcolor{comment}{! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1]}
1149   \textcolor{keywordtype}{logical} :: do\_shelf, do\_OBCs
1150   \textcolor{keywordtype}{integer} :: i, k, is, ie, max\_nk
1151   \textcolor{keywordtype}{integer} :: nz
1152 
1153   a\_cpl(:,:) = 0.0
1154   kv\_tot(:,:) = 0.0
1155 
1156   \textcolor{keywordflow}{if} (work\_on\_u) \textcolor{keywordflow}{then} ; is = g%IscB ; ie = g%IecB
1157   \textcolor{keywordflow}{else} ; is = g%isc ; ie = g%iec ;\textcolor{keywordflow}{ endif}
1158   nz = g%ke
1159   h\_neglect = gv%H\_subroundoff
1160 
1161   \textcolor{keywordflow}{if} (cs%answers\_2018) \textcolor{keywordflow}{then}
1162     \textcolor{comment}{!   The maximum coupling coefficent was originally introduced to avoid}
1163     \textcolor{comment}{! truncation error problems in the tridiagonal solver. Effectively, the 1e-10}
1164     \textcolor{comment}{! sets the maximum coupling coefficient increment to 1e10 m per timestep.}
1165     i\_amax = (1.0e-10*us%Z\_to\_m) * dt
1166   \textcolor{keywordflow}{else}
1167     i\_amax = 0.0
1168 \textcolor{keywordflow}{  endif}
1169 
1170   do\_shelf = .false. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(shelf)) do\_shelf = shelf
1171   do\_obcs = .false.
1172   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; do\_obcs = (obc%number\_of\_segments > 0) ;\textcolor{keywordflow}{ endif}
1173   h\_ml(:) = 0.0
1174 
1175 \textcolor{comment}{!    The following loop calculates the vertical average velocity and}
1176 \textcolor{comment}{!  surface mixed layer contributions to the vertical viscosity.}
1177   \textcolor{keywordflow}{do} i=is,ie ; kv\_tot(i,1) = 0.0 ;\textcolor{keywordflow}{ enddo}
1178   \textcolor{keywordflow}{if} ((gv%nkml>0) .or. do\_shelf) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie
1179     \textcolor{keywordflow}{if} (do\_i(i)) kv\_tot(i,k) = cs%Kv
1180 \textcolor{keywordflow}{  enddo} ;\textcolor{keywordflow}{ enddo} ; \textcolor{keywordflow}{else}
1181     i\_hmix = 1.0 / (cs%Hmix + h\_neglect)
1182     \textcolor{keywordflow}{do} i=is,ie ; z\_t(i) = h\_neglect*i\_hmix ;\textcolor{keywordflow}{ enddo}
1183     \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1184       z\_t(i) = z\_t(i) + h\_harm(i,k-1)*i\_hmix
1185       kv\_tot(i,k) = cs%Kv + cs%Kvml / ((z\_t(i)*z\_t(i)) *  &
1186                (1.0 + 0.09*z\_t(i)*z\_t(i)*z\_t(i)*z\_t(i)*z\_t(i)*z\_t(i)))
1187 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1188 \textcolor{keywordflow}{  endif}
1189 
1190   \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1191     \textcolor{keywordflow}{if} (cs%bottomdraglaw) \textcolor{keywordflow}{then}
1192       r = hvel(i,nz)*0.5
1193       \textcolor{keywordflow}{if} (r < bbl\_thick(i)) \textcolor{keywordflow}{then}
1194         a\_cpl(i,nz+1) = kv\_bbl(i) / (i\_amax*kv\_bbl(i) + (r+h\_neglect)*gv%H\_to\_Z)
1195       \textcolor{keywordflow}{else}
1196         a\_cpl(i,nz+1) = kv\_bbl(i) / (i\_amax*kv\_bbl(i) + (bbl\_thick(i)+h\_neglect)*gv%H\_to\_Z)
1197 \textcolor{keywordflow}{      endif}
1198     \textcolor{keywordflow}{else}
1199       a\_cpl(i,nz+1) = cs%Kvbbl / ((0.5*hvel(i,nz)+h\_neglect)*gv%H\_to\_Z + i\_amax*cs%Kvbbl)
1200 \textcolor{keywordflow}{    endif}
1201 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ enddo}
1202 
1203   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_shear)) \textcolor{keywordflow}{then}
1204     \textcolor{comment}{! The factor of 2 that used to be required in the viscosities is no longer needed.}
1205     \textcolor{keywordflow}{if} (work\_on\_u) \textcolor{keywordflow}{then}
1206       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1207         kv\_add(i,k) = 0.5*(visc%Kv\_shear(i,j,k) + visc%Kv\_shear(i+1,j,k))
1208 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1209       \textcolor{keywordflow}{if} (do\_obcs) \textcolor{keywordflow}{then}
1210         \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_u(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
1211           \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_e) \textcolor{keywordflow}{then}
1212             \textcolor{keywordflow}{do} k=2,nz ; kv\_add(i,k) = visc%Kv\_shear(i,j,k) ;\textcolor{keywordflow}{ enddo}
1213           \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_w) \textcolor{keywordflow}{then}
1214             \textcolor{keywordflow}{do} k=2,nz ; kv\_add(i,k) = visc%Kv\_shear(i+1,j,k) ;\textcolor{keywordflow}{ enddo}
1215 \textcolor{keywordflow}{          endif}
1216 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo}
1217 \textcolor{keywordflow}{      endif}
1218       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1219         kv\_tot(i,k) = kv\_tot(i,k) + kv\_add(i,k)
1220 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1221     \textcolor{keywordflow}{else}
1222       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1223         kv\_add(i,k) = 0.5*(visc%Kv\_shear(i,j,k) + visc%Kv\_shear(i,j+1,k))
1224 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1225       \textcolor{keywordflow}{if} (do\_obcs) \textcolor{keywordflow}{then}
1226         \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_v(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
1227           \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_n) \textcolor{keywordflow}{then}
1228             \textcolor{keywordflow}{do} k=2,nz ; kv\_add(i,k) = visc%Kv\_shear(i,j,k) ;\textcolor{keywordflow}{ enddo}
1229           \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_s) \textcolor{keywordflow}{then}
1230             \textcolor{keywordflow}{do} k=2,nz ; kv\_add(i,k) = visc%Kv\_shear(i,j+1,k) ;\textcolor{keywordflow}{ enddo}
1231 \textcolor{keywordflow}{          endif}
1232 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo}
1233 \textcolor{keywordflow}{      endif}
1234       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1235         kv\_tot(i,k) = kv\_tot(i,k) + kv\_add(i,k)
1236 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1237 \textcolor{keywordflow}{    endif}
1238 \textcolor{keywordflow}{  endif}
1239 
1240   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_shear\_Bu)) \textcolor{keywordflow}{then}
1241     \textcolor{keywordflow}{if} (work\_on\_u) \textcolor{keywordflow}{then}
1242       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{If} (do\_i(i)) \textcolor{keywordflow}{then}
1243         kv\_tot(i,k) = kv\_tot(i,k) + (0.5)*(visc%Kv\_shear\_Bu(i,j-1,k) + visc%Kv\_shear\_Bu(i,j,k))
1244 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1245     \textcolor{keywordflow}{else}
1246       \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1247         kv\_tot(i,k) = kv\_tot(i,k) + (0.5)*(visc%Kv\_shear\_Bu(i-1,j,k) + visc%Kv\_shear\_Bu(i,j,k))
1248 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1249 \textcolor{keywordflow}{    endif}
1250 \textcolor{keywordflow}{  endif}
1251 
1252   \textcolor{keywordflow}{do} k=nz,2,-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1253     \textcolor{comment}{!    botfn determines when a point is within the influence of the bottom}
1254     \textcolor{comment}{!  boundary layer, going from 1 at the bottom to 0 in the interior.}
1255     z2 = z\_i(i,k)
1256     botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
1257 
1258     \textcolor{keywordflow}{if} (cs%bottomdraglaw) \textcolor{keywordflow}{then}
1259       kv\_tot(i,k) = kv\_tot(i,k) + (kv\_bbl(i) - cs%Kv)*botfn
1260       r = 0.5*(hvel(i,k) + hvel(i,k-1))
1261       \textcolor{keywordflow}{if} (r > bbl\_thick(i)) \textcolor{keywordflow}{then}
1262         h\_shear = ((1.0 - botfn) * r + botfn*bbl\_thick(i)) + h\_neglect
1263       \textcolor{keywordflow}{else}
1264         h\_shear = r + h\_neglect
1265 \textcolor{keywordflow}{      endif}
1266     \textcolor{keywordflow}{else}
1267       kv\_tot(i,k) = kv\_tot(i,k) + (cs%Kvbbl-cs%Kv)*botfn
1268       h\_shear = 0.5*(hvel(i,k) + hvel(i,k-1) + h\_neglect)
1269 \textcolor{keywordflow}{    endif}
1270 
1271     \textcolor{comment}{! Calculate the coupling coefficients from the viscosities.}
1272     a\_cpl(i,k) = kv\_tot(i,k) / (h\_shear*gv%H\_to\_Z + i\_amax*kv\_tot(i,k))
1273 \textcolor{keywordflow}{  endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i & k loops}
1274 
1275   \textcolor{keywordflow}{if} (do\_shelf) \textcolor{keywordflow}{then}
1276     \textcolor{comment}{! Set the coefficients to include the no-slip surface stress.}
1277     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1278       \textcolor{keywordflow}{if} (work\_on\_u) \textcolor{keywordflow}{then}
1279         kv\_tbl(i) = visc%Kv\_tbl\_shelf\_u(i,j)
1280         tbl\_thick(i) = visc%tbl\_thick\_shelf\_u(i,j) * gv%Z\_to\_H + h\_neglect
1281       \textcolor{keywordflow}{else}
1282         kv\_tbl(i) = visc%Kv\_tbl\_shelf\_v(i,j)
1283         tbl\_thick(i) = visc%tbl\_thick\_shelf\_v(i,j) * gv%Z\_to\_H + h\_neglect
1284 \textcolor{keywordflow}{      endif}
1285       z\_t(i) = 0.0
1286 
1287       \textcolor{comment}{! If a\_cpl(i,1) were not already 0, it would be added here.}
1288       \textcolor{keywordflow}{if} (0.5*hvel(i,1) > tbl\_thick(i)) \textcolor{keywordflow}{then}
1289         a\_cpl(i,1) = kv\_tbl(i) / (tbl\_thick(i)*gv%H\_to\_Z + i\_amax*kv\_tbl(i))
1290       \textcolor{keywordflow}{else}
1291         a\_cpl(i,1) = kv\_tbl(i) / ((0.5*hvel(i,1)+h\_neglect)*gv%H\_to\_Z + i\_amax*kv\_tbl(i))
1292 \textcolor{keywordflow}{      endif}
1293 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo}
1294 
1295     \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ;  \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1296       z\_t(i) = z\_t(i) + hvel(i,k-1) / tbl\_thick(i)
1297       topfn = 1.0 / (1.0 + 0.09 * z\_t(i)**6)
1298 
1299       r = 0.5*(hvel(i,k)+hvel(i,k-1))
1300       \textcolor{keywordflow}{if} (r > tbl\_thick(i)) \textcolor{keywordflow}{then}
1301         h\_shear = ((1.0 - topfn) * r + topfn*tbl\_thick(i)) + h\_neglect
1302       \textcolor{keywordflow}{else}
1303         h\_shear = r + h\_neglect
1304 \textcolor{keywordflow}{      endif}
1305 
1306       kv\_top = topfn * kv\_tbl(i)
1307       a\_cpl(i,k) = a\_cpl(i,k) + kv\_top / (h\_shear*gv%H\_to\_Z + i\_amax*kv\_top)
1308 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1309   \textcolor{keywordflow}{elseif} (cs%dynamic\_viscous\_ML .or. (gv%nkml>0)) \textcolor{keywordflow}{then}
1310     max\_nk = 0
1311     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1312       \textcolor{keywordflow}{if} (gv%nkml>0) nk\_visc(i) = \textcolor{keywordtype}{real}(GV%nkml+1)
1313       \textcolor{keywordflow}{if} (work\_on\_u) \textcolor{keywordflow}{then}
1314         u\_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
1315         absf(i) = 0.5*(abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))
1316         \textcolor{keywordflow}{if} (cs%dynamic\_viscous\_ML) nk\_visc(i) = visc%nkml\_visc\_u(i,j) + 1
1317       \textcolor{keywordflow}{else}
1318         u\_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
1319         absf(i) = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1320         \textcolor{keywordflow}{if} (cs%dynamic\_viscous\_ML) nk\_visc(i) = visc%nkml\_visc\_v(i,j) + 1
1321 \textcolor{keywordflow}{      endif}
1322       h\_ml(i) = h\_neglect ; z\_t(i) = 0.0
1323       max\_nk = max(max\_nk,ceiling(nk\_visc(i) - 1.0))
1324 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo}
1325 
1326     \textcolor{keywordflow}{if} (do\_obcs) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (work\_on\_u) \textcolor{keywordflow}{then}
1327       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_u(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
1328         \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_e) &
1329           u\_star(i) = forces%ustar(i,j)
1330         \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_w) &
1331           u\_star(i) = forces%ustar(i+1,j)
1332 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
1333     \textcolor{keywordflow}{else}
1334       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_v(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
1335         \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_n) &
1336           u\_star(i) = forces%ustar(i,j)
1337         \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_s) &
1338           u\_star(i) = forces%ustar(i,j+1)
1339 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
1340 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
1341 
1342     \textcolor{keywordflow}{do} k=1,max\_nk ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
1343       \textcolor{keywordflow}{if} (k+1 <= nk\_visc(i)) \textcolor{keywordflow}{then} \textcolor{comment}{! This layer is all in the ML.}
1344         h\_ml(i) = h\_ml(i) + hvel(i,k)
1345       \textcolor{keywordflow}{elseif} (k < nk\_visc(i)) \textcolor{keywordflow}{then} \textcolor{comment}{! Part of this layer is in the ML.}
1346         h\_ml(i) = h\_ml(i) + (nk\_visc(i) - k) * hvel(i,k)
1347 \textcolor{keywordflow}{      endif}
1348 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1349 
1350     \textcolor{keywordflow}{do} k=2,max\_nk ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (k < nk\_visc(i)) \textcolor{keywordflow}{then}
1351       \textcolor{comment}{! Set the viscosity at the interfaces.}
1352       z\_t(i) = z\_t(i) + hvel(i,k-1)
1353       temp1 = (z\_t(i)*h\_ml(i) - z\_t(i)*z\_t(i))*gv%H\_to\_Z
1354       \textcolor{comment}{!   This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)}
1355       \textcolor{comment}{! and be further limited by rotation to give the natural Ekman length.}
1356       visc\_ml = u\_star(i) * 0.41 * (temp1*u\_star(i)) / (absf(i)*temp1 + (h\_ml(i)+h\_neglect)*u\_star(i))
1357       a\_ml = visc\_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h\_neglect) * gv%H\_to\_Z + 0.5*i\_amax*visc\_ml)
1358       \textcolor{comment}{! Choose the largest estimate of a.}
1359       \textcolor{keywordflow}{if} (a\_ml > a\_cpl(i,k)) a\_cpl(i,k) = a\_ml
1360 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1361 \textcolor{keywordflow}{  endif}
1362 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__vert__friction_a62e586a80ed4bdd3fd27ab62ca4c054f}\label{namespacemom__vert__friction_a62e586a80ed4bdd3fd27ab62ca4c054f}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!updatecfltruncationvalue@{updatecfltruncationvalue}}
\index{updatecfltruncationvalue@{updatecfltruncationvalue}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{updatecfltruncationvalue()}{updatecfltruncationvalue()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+vert\+\_\+friction\+::updatecfltruncationvalue (\begin{DoxyParamCaption}\item[{type(time\+\_\+type), intent(in), target}]{Time,  }\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS,  }\item[{logical, intent(in), optional}]{activate }\end{DoxyParamCaption})}



Update the C\+FL truncation value as a function of time. If called with the optional argument activate=.true., record the value of Time as the beginning of the ramp period. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em time} & Current model time\\
\hline
 & {\em cs} & Vertical viscosity control structure\\
\hline
\mbox{\tt in}  & {\em activate} & Specifiy whether to record the value of Time as the beginning of the ramp period \\
\hline
\end{DoxyParams}


Definition at line 1852 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
1852   \textcolor{keywordtype}{type}(time\_type), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)}    :: Time\textcolor{comment}{     !< Current model time}
1853   \textcolor{keywordtype}{type}(vertvisc\_CS),       \textcolor{keywordtype}{pointer}       :: CS\textcolor{comment}{       !< Vertical viscosity control structure}
1854   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{optional},       \textcolor{keywordtype}{intent(in)}    :: activate\textcolor{comment}{ !< Specifiy whether to record the value of}
1855 \textcolor{comment}{                                                     !! Time as the beginning of the ramp period}
1856 
1857   \textcolor{comment}{! Local variables}
1858   \textcolor{keywordtype}{real} :: deltaTime, wghtA
1859   \textcolor{keywordtype}{character(len=12)} :: msg
1860 
1861   \textcolor{keywordflow}{if} (cs%truncRampTime==0.) \textcolor{keywordflow}{return} \textcolor{comment}{! This indicates to ramping is turned off}
1862 
1863   \textcolor{comment}{! We use the optional argument to indicate this Time should be recorded as the}
1864   \textcolor{comment}{! beginning of the ramp-up period.}
1865   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(activate)) \textcolor{keywordflow}{then}
1866     \textcolor{keywordflow}{if} (activate) \textcolor{keywordflow}{then}
1867       cs%rampStartTime = time \textcolor{comment}{! Record the current time}
1868       cs%CFLrampingIsActivated = .true.
1869 \textcolor{keywordflow}{    endif}
1870 \textcolor{keywordflow}{  endif}
1871   \textcolor{keywordflow}{if} (.not.cs%CFLrampingIsActivated) \textcolor{keywordflow}{return}
1872   deltatime = max( 0., time\_type\_to\_real( time - cs%rampStartTime ) )
1873   \textcolor{keywordflow}{if} (deltatime >= cs%truncRampTime) \textcolor{keywordflow}{then}
1874     cs%CFL\_trunc = cs%CFL\_truncE
1875     cs%truncRampTime = 0. \textcolor{comment}{! This turns off ramping after this call}
1876   \textcolor{keywordflow}{else}
1877     wghta = min( 1., deltatime / cs%truncRampTime ) \textcolor{comment}{! Linear profile in time}
1878     \textcolor{comment}{!wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time}
1879     \textcolor{comment}{!wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile}
1880     wghta = 1. - ( (1. - wghta)**2 ) \textcolor{comment}{! Convert linear profiel to nverted parabolic profile}
1881     cs%CFL\_trunc = cs%CFL\_truncS + wghta * ( cs%CFL\_truncE - cs%CFL\_truncS )
1882 \textcolor{keywordflow}{  endif}
1883   \textcolor{keyword}{write}(msg(1:12),\textcolor{stringliteral}{'(es12.3)'}) cs%CFL\_trunc
1884   \textcolor{keyword}{call }mom\_error(note, \textcolor{stringliteral}{"MOM\_vert\_friction: updateCFLtruncationValue set CFL"}// &
1885                        \textcolor{stringliteral}{" limit to "}//trim(msg))
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__vert__friction_a8f1a390fa24fbe985068fed9ac26873c}\label{namespacemom__vert__friction_a8f1a390fa24fbe985068fed9ac26873c}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!vertvisc@{vertvisc}}
\index{vertvisc@{vertvisc}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{vertvisc()}{vertvisc()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+vert\+\_\+friction\+::vertvisc (\begin{DoxyParamCaption}\item[{real, dimension( g \%isdb\+: g \%iedb, g \%jsd\+: g \%jed, gv \%ke), intent(inout)}]{u,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsdb\+: g \%jedb, gv \%ke), intent(inout)}]{v,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsd\+: g \%jed, gv \%ke), intent(in)}]{h,  }\item[{type(mech\+\_\+forcing), intent(in)}]{forces,  }\item[{type(vertvisc\+\_\+type), intent(inout)}]{visc,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+obc\+\_\+type), pointer}]{O\+BC,  }\item[{type(accel\+\_\+diag\+\_\+ptrs), intent(inout)}]{A\+Dp,  }\item[{type(cont\+\_\+diag\+\_\+ptrs), intent(inout)}]{C\+Dp,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS,  }\item[{real, dimension( g \%isdb\+: g \%iedb, g \%jsd\+: g \%jed), intent(out), optional}]{taux\+\_\+bot,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsdb\+: g \%jedb), intent(out), optional}]{tauy\+\_\+bot,  }\item[{type(wave\+\_\+parameters\+\_\+cs), optional, pointer}]{Waves }\end{DoxyParamCaption})}



Perform a fully implicit vertical diffusion of momentum. Stress top and bottom boundary conditions are used. 

This is solving the tridiagonal system \[ \left(h_k + a_{k + 1/2} + a_{k - 1/2} + r_k\right) u_k^{n+1} = h_k u_k^n + a_{k + 1/2} u_{k+1}^{n+1} + a_{k - 1/2} u_{k-1}^{n+1} \] where $a_{k + 1/2} = \Delta t \nu_{k + 1/2} / h_{k + 1/2}$ is the {\itshape interfacial coupling thickness per time step}, encompassing background viscosity as well as contributions from enhanced mixed and bottom layer viscosities. \$r\+\_\+k\$ is a Rayleight drag term due to channel drag. There is an additional stress term on the right-\/hand side if D\+I\+R\+E\+C\+T\+\_\+\+S\+T\+R\+E\+SS is true, applied to the surface layer.


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in,out}  & {\em u} & Zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em v} & Meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em forces} & A structure with the driving mechanical forces\\
\hline
\mbox{\tt in,out}  & {\em visc} & Viscosities and bottom drag\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
 & {\em obc} & Open boundary condition structure\\
\hline
\mbox{\tt in,out}  & {\em adp} & Accelerations in the momentum equations for diagnostics\\
\hline
\mbox{\tt in,out}  & {\em cdp} & Continuity equation terms\\
\hline
 & {\em cs} & Vertical viscosity control structure\\
\hline
\mbox{\tt out}  & {\em taux\+\_\+bot} & Zonal bottom stress from ocean to\\
\hline
\mbox{\tt out}  & {\em tauy\+\_\+bot} & Meridional bottom stress from ocean to\\
\hline
 & {\em waves} & Container for wave/\+Stokes information \\
\hline
\end{DoxyParams}


Definition at line 159 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
159   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: G\textcolor{comment}{      !< Ocean grid structure}
160   \textcolor{keywordtype}{type}(verticalGrid\_type), \textcolor{keywordtype}{intent(in)}    :: GV\textcolor{comment}{     !< Ocean vertical grid structure}
161   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: US\textcolor{comment}{     !< A dimensional unit scaling type}
162   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(GV))}, &
163                            \textcolor{keywordtype}{intent(inout)} :: u\textcolor{comment}{      !< Zonal velocity [L T-1 ~> m s-1]}
164   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(GV))}, &
165                            \textcolor{keywordtype}{intent(inout)} :: v\textcolor{comment}{      !< Meridional velocity [L T-1 ~> m s-1]}
166   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))}, &
167                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{      !< Layer thickness [H ~> m or kg m-2]}
168   \textcolor{keywordtype}{type}(mech\_forcing),    \textcolor{keywordtype}{intent(in)}      :: forces\textcolor{comment}{ !< A structure with the driving mechanical forces}
169   \textcolor{keywordtype}{type}(vertvisc\_type),   \textcolor{keywordtype}{intent(inout)}   :: visc\textcolor{comment}{   !< Viscosities and bottom drag}
170   \textcolor{keywordtype}{real},                  \textcolor{keywordtype}{intent(in)}      :: dt\textcolor{comment}{     !< Time increment [T ~> s]}
171   \textcolor{keywordtype}{type}(ocean\_OBC\_type),  \textcolor{keywordtype}{pointer}         :: OBC\textcolor{comment}{    !< Open boundary condition structure}
172   \textcolor{keywordtype}{type}(accel\_diag\_ptrs), \textcolor{keywordtype}{intent(inout)}   :: ADp\textcolor{comment}{    !< Accelerations in the momentum}
173 \textcolor{comment}{                                                   !! equations for diagnostics}
174   \textcolor{keywordtype}{type}(cont\_diag\_ptrs),  \textcolor{keywordtype}{intent(inout)}   :: CDp\textcolor{comment}{    !< Continuity equation terms}
175   \textcolor{keywordtype}{type}(vertvisc\_CS),     \textcolor{keywordtype}{pointer}         :: CS\textcolor{comment}{     !< Vertical viscosity control structure}
176   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G))}, &
177                    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: taux\_bot\textcolor{comment}{ !< Zonal bottom stress from ocean to}
178 \textcolor{comment}{                                                     !! rock [R L Z T-2 ~> Pa]}
179   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G))}, &
180                    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(out)} :: tauy\_bot\textcolor{comment}{ !< Meridional bottom stress from ocean to}
181 \textcolor{comment}{                                                     !! rock [R L Z T-2 ~> Pa]}
182   \textcolor{keywordtype}{type}(wave\_parameters\_CS), &
183                    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{pointer}     :: Waves\textcolor{comment}{ !< Container for wave/Stokes information}
184 
185   \textcolor{comment}{! Fields from forces used in this subroutine:}
186   \textcolor{comment}{!   taux: Zonal wind stress [R L Z T-2 ~> Pa].}
187   \textcolor{comment}{!   tauy: Meridional wind stress [R L Z T-2 ~> Pa].}
188 
189   \textcolor{comment}{! Local variables}
190 
191   \textcolor{keywordtype}{real} :: b1(SZIB\_(G))          \textcolor{comment}{! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].}
192   \textcolor{keywordtype}{real} :: c1(SZIB\_(G),SZK\_(G))  \textcolor{comment}{! A variable used by the tridiagonal solver [nondim].}
193   \textcolor{keywordtype}{real} :: d1(SZIB\_(G))          \textcolor{comment}{! d1=1-c1 is used by the tridiagonal solver [nondim].}
194   \textcolor{keywordtype}{real} :: Ray(SZIB\_(G),SZK\_(G)) \textcolor{comment}{! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].}
195   \textcolor{keywordtype}{real} :: b\_denom\_1             \textcolor{comment}{! The first term in the denominator of b1 [H ~> m or kg m-2].}
196 
197   \textcolor{keywordtype}{real} :: Hmix             \textcolor{comment}{! The mixed layer thickness over which stress}
198                            \textcolor{comment}{! is applied with direct\_stress [H ~> m or kg m-2].}
199   \textcolor{keywordtype}{real} :: I\_Hmix           \textcolor{comment}{! The inverse of Hmix [H-1 ~> m-1 or m2 kg-1].}
200   \textcolor{keywordtype}{real} :: Idt              \textcolor{comment}{! The inverse of the time step [T-1 ~> s-1].}
201   \textcolor{keywordtype}{real} :: dt\_Rho0          \textcolor{comment}{! The time step divided by the mean density [T H Z-1 R-1 ~> s m3 kg-1 or s].}
202   \textcolor{keywordtype}{real} :: dt\_Z\_to\_H        \textcolor{comment}{! The time step times the conversion from Z to the}
203                            \textcolor{comment}{! units of thickness - [T H Z-1 ~> s or s kg m-3].}
204   \textcolor{keywordtype}{real} :: h\_neglect        \textcolor{comment}{! A thickness that is so small it is usually lost}
205                            \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
206 
207   \textcolor{keywordtype}{real} :: stress           \textcolor{comment}{!   The surface stress times the time step, divided}
208                            \textcolor{comment}{! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1].}
209   \textcolor{keywordtype}{real} :: zDS, hfr, h\_a    \textcolor{comment}{! Temporary variables used with direct\_stress.}
210   \textcolor{keywordtype}{real} :: surface\_stress(SZIB\_(G))\textcolor{comment}{! The same as stress, unless the wind stress}
211                            \textcolor{comment}{! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1].}
212 
213   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:,:)} :: hf\_du\_dt\_visc\_2d \textcolor{comment}{! Depth sum of hf\_du\_dt\_visc [L T-2 ~> m s-2]}
214   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:,:)} :: hf\_dv\_dt\_visc\_2d \textcolor{comment}{! Depth sum of hf\_dv\_dt\_visc [L T-2 ~> m s-2]}
215 
216   \textcolor{keywordtype}{logical} :: do\_i(SZIB\_(G))
217   \textcolor{keywordtype}{logical} :: DoStokesMixing
218 
219   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
220   is = g%isc ; ie = g%iec; js = g%jsc; je = g%jec
221   isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
222 
223   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal,\textcolor{stringliteral}{"MOM\_vert\_friction(visc): "}// &
224          \textcolor{stringliteral}{"Module must be initialized before it is used."})
225 
226   \textcolor{keywordflow}{if} (cs%direct\_stress) \textcolor{keywordflow}{then}
227     hmix = cs%Hmix\_stress
228     i\_hmix = 1.0 / hmix
229 \textcolor{keywordflow}{  endif}
230   dt\_rho0 = dt / gv%H\_to\_RZ
231   dt\_z\_to\_h = dt*gv%Z\_to\_H
232   h\_neglect = gv%H\_subroundoff
233   idt = 1.0 / dt
234 
235   \textcolor{comment}{!Check if Stokes mixing allowed if requested (present and associated)}
236   dostokesmixing=.false.
237   \textcolor{keywordflow}{if} (cs%StokesMixing) \textcolor{keywordflow}{then}
238     \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(waves)) dostokesmixing = \textcolor{keyword}{associated}(waves)
239     \textcolor{keywordflow}{if} (.not. dostokesmixing) &
240       \textcolor{keyword}{call }mom\_error(fatal,\textcolor{stringliteral}{"Stokes Mixing called without allocated"}//&
241                      \textcolor{stringliteral}{"Waves Control Structure"})
242 \textcolor{keywordflow}{  endif}
243 
244   \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; ray(i,k) = 0.0 ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
245 
246   \textcolor{comment}{!   Update the zonal velocity component using a modification of a standard}
247   \textcolor{comment}{! tridagonal solver.}
248 
249   \textcolor{comment}{!$OMP parallel do default(shared) firstprivate(Ray) &}
250   \textcolor{comment}{!$OMP                 private(do\_i,surface\_stress,zDS,stress,h\_a,hfr, &}
251   \textcolor{comment}{!$OMP                         b\_denom\_1,b1,d1,c1)}
252   \textcolor{keywordflow}{do} j=g%jsc,g%jec
253     \textcolor{keywordflow}{do} i=isq,ieq ; do\_i(i) = (g%mask2dCu(i,j) > 0) ;\textcolor{keywordflow}{ enddo}
254 
255     \textcolor{comment}{! When mixing down Eulerian current + Stokes drift add before calling solver}
256     \textcolor{keywordflow}{if} (dostokesmixing) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
257       \textcolor{keywordflow}{if} (do\_i(i)) u(i,j,k) = u(i,j,k) + us%m\_s\_to\_L\_T*waves%Us\_x(i,j,k)
258 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
259 
260     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(adp%du\_dt\_visc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
261       adp%du\_dt\_visc(i,j,k) = u(i,j,k)
262 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
263 
264 \textcolor{comment}{!   One option is to have the wind stress applied as a body force}
265 \textcolor{comment}{! over the topmost Hmix fluid.  If DIRECT\_STRESS is not defined,}
266 \textcolor{comment}{! the wind stress is applied as a stress boundary condition.}
267     \textcolor{keywordflow}{if} (cs%direct\_stress) \textcolor{keywordflow}{then}
268       \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
269         surface\_stress(i) = 0.0
270         zds = 0.0
271         stress = dt\_rho0 * forces%taux(i,j)
272         \textcolor{keywordflow}{do} k=1,nz
273           h\_a = 0.5 * (h(i,j,k) + h(i+1,j,k)) + h\_neglect
274           hfr = 1.0 ; \textcolor{keywordflow}{if} ((zds+h\_a) > hmix) hfr = (hmix - zds) / h\_a
275           u(i,j,k) = u(i,j,k) + i\_hmix * hfr * stress
276           zds = zds + h\_a ; \textcolor{keywordflow}{if} (zds >= hmix) \textcolor{keywordflow}{exit}
277 \textcolor{keywordflow}{        enddo}
278 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end of i loop}
279     \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{do} i=isq,ieq
280       surface\_stress(i) = dt\_rho0 * (g%mask2dCu(i,j)*forces%taux(i,j))
281 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ endif} \textcolor{comment}{! direct\_stress}
282 
283     \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
284       ray(i,k) = visc%Ray\_u(i,j,k)
285 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
286 
287     \textcolor{comment}{! perform forward elimination on the tridiagonal system}
288     \textcolor{comment}{!}
289     \textcolor{comment}{! denote the diagonal of the system as b\_k, the subdiagonal as a\_k}
290     \textcolor{comment}{! and the superdiagonal as c\_k. The right-hand side terms are d\_k.}
291     \textcolor{comment}{!}
292     \textcolor{comment}{! ignoring the rayleigh drag contribution,}
293     \textcolor{comment}{! we have a\_k = -dt\_Z\_to\_H * a\_u(k)}
294     \textcolor{comment}{!         b\_k = h\_u(k) + dt\_Z\_to\_H * (a\_u(k) + a\_u(k+1))}
295     \textcolor{comment}{!         c\_k = -dt\_Z\_to\_H * a\_u(k+1)}
296     \textcolor{comment}{!}
297     \textcolor{comment}{! for forward elimination, we want to:}
298     \textcolor{comment}{! calculate c'\_k = - c\_k                / (b\_k + a\_k c'\_(k-1))}
299     \textcolor{comment}{! and       d'\_k = (d\_k - a\_k d'\_(k-1)) / (b\_k + a\_k c'\_(k-1))}
300     \textcolor{comment}{! where c'\_1 = c\_1/b\_1 and d'\_1 = d\_1/b\_1}
301     \textcolor{comment}{! (see Thomas' tridiagonal matrix algorithm)}
302     \textcolor{comment}{!}
303     \textcolor{comment}{! b1 is the denominator term 1 / (b\_k + a\_k c'\_(k-1))}
304     \textcolor{comment}{! b\_denom\_1 is (b\_k + a\_k + c\_k) - a\_k(1 - c'\_(k-1))}
305     \textcolor{comment}{!            = (b\_k + c\_k + c'\_(k-1))}
306     \textcolor{comment}{! this is done so that d1 = b1 * b\_denom\_1 = 1 - c'\_(k-1)}
307     \textcolor{comment}{! c1(k) is -c'\_(k - 1)}
308     \textcolor{comment}{! and the right-hand-side is destructively updated to be d'\_k}
309     \textcolor{comment}{!}
310     \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
311       b\_denom\_1 = cs%h\_u(i,j,1) + dt\_z\_to\_h * (ray(i,1) + cs%a\_u(i,j,1))
312       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h*cs%a\_u(i,j,2))
313       d1(i) = b\_denom\_1 * b1(i)
314       u(i,j,1) = b1(i) * (cs%h\_u(i,j,1) * u(i,j,1) + surface\_stress(i))
315 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo}
316     \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
317       c1(i,k) = dt\_z\_to\_h * cs%a\_u(i,j,k) * b1(i)
318       b\_denom\_1 = cs%h\_u(i,j,k) + dt\_z\_to\_h * (ray(i,k) + cs%a\_u(i,j,k)*d1(i))
319       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h * cs%a\_u(i,j,k+1))
320       d1(i) = b\_denom\_1 * b1(i)
321       u(i,j,k) = (cs%h\_u(i,j,k) * u(i,j,k) + &
322                   dt\_z\_to\_h * cs%a\_u(i,j,k) * u(i,j,k-1)) * b1(i)
323 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
324 
325     \textcolor{comment}{! back substitute to solve for the new velocities}
326     \textcolor{comment}{! u\_k = d'\_k - c'\_k x\_(k+1)}
327     \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
328       u(i,j,k) = u(i,j,k) + c1(i,k+1) * u(i,j,k+1)
329 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i and k loops}
330 
331     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(adp%du\_dt\_visc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
332       adp%du\_dt\_visc(i,j,k) = (u(i,j,k) - adp%du\_dt\_visc(i,j,k))*idt
333 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
334 
335     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%taux\_shelf)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=isq,ieq
336       visc%taux\_shelf(i,j) = -gv%Rho0*cs%a1\_shelf\_u(i,j)*u(i,j,1) \textcolor{comment}{! - u\_shelf?}
337 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ endif}
338 
339     \textcolor{keywordflow}{if} (\textcolor{keyword}{PRESENT}(taux\_bot)) \textcolor{keywordflow}{then}
340       \textcolor{keywordflow}{do} i=isq,ieq
341         taux\_bot(i,j) = gv%Rho0 * (u(i,j,nz)*cs%a\_u(i,j,nz+1))
342 \textcolor{keywordflow}{      enddo}
343       \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
344         taux\_bot(i,j) = taux\_bot(i,j) + gv%Rho0 * (ray(i,k)*u(i,j,k))
345 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
346 \textcolor{keywordflow}{    endif}
347 
348     \textcolor{comment}{! When mixing down Eulerian current + Stokes drift subtract after calling solver}
349     \textcolor{keywordflow}{if} (dostokesmixing) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
350       \textcolor{keywordflow}{if} (do\_i(i)) u(i,j,k) = u(i,j,k) - us%m\_s\_to\_L\_T*waves%Us\_x(i,j,k)
351 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
352 
353 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end u-component j loop}
354 
355   \textcolor{comment}{! Now work on the meridional velocity component.}
356 
357   \textcolor{comment}{!$OMP parallel do default(shared) firstprivate(Ray) &}
358   \textcolor{comment}{!$OMP               private(do\_i,surface\_stress,zDS,stress,h\_a,hfr, &}
359   \textcolor{comment}{!$OMP                       b\_denom\_1,b1,d1,c1)}
360   \textcolor{keywordflow}{do} j=jsq,jeq
361     \textcolor{keywordflow}{do} i=is,ie ; do\_i(i) = (g%mask2dCv(i,j) > 0) ;\textcolor{keywordflow}{ enddo}
362 
363     \textcolor{comment}{! When mixing down Eulerian current + Stokes drift add before calling solver}
364     \textcolor{keywordflow}{if} (dostokesmixing) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
365       \textcolor{keywordflow}{if} (do\_i(i)) v(i,j,k) = v(i,j,k) + us%m\_s\_to\_L\_T*waves%Us\_y(i,j,k)
366 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
367 
368     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(adp%dv\_dt\_visc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
369       adp%dv\_dt\_visc(i,j,k) = v(i,j,k)
370 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
371 
372 \textcolor{comment}{!   One option is to have the wind stress applied as a body force}
373 \textcolor{comment}{! over the topmost Hmix fluid.  If DIRECT\_STRESS is not defined,}
374 \textcolor{comment}{! the wind stress is applied as a stress boundary condition.}
375     \textcolor{keywordflow}{if} (cs%direct\_stress) \textcolor{keywordflow}{then}
376       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
377         surface\_stress(i) = 0.0
378         zds = 0.0
379         stress = dt\_rho0 * forces%tauy(i,j)
380         \textcolor{keywordflow}{do} k=1,nz
381           h\_a = 0.5 * (h(i,j,k) + h(i,j+1,k)) + h\_neglect
382           hfr = 1.0 ; \textcolor{keywordflow}{if} ((zds+h\_a) > hmix) hfr = (hmix - zds) / h\_a
383           v(i,j,k) = v(i,j,k) + i\_hmix * hfr * stress
384           zds = zds + h\_a ; \textcolor{keywordflow}{if} (zds >= hmix) \textcolor{keywordflow}{exit}
385 \textcolor{keywordflow}{        enddo}
386 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end of i loop}
387     \textcolor{keywordflow}{else} ; \textcolor{keywordflow}{do} i=is,ie
388       surface\_stress(i) = dt\_rho0 * (g%mask2dCv(i,j)*forces%tauy(i,j))
389 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ endif} \textcolor{comment}{! direct\_stress}
390 
391     \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
392       ray(i,k) = visc%Ray\_v(i,j,k)
393 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
394 
395     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
396       b\_denom\_1 = cs%h\_v(i,j,1) + dt\_z\_to\_h * (ray(i,1) + cs%a\_v(i,j,1))
397       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h*cs%a\_v(i,j,2))
398       d1(i) = b\_denom\_1 * b1(i)
399       v(i,j,1) = b1(i) * (cs%h\_v(i,j,1) * v(i,j,1) + surface\_stress(i))
400 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo}
401     \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
402       c1(i,k) = dt\_z\_to\_h * cs%a\_v(i,j,k) * b1(i)
403       b\_denom\_1 = cs%h\_v(i,j,k) + dt\_z\_to\_h * (ray(i,k) + cs%a\_v(i,j,k)*d1(i))
404       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h * cs%a\_v(i,j,k+1))
405       d1(i) = b\_denom\_1 * b1(i)
406       v(i,j,k) = (cs%h\_v(i,j,k) * v(i,j,k) + dt\_z\_to\_h * cs%a\_v(i,j,k) * v(i,j,k-1)) * b1(i)
407 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
408     \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
409       v(i,j,k) = v(i,j,k) + c1(i,k+1) * v(i,j,k+1)
410 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i and k loops}
411 
412     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(adp%dv\_dt\_visc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
413       adp%dv\_dt\_visc(i,j,k) = (v(i,j,k) - adp%dv\_dt\_visc(i,j,k))*idt
414 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
415 
416     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%tauy\_shelf)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie
417       visc%tauy\_shelf(i,j) = -gv%Rho0*cs%a1\_shelf\_v(i,j)*v(i,j,1) \textcolor{comment}{! - v\_shelf?}
418 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ endif}
419 
420     \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(tauy\_bot)) \textcolor{keywordflow}{then}
421       \textcolor{keywordflow}{do} i=is,ie
422         tauy\_bot(i,j) = gv%Rho0 * (v(i,j,nz)*cs%a\_v(i,j,nz+1))
423 \textcolor{keywordflow}{      enddo}
424       \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
425         tauy\_bot(i,j) = tauy\_bot(i,j) + gv%Rho0 * (ray(i,k)*v(i,j,k))
426 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
427 \textcolor{keywordflow}{    endif}
428 
429     \textcolor{comment}{! When mixing down Eulerian current + Stokes drift subtract after calling solver}
430     \textcolor{keywordflow}{if} (dostokesmixing) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
431       \textcolor{keywordflow}{if} (do\_i(i)) v(i,j,k) = v(i,j,k) - us%m\_s\_to\_L\_T*waves%Us\_y(i,j,k)
432 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
433 
434 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end of v-component J loop}
435 
436   \textcolor{keyword}{call }vertvisc\_limit\_vel(u, v, h, adp, cdp, forces, visc, dt, g, gv, us, cs)
437 
438   \textcolor{comment}{! Here the velocities associated with open boundary conditions are applied.}
439   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then}
440     \textcolor{keywordflow}{do} n=1,obc%number\_of\_segments
441       \textcolor{keywordflow}{if} (obc%segment(n)%specified) \textcolor{keywordflow}{then}
442         \textcolor{keywordflow}{if} (obc%segment(n)%is\_N\_or\_S) \textcolor{keywordflow}{then}
443           j = obc%segment(n)%HI%JsdB
444           \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
445             v(i,j,k) = obc%segment(n)%normal\_vel(i,j,k)
446 \textcolor{keywordflow}{          enddo} ;\textcolor{keywordflow}{ enddo}
447         \textcolor{keywordflow}{elseif} (obc%segment(n)%is\_E\_or\_W) \textcolor{keywordflow}{then}
448           i = obc%segment(n)%HI%IsdB
449           \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
450             u(i,j,k) = obc%segment(n)%normal\_vel(i,j,k)
451 \textcolor{keywordflow}{          enddo} ;\textcolor{keywordflow}{ enddo}
452 \textcolor{keywordflow}{        endif}
453 \textcolor{keywordflow}{      endif}
454 \textcolor{keywordflow}{    enddo}
455 \textcolor{keywordflow}{  endif}
456 
457 \textcolor{comment}{! Offer diagnostic fields for averaging.}
458   \textcolor{keywordflow}{if} (cs%id\_du\_dt\_visc > 0) &
459     \textcolor{keyword}{call }post\_data(cs%id\_du\_dt\_visc, adp%du\_dt\_visc, cs%diag)
460   \textcolor{keywordflow}{if} (cs%id\_dv\_dt\_visc > 0) &
461     \textcolor{keyword}{call }post\_data(cs%id\_dv\_dt\_visc, adp%dv\_dt\_visc, cs%diag)
462   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(taux\_bot) .and. (cs%id\_taux\_bot > 0)) &
463     \textcolor{keyword}{call }post\_data(cs%id\_taux\_bot, taux\_bot, cs%diag)
464   \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(tauy\_bot) .and. (cs%id\_tauy\_bot > 0)) &
465     \textcolor{keyword}{call }post\_data(cs%id\_tauy\_bot, tauy\_bot, cs%diag)
466 
467   \textcolor{comment}{! Diagnostics for terms multiplied by fractional thicknesses}
468 
469   \textcolor{comment}{! 3D diagnostics hf\_du(dv)\_dt\_visc are commented because there is no clarity on proper remapping grid
       option.}
470   \textcolor{comment}{! The code is retained for degugging purposes in the future.}
471   \textcolor{comment}{!if (CS%id\_hf\_du\_dt\_visc > 0) then}
472   \textcolor{comment}{!  do k=1,nz ; do j=js,je ; do I=Isq,Ieq}
473   \textcolor{comment}{!    CS%hf\_du\_dt\_visc(I,j,k) = ADp%du\_dt\_visc(I,j,k) * ADp%diag\_hfrac\_u(I,j,k)}
474   \textcolor{comment}{!  enddo ; enddo ; enddo}
475   \textcolor{comment}{!  call post\_data(CS%id\_hf\_du\_dt\_visc, CS%hf\_du\_dt\_visc, CS%diag)}
476   \textcolor{comment}{!endif}
477   \textcolor{comment}{!if (CS%id\_hf\_dv\_dt\_visc > 0) then}
478   \textcolor{comment}{!  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie}
479   \textcolor{comment}{!    CS%hf\_dv\_dt\_visc(i,J,k) = ADp%dv\_dt\_visc(i,J,k) * ADp%diag\_hfrac\_v(i,J,k)}
480   \textcolor{comment}{!  enddo ; enddo ; enddo}
481   \textcolor{comment}{!  call post\_data(CS%id\_hf\_dv\_dt\_visc, CS%hf\_dv\_dt\_visc, CS%diag)}
482   \textcolor{comment}{!endif}
483   \textcolor{keywordflow}{if} (cs%id\_hf\_du\_dt\_visc\_2d > 0) \textcolor{keywordflow}{then}
484     \textcolor{keyword}{allocate}(hf\_du\_dt\_visc\_2d(g%IsdB:g%IedB,g%jsd:g%jed))
485     hf\_du\_dt\_visc\_2d(:,:) = 0.0
486     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=isq,ieq
487       hf\_du\_dt\_visc\_2d(i,j) = hf\_du\_dt\_visc\_2d(i,j) + adp%du\_dt\_visc(i,j,k) * adp%diag\_hfrac\_u(i,j,k)
488 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
489     \textcolor{keyword}{call }post\_data(cs%id\_hf\_du\_dt\_visc\_2d, hf\_du\_dt\_visc\_2d, cs%diag)
490     \textcolor{keyword}{deallocate}(hf\_du\_dt\_visc\_2d)
491 \textcolor{keywordflow}{  endif}
492   \textcolor{keywordflow}{if} (cs%id\_hf\_dv\_dt\_visc\_2d > 0) \textcolor{keywordflow}{then}
493     \textcolor{keyword}{allocate}(hf\_dv\_dt\_visc\_2d(g%isd:g%ied,g%JsdB:g%JedB))
494     hf\_dv\_dt\_visc\_2d(:,:) = 0.0
495     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} i=is,ie
496       hf\_dv\_dt\_visc\_2d(i,j) = hf\_dv\_dt\_visc\_2d(i,j) + adp%dv\_dt\_visc(i,j,k) * adp%diag\_hfrac\_v(i,j,k)
497 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
498     \textcolor{keyword}{call }post\_data(cs%id\_hf\_dv\_dt\_visc\_2d, hf\_dv\_dt\_visc\_2d, cs%diag)
499     \textcolor{keyword}{deallocate}(hf\_dv\_dt\_visc\_2d)
500 \textcolor{keywordflow}{  endif}
501 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__vert__friction_ac281f6595593b33436594112785e982b}\label{namespacemom__vert__friction_ac281f6595593b33436594112785e982b}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!vertvisc\+\_\+coef@{vertvisc\+\_\+coef}}
\index{vertvisc\+\_\+coef@{vertvisc\+\_\+coef}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{vertvisc\+\_\+coef()}{vertvisc\_coef()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+vert\+\_\+friction\+::vertvisc\+\_\+coef (\begin{DoxyParamCaption}\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(gv)), intent(in)}]{u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(gv)), intent(in)}]{v,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(gv)), intent(in)}]{h,  }\item[{type(mech\+\_\+forcing), intent(in)}]{forces,  }\item[{type(vertvisc\+\_\+type), intent(in)}]{visc,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS,  }\item[{type(ocean\+\_\+obc\+\_\+type), pointer}]{O\+BC }\end{DoxyParamCaption})}



Calculate the coupling coefficients (CSa\+\_\+u and CSa\+\_\+v) and effective layer thicknesses (CSh\+\_\+u and CSh\+\_\+v) for later use in the applying the implicit vertical viscosity via \mbox{\hyperlink{namespacemom__vert__friction_a8f1a390fa24fbe985068fed9ac26873c}{vertvisc()}}. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em u} & Zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em v} & Meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em forces} & A structure with the driving mechanical forces\\
\hline
\mbox{\tt in}  & {\em visc} & Viscosities and bottom drag\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
 & {\em cs} & Vertical viscosity control structure\\
\hline
 & {\em obc} & Open boundary condition structure \\
\hline
\end{DoxyParams}


Definition at line 618 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
618   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: G\textcolor{comment}{      !< Ocean grid structure}
619   \textcolor{keywordtype}{type}(verticalGrid\_type), \textcolor{keywordtype}{intent(in)}    :: GV\textcolor{comment}{     !< Ocean vertical grid structure}
620   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: US\textcolor{comment}{     !< A dimensional unit scaling type}
621   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(GV))}, &
622                            \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{      !< Zonal velocity [L T-1 ~> m s-1]}
623   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(GV))}, &
624                            \textcolor{keywordtype}{intent(in)}    :: v\textcolor{comment}{      !< Meridional velocity [L T-1 ~> m s-1]}
625   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))}, &
626                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{      !< Layer thickness [H ~> m or kg m-2]}
627   \textcolor{keywordtype}{type}(mech\_forcing),      \textcolor{keywordtype}{intent(in)}    :: forces\textcolor{comment}{ !< A structure with the driving mechanical forces}
628   \textcolor{keywordtype}{type}(vertvisc\_type),     \textcolor{keywordtype}{intent(in)}    :: visc\textcolor{comment}{   !< Viscosities and bottom drag}
629   \textcolor{keywordtype}{real},                    \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{     !< Time increment [T ~> s]}
630   \textcolor{keywordtype}{type}(vertvisc\_CS),       \textcolor{keywordtype}{pointer}       :: CS\textcolor{comment}{     !< Vertical viscosity control structure}
631   \textcolor{keywordtype}{type}(ocean\_OBC\_type),    \textcolor{keywordtype}{pointer}       :: OBC\textcolor{comment}{    !< Open boundary condition structure}
632 
633   \textcolor{comment}{! Field from forces used in this subroutine:}
634   \textcolor{comment}{!   ustar: the friction velocity [Z T-1 ~> m s-1], used here as the mixing}
635   \textcolor{comment}{!     velocity in the mixed layer if NKML > 1 in a bulk mixed layer.}
636 
637   \textcolor{comment}{! Local variables}
638 
639   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(G))} :: &
640     h\_harm, &   \textcolor{comment}{! Harmonic mean of the thicknesses around a velocity grid point,}
641                 \textcolor{comment}{! given by 2*(h+ * h-)/(h+ + h-) [H ~> m or kg m-2].}
642     h\_arith, &  \textcolor{comment}{! The arithmetic mean thickness [H ~> m or kg m-2].}
643     h\_delta, &  \textcolor{comment}{! The lateral difference of thickness [H ~> m or kg m-2].}
644     hvel, &     \textcolor{comment}{! hvel is the thickness used at a velocity grid point [H ~> m or kg m-2].}
645     hvel\_shelf  \textcolor{comment}{! The equivalent of hvel under shelves [H ~> m or kg m-2].}
646   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZK\_(G)+1)} :: &
647     a\_cpl, &    \textcolor{comment}{! The drag coefficients across interfaces [Z T-1 ~> m s-1].  a\_cpl times}
648                 \textcolor{comment}{! the velocity difference gives the stress across an interface.}
649     a\_shelf, &  \textcolor{comment}{! The drag coefficients across interfaces in water columns under}
650                 \textcolor{comment}{! ice shelves [Z T-1 ~> m s-1].}
651     z\_i         \textcolor{comment}{! An estimate of each interface's height above the bottom,}
652                 \textcolor{comment}{! normalized by the bottom boundary layer thickness, nondim.}
653   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
654     kv\_bbl, &     \textcolor{comment}{! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].}
655     bbl\_thick, &  \textcolor{comment}{! The bottom boundary layer thickness [H ~> m or kg m-2].}
656     I\_Hbbl, &     \textcolor{comment}{! The inverse of the bottom boundary layer thickness [H-1 ~> m-1 or m2 kg-1].}
657     I\_Htbl, &     \textcolor{comment}{! The inverse of the top boundary layer thickness [H-1 ~> m-1 or m2 kg-1].}
658     zcol1, &      \textcolor{comment}{! The height of the interfaces to the north and south of a}
659     zcol2, &      \textcolor{comment}{! v-point [H ~> m or kg m-2].}
660     Ztop\_min, &   \textcolor{comment}{! The deeper of the two adjacent surface heights [H ~> m or kg m-2].}
661     Dmin, &       \textcolor{comment}{! The shallower of the two adjacent bottom depths converted to}
662                   \textcolor{comment}{! thickness units [H ~> m or kg m-2].}
663     zh, &         \textcolor{comment}{! An estimate of the interface's distance from the bottom}
664                   \textcolor{comment}{! based on harmonic mean thicknesses [H ~> m or kg m-2].}
665     h\_ml          \textcolor{comment}{! The mixed layer depth [H ~> m or kg m-2].}
666   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:,:)} :: hML\_u \textcolor{comment}{! Diagnostic of the mixed layer depth at u points [H ~> m or
       kg m-2].}
667   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:,:)} :: hML\_v \textcolor{comment}{! Diagnostic of the mixed layer depth at v points [H ~> m or
       kg m-2].}
668   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:,:,:)} :: Kv\_u\textcolor{comment}{ !< Total vertical viscosity at u-points [Z2 T-1 ~> m2 s-1].}
669   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{allocatable}, \textcolor{keywordtype}{dimension(:,:,:)} :: Kv\_v\textcolor{comment}{ !< Total vertical viscosity at v-points [Z2 T-1 ~> m2 s-1].}
670   \textcolor{keywordtype}{real} :: zcol(SZI\_(G)) \textcolor{comment}{! The height of an interface at h-points [H ~> m or kg m-2].}
671   \textcolor{keywordtype}{real} :: botfn   \textcolor{comment}{! A function which goes from 1 at the bottom to 0 much more}
672                   \textcolor{comment}{! than Hbbl into the interior.}
673   \textcolor{keywordtype}{real} :: topfn   \textcolor{comment}{! A function which goes from 1 at the top to 0 much more}
674                   \textcolor{comment}{! than Htbl into the interior.}
675   \textcolor{keywordtype}{real} :: z2      \textcolor{comment}{! The distance from the bottom, normalized by Hbbl, nondim.}
676   \textcolor{keywordtype}{real} :: z2\_wt   \textcolor{comment}{! A nondimensional (0-1) weight used when calculating z2.}
677   \textcolor{keywordtype}{real} :: z\_clear \textcolor{comment}{! The clearance of an interface above the surrounding topography [H ~> m or kg m-2].}
678   \textcolor{keywordtype}{real} :: a\_cpl\_max  \textcolor{comment}{! The maximum drag doefficient across interfaces, set so that it will be}
679                      \textcolor{comment}{! representable as a 32-bit float in MKS units  [Z T-1 ~> m s-1]}
680   \textcolor{keywordtype}{real} :: h\_neglect  \textcolor{comment}{! A thickness that is so small it is usually lost}
681                      \textcolor{comment}{! in roundoff and can be neglected [H ~> m or kg m-2].}
682 
683   \textcolor{keywordtype}{real} :: I\_valBL \textcolor{comment}{! The inverse of a scaling factor determining when water is}
684                   \textcolor{comment}{! still within the boundary layer, as determined by the sum}
685                   \textcolor{comment}{! of the harmonic mean thicknesses.}
686   \textcolor{keywordtype}{logical}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: do\_i, do\_i\_shelf
687   \textcolor{keywordtype}{logical} :: do\_any\_shelf
688   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{dimension(SZIB\_(G))} :: &
689     zi\_dir   \textcolor{comment}{!  A trinary logical array indicating which thicknesses to use for}
690              \textcolor{comment}{!  finding z\_clear.}
691   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
692   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
693   isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
694 
695   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal,\textcolor{stringliteral}{"MOM\_vert\_friction(coef): "}// &
696          \textcolor{stringliteral}{"Module must be initialized before it is used."})
697 
698   h\_neglect = gv%H\_subroundoff
699   a\_cpl\_max = 1.0e37 * us%m\_to\_Z * us%T\_to\_s
700   i\_hbbl(:) = 1.0 / (cs%Hbbl + h\_neglect)
701   i\_valbl = 0.0 ; \textcolor{keywordflow}{if} (cs%harm\_BL\_val > 0.0) i\_valbl = 1.0 / cs%harm\_BL\_val
702 
703   \textcolor{keywordflow}{if} (cs%id\_Kv\_u > 0) \textcolor{keywordflow}{then}
704     \textcolor{keyword}{allocate}(kv\_u(g%IsdB:g%IedB,g%jsd:g%jed,g%ke)) ; kv\_u(:,:,:) = 0.0
705 \textcolor{keywordflow}{  endif}
706 
707   \textcolor{keywordflow}{if} (cs%id\_Kv\_v > 0) \textcolor{keywordflow}{then}
708     \textcolor{keyword}{allocate}(kv\_v(g%isd:g%ied,g%JsdB:g%JedB,g%ke)) ; kv\_v(:,:,:) = 0.0
709 \textcolor{keywordflow}{  endif}
710 
711   \textcolor{keywordflow}{if} (cs%debug .or. (cs%id\_hML\_u > 0)) \textcolor{keywordflow}{then}
712     \textcolor{keyword}{allocate}(hml\_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; hml\_u(:,:) = 0.0
713 \textcolor{keywordflow}{  endif}
714   \textcolor{keywordflow}{if} (cs%debug .or. (cs%id\_hML\_v > 0)) \textcolor{keywordflow}{then}
715     \textcolor{keyword}{allocate}(hml\_v(g%isd:g%ied,g%JsdB:g%JedB)) ; hml\_v(:,:) = 0.0
716 \textcolor{keywordflow}{  endif}
717 
718   \textcolor{keywordflow}{if} ((\textcolor{keyword}{associated}(visc%taux\_shelf) .or. \textcolor{keyword}{associated}(forces%frac\_shelf\_u)) .and. &
719       .not.\textcolor{keyword}{associated}(cs%a1\_shelf\_u)) \textcolor{keywordflow}{then}
720     \textcolor{keyword}{allocate}(cs%a1\_shelf\_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; cs%a1\_shelf\_u(:,:)=0.0
721 \textcolor{keywordflow}{  endif}
722   \textcolor{keywordflow}{if} ((\textcolor{keyword}{associated}(visc%tauy\_shelf) .or. \textcolor{keyword}{associated}(forces%frac\_shelf\_v)) .and. &
723       .not.\textcolor{keyword}{associated}(cs%a1\_shelf\_v)) \textcolor{keywordflow}{then}
724     \textcolor{keyword}{allocate}(cs%a1\_shelf\_v(g%isd:g%ied,g%JsdB:g%JedB)) ; cs%a1\_shelf\_v(:,:)=0.0
725 \textcolor{keywordflow}{  endif}
726 
727   \textcolor{comment}{!$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML\_u, &}
728   \textcolor{comment}{!$OMP                                     OBC,h\_neglect,dt,I\_valBL,Kv\_u,a\_cpl\_max) &}
729   \textcolor{comment}{!$OMP                     firstprivate(i\_hbbl)}
730   \textcolor{keywordflow}{do} j=g%Jsc,g%Jec
731     \textcolor{keywordflow}{do} i=isq,ieq ; do\_i(i) = (g%mask2dCu(i,j) > 0) ;\textcolor{keywordflow}{ enddo}
732 
733     \textcolor{keywordflow}{if} (cs%bottomdraglaw) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=isq,ieq
734       kv\_bbl(i) = visc%Kv\_bbl\_u(i,j)
735       bbl\_thick(i) = visc%bbl\_thick\_u(i,j) * gv%Z\_to\_H + h\_neglect
736       \textcolor{keywordflow}{if} (do\_i(i)) i\_hbbl(i) = 1.0 / bbl\_thick(i)
737 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ endif}
738 
739     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
740       h\_harm(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h\_neglect)
741       h\_arith(i,k) = 0.5*(h(i+1,j,k)+h(i,j,k))
742       h\_delta(i,k) = h(i+1,j,k) - h(i,j,k)
743 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
744     \textcolor{keywordflow}{do} i=isq,ieq
745       dmin(i) = min(g%bathyT(i,j), g%bathyT(i+1,j)) * gv%Z\_to\_H
746       zi\_dir(i) = 0
747 \textcolor{keywordflow}{    enddo}
748 
749     \textcolor{comment}{! Project thickness outward across OBCs using a zero-gradient condition.}
750     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (obc%number\_of\_segments > 0) \textcolor{keywordflow}{then}
751       \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_u(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
752         \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_e) \textcolor{keywordflow}{then}
753           \textcolor{keywordflow}{do} k=1,nz ; h\_harm(i,k) = h(i,j,k) ; h\_arith(i,k) = h(i,j,k) ; h\_delta(i,k) = 0. ;\textcolor{keywordflow}{ enddo}
754           dmin(i) = g%bathyT(i,j) * gv%Z\_to\_H
755           zi\_dir(i) = -1
756         \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_u(i,j))%direction == obc\_direction\_w) \textcolor{keywordflow}{then}
757           \textcolor{keywordflow}{do} k=1,nz ; h\_harm(i,k) = h(i+1,j,k) ; h\_arith(i,k) = h(i+1,j,k) ; h\_delta(i,k) = 0. ;\textcolor{keywordflow}{ enddo}
758           dmin(i) = g%bathyT(i+1,j) * gv%Z\_to\_H
759           zi\_dir(i) = 1
760 \textcolor{keywordflow}{        endif}
761 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
762 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
763 
764 \textcolor{comment}{!    The following block calculates the thicknesses at velocity}
765 \textcolor{comment}{!  grid points for the vertical viscosity (hvel).  Near the}
766 \textcolor{comment}{!  bottom an upwind biased thickness is used to control the effect}
767 \textcolor{comment}{!  of spurious Montgomery potential gradients at the bottom where}
768 \textcolor{comment}{!  nearly massless layers layers ride over the topography.}
769     \textcolor{keywordflow}{if} (cs%harmonic\_visc) \textcolor{keywordflow}{then}
770       \textcolor{keywordflow}{do} i=isq,ieq ; z\_i(i,nz+1) = 0.0 ;\textcolor{keywordflow}{ enddo}
771       \textcolor{keywordflow}{do} k=nz,1,-1 ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
772         hvel(i,k) = h\_harm(i,k)
773         \textcolor{keywordflow}{if} (u(i,j,k) * h\_delta(i,k) < 0) \textcolor{keywordflow}{then}
774           z2 = z\_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
775           hvel(i,k) = (1.0-botfn)*h\_harm(i,k) + botfn*h\_arith(i,k)
776 \textcolor{keywordflow}{        endif}
777         z\_i(i,k) =  z\_i(i,k+1) + h\_harm(i,k)*i\_hbbl(i)
778 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i & k loops}
779     \textcolor{keywordflow}{else} \textcolor{comment}{! Not harmonic\_visc}
780       \textcolor{keywordflow}{do} i=isq,ieq ; zh(i) = 0.0 ; z\_i(i,nz+1) = 0.0 ;\textcolor{keywordflow}{ enddo}
781       \textcolor{keywordflow}{do} i=isq,ieq+1 ; zcol(i) = -g%bathyT(i,j) * gv%Z\_to\_H ;\textcolor{keywordflow}{ enddo}
782       \textcolor{keywordflow}{do} k=nz,1,-1
783         \textcolor{keywordflow}{do} i=isq,ieq+1 ; zcol(i) = zcol(i) + h(i,j,k) ;\textcolor{keywordflow}{ enddo}
784         \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
785           zh(i) = zh(i) + h\_harm(i,k)
786 
787           z\_clear = max(zcol(i),zcol(i+1)) + dmin(i)
788           \textcolor{keywordflow}{if} (zi\_dir(i) < 0) z\_clear = zcol(i) + dmin(i)
789           \textcolor{keywordflow}{if} (zi\_dir(i) > 0) z\_clear = zcol(i+1) + dmin(i)
790 
791           z\_i(i,k) = max(zh(i), z\_clear) * i\_hbbl(i)
792 
793           hvel(i,k) = h\_arith(i,k)
794           \textcolor{keywordflow}{if} (u(i,j,k) * h\_delta(i,k) > 0) \textcolor{keywordflow}{then}
795             \textcolor{keywordflow}{if} (zh(i) * i\_hbbl(i) < cs%harm\_BL\_val) \textcolor{keywordflow}{then}
796               hvel(i,k) = h\_harm(i,k)
797             \textcolor{keywordflow}{else}
798               z2\_wt = 1.0  ; \textcolor{keywordflow}{if} (zh(i) * i\_hbbl(i) < 2.0*cs%harm\_BL\_val) &
799                 z2\_wt = max(0.0, min(1.0, zh(i) * i\_hbbl(i) * i\_valbl - 1.0))
800               z2 = z2\_wt * (max(zh(i), z\_clear) * i\_hbbl(i))
801               botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
802               hvel(i,k) = (1.0-botfn)*h\_arith(i,k) + botfn*h\_harm(i,k)
803 \textcolor{keywordflow}{            endif}
804 \textcolor{keywordflow}{          endif}
805 
806 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i  loop}
807 \textcolor{keywordflow}{      enddo} \textcolor{comment}{! k loop}
808 \textcolor{keywordflow}{    endif}
809 
810     \textcolor{keyword}{call }find\_coupling\_coef(a\_cpl, hvel, do\_i, h\_harm, bbl\_thick, kv\_bbl, z\_i, h\_ml, &
811                             dt, j, g, gv, us, cs, visc, forces, work\_on\_u=.true., obc=obc)
812     \textcolor{keywordflow}{if} (\textcolor{keyword}{allocated}(hml\_u)) \textcolor{keywordflow}{then}
813       \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then} ; hml\_u(i,j) = h\_ml(i) ;\textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo}
814 \textcolor{keywordflow}{    endif}
815 
816     do\_any\_shelf = .false.
817     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(forces%frac\_shelf\_u)) \textcolor{keywordflow}{then}
818       \textcolor{keywordflow}{do} i=isq,ieq
819         cs%a1\_shelf\_u(i,j) = 0.0
820         do\_i\_shelf(i) = (do\_i(i) .and. forces%frac\_shelf\_u(i,j) > 0.0)
821         \textcolor{keywordflow}{if} (do\_i\_shelf(i)) do\_any\_shelf = .true.
822 \textcolor{keywordflow}{      enddo}
823       \textcolor{keywordflow}{if} (do\_any\_shelf) \textcolor{keywordflow}{then}
824         \textcolor{keywordflow}{if} (cs%harmonic\_visc) \textcolor{keywordflow}{then}
825           \textcolor{keyword}{call }find\_coupling\_coef(a\_shelf, hvel, do\_i\_shelf, h\_harm, bbl\_thick, &
826                                   kv\_bbl, z\_i, h\_ml, dt, j, g, gv, us, cs, &
827                                   visc, forces, work\_on\_u=.true., obc=obc, shelf=.true.)
828         \textcolor{keywordflow}{else}  \textcolor{comment}{! Find upwind-biased thickness near the surface.}
829           \textcolor{comment}{! Perhaps this needs to be done more carefully, via find\_eta.}
830           \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
831             zh(i) = 0.0 ; ztop\_min(i) = min(zcol(i), zcol(i+1))
832             i\_htbl(i) = 1.0 / (visc%tbl\_thick\_shelf\_u(i,j)*gv%Z\_to\_H + h\_neglect)
833 \textcolor{keywordflow}{          endif} ;\textcolor{keywordflow}{ enddo}
834           \textcolor{keywordflow}{do} k=1,nz
835             \textcolor{keywordflow}{do} i=isq,ieq+1 ; zcol(i) = zcol(i) - h(i,j,k) ;\textcolor{keywordflow}{ enddo}
836             \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
837               zh(i) = zh(i) + h\_harm(i,k)
838 
839               hvel\_shelf(i,k) = hvel(i,k)
840               \textcolor{keywordflow}{if} (u(i,j,k) * h\_delta(i,k) > 0) \textcolor{keywordflow}{then}
841                 \textcolor{keywordflow}{if} (zh(i) * i\_htbl(i) < cs%harm\_BL\_val) \textcolor{keywordflow}{then}
842                   hvel\_shelf(i,k) = min(hvel(i,k), h\_harm(i,k))
843                 \textcolor{keywordflow}{else}
844                   z2\_wt = 1.0  ; \textcolor{keywordflow}{if} (zh(i) * i\_htbl(i) < 2.0*cs%harm\_BL\_val) &
845                     z2\_wt = max(0.0, min(1.0, zh(i) * i\_htbl(i) * i\_valbl - 1.0))
846                   z2 = z2\_wt * (max(zh(i), ztop\_min(i) - min(zcol(i),zcol(i+1))) * i\_htbl(i))
847                   topfn = 1.0 / (1.0 + 0.09*z2**6)
848                   hvel\_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h\_arith(i,k) + topfn*h\_harm(i,k))
849 \textcolor{keywordflow}{                endif}
850 \textcolor{keywordflow}{              endif}
851 \textcolor{keywordflow}{            endif} ;\textcolor{keywordflow}{ enddo}
852 \textcolor{keywordflow}{          enddo}
853           \textcolor{keyword}{call }find\_coupling\_coef(a\_shelf, hvel\_shelf, do\_i\_shelf, h\_harm, &
854                                   bbl\_thick, kv\_bbl, z\_i, h\_ml, dt, j, g, gv, us, cs, &
855                                   visc, forces, work\_on\_u=.true., obc=obc, shelf=.true.)
856 \textcolor{keywordflow}{        endif}
857         \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) cs%a1\_shelf\_u(i,j) = a\_shelf(i,1) ;\textcolor{keywordflow}{ enddo}
858 \textcolor{keywordflow}{      endif}
859 \textcolor{keywordflow}{    endif}
860 
861     \textcolor{keywordflow}{if} (do\_any\_shelf) \textcolor{keywordflow}{then}
862       \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
863         cs%a\_u(i,j,k) = min(a\_cpl\_max, forces%frac\_shelf\_u(i,j)  * a\_shelf(i,k) + &
864                                        (1.0-forces%frac\_shelf\_u(i,j)) * a\_cpl(i,k))
865 \textcolor{comment}{! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH}
866 \textcolor{comment}{!        CS%a\_u(I,j,K) = min(a\_cpl\_max, forces%frac\_shelf\_u(I,j)  * max(a\_shelf(I,K), a\_cpl(I,K)) + &}
867 \textcolor{comment}{!                                       (1.0-forces%frac\_shelf\_u(I,j)) * a\_cpl(I,K))}
868       \textcolor{keywordflow}{elseif} (do\_i(i)) \textcolor{keywordflow}{then}
869         cs%a\_u(i,j,k) = min(a\_cpl\_max, a\_cpl(i,k))
870 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
871       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
872         \textcolor{comment}{! Should we instead take the inverse of the average of the inverses?}
873         cs%h\_u(i,j,k) = forces%frac\_shelf\_u(i,j)  * hvel\_shelf(i,k) + &
874                    (1.0-forces%frac\_shelf\_u(i,j)) * hvel(i,k) + h\_neglect
875       \textcolor{keywordflow}{elseif} (do\_i(i)) \textcolor{keywordflow}{then}
876         cs%h\_u(i,j,k) = hvel(i,k) + h\_neglect
877 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
878     \textcolor{keywordflow}{else}
879       \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) cs%a\_u(i,j,k) = min(a\_cpl\_max, a\_cpl(i,k)) ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
880       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) cs%h\_u(i,j,k) = hvel(i,k) + h\_neglect ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
881 \textcolor{keywordflow}{    endif}
882 
883     \textcolor{comment}{! Diagnose total Kv at u-points}
884     \textcolor{keywordflow}{if} (cs%id\_Kv\_u > 0) \textcolor{keywordflow}{then}
885       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
886         \textcolor{keywordflow}{if} (do\_i(i)) kv\_u(i,j,k) = 0.5 * gv%H\_to\_Z*(cs%a\_u(i,j,k)+cs%a\_u(i,j,k+1)) * cs%h\_u(i,j,k)
887 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
888 \textcolor{keywordflow}{    endif}
889 
890 \textcolor{keywordflow}{  enddo}
891 
892 
893   \textcolor{comment}{! Now work on v-points.}
894   \textcolor{comment}{!$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML\_v, &}
895   \textcolor{comment}{!$OMP                                  OBC,h\_neglect,dt,I\_valBL,Kv\_v,a\_cpl\_max) &}
896   \textcolor{comment}{!$OMP                     firstprivate(i\_hbbl)}
897   \textcolor{keywordflow}{do} j=jsq,jeq
898     \textcolor{keywordflow}{do} i=is,ie ; do\_i(i) = (g%mask2dCv(i,j) > 0) ;\textcolor{keywordflow}{ enddo}
899 
900     \textcolor{keywordflow}{if} (cs%bottomdraglaw) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} i=is,ie
901       kv\_bbl(i) = visc%Kv\_bbl\_v(i,j)
902       bbl\_thick(i) = visc%bbl\_thick\_v(i,j) * gv%Z\_to\_H + h\_neglect
903       \textcolor{keywordflow}{if} (do\_i(i)) i\_hbbl(i) = 1.0 / bbl\_thick(i)
904 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ endif}
905 
906     \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
907       h\_harm(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / (h(i,j,k)+h(i,j+1,k)+h\_neglect)
908       h\_arith(i,k) = 0.5*(h(i,j+1,k)+h(i,j,k))
909       h\_delta(i,k) = h(i,j+1,k) - h(i,j,k)
910 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
911     \textcolor{keywordflow}{do} i=is,ie
912       dmin(i) = min(g%bathyT(i,j), g%bathyT(i,j+1)) * gv%Z\_to\_H
913       zi\_dir(i) = 0
914 \textcolor{keywordflow}{    enddo}
915 
916     \textcolor{comment}{! Project thickness outward across OBCs using a zero-gradient condition.}
917     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(obc)) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (obc%number\_of\_segments > 0) \textcolor{keywordflow}{then}
918       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i) .and. (obc%segnum\_v(i,j) /= obc\_none)) \textcolor{keywordflow}{then}
919         \textcolor{keywordflow}{if} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_n) \textcolor{keywordflow}{then}
920           \textcolor{keywordflow}{do} k=1,nz ; h\_harm(i,k) = h(i,j,k) ; h\_arith(i,k) = h(i,j,k) ; h\_delta(i,k) = 0. ;\textcolor{keywordflow}{ enddo}
921           dmin(i) = g%bathyT(i,j) * gv%Z\_to\_H
922           zi\_dir(i) = -1
923         \textcolor{keywordflow}{elseif} (obc%segment(obc%segnum\_v(i,j))%direction == obc\_direction\_s) \textcolor{keywordflow}{then}
924           \textcolor{keywordflow}{do} k=1,nz ; h\_harm(i,k) = h(i,j+1,k) ; h\_arith(i,k) = h(i,j+1,k) ; h\_delta(i,k) = 0. ;\textcolor{keywordflow}{ enddo}
925           dmin(i) = g%bathyT(i,j+1) * gv%Z\_to\_H
926           zi\_dir(i) = 1
927 \textcolor{keywordflow}{        endif}
928 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
929 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ endif}
930 
931 \textcolor{comment}{!    The following block calculates the thicknesses at velocity}
932 \textcolor{comment}{!  grid points for the vertical viscosity (hvel).  Near the}
933 \textcolor{comment}{!  bottom an upwind biased thickness is used to control the effect}
934 \textcolor{comment}{!  of spurious Montgomery potential gradients at the bottom where}
935 \textcolor{comment}{!  nearly massless layers layers ride over the topography.}
936     \textcolor{keywordflow}{if} (cs%harmonic\_visc) \textcolor{keywordflow}{then}
937       \textcolor{keywordflow}{do} i=is,ie ; z\_i(i,nz+1) = 0.0 ;\textcolor{keywordflow}{ enddo}
938 
939       \textcolor{keywordflow}{do} k=nz,1,-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
940         hvel(i,k) = h\_harm(i,k)
941         \textcolor{keywordflow}{if} (v(i,j,k) * h\_delta(i,k) < 0) \textcolor{keywordflow}{then}
942           z2 = z\_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
943           hvel(i,k) = (1.0-botfn)*h\_harm(i,k) + botfn*h\_arith(i,k)
944 \textcolor{keywordflow}{        endif}
945         z\_i(i,k) = z\_i(i,k+1)  + h\_harm(i,k)*i\_hbbl(i)
946 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i & k loops}
947     \textcolor{keywordflow}{else} \textcolor{comment}{! Not harmonic\_visc}
948       \textcolor{keywordflow}{do} i=is,ie
949         zh(i) = 0.0 ; z\_i(i,nz+1) = 0.0
950         zcol1(i) = -g%bathyT(i,j) * gv%Z\_to\_H
951         zcol2(i) = -g%bathyT(i,j+1) * gv%Z\_to\_H
952 \textcolor{keywordflow}{      enddo}
953       \textcolor{keywordflow}{do} k=nz,1,-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
954         zh(i) = zh(i) + h\_harm(i,k)
955         zcol1(i) = zcol1(i) + h(i,j,k) ; zcol2(i) = zcol2(i) + h(i,j+1,k)
956 
957         z\_clear = max(zcol1(i),zcol2(i)) + dmin(i)
958         \textcolor{keywordflow}{if} (zi\_dir(i) < 0) z\_clear = zcol1(i) + dmin(i)
959         \textcolor{keywordflow}{if} (zi\_dir(i) > 0) z\_clear = zcol2(i) + dmin(i)
960 
961         z\_i(i,k) = max(zh(i), z\_clear) * i\_hbbl(i)
962 
963         hvel(i,k) = h\_arith(i,k)
964         \textcolor{keywordflow}{if} (v(i,j,k) * h\_delta(i,k) > 0) \textcolor{keywordflow}{then}
965           \textcolor{keywordflow}{if} (zh(i) * i\_hbbl(i) < cs%harm\_BL\_val) \textcolor{keywordflow}{then}
966             hvel(i,k) = h\_harm(i,k)
967           \textcolor{keywordflow}{else}
968             z2\_wt = 1.0  ; \textcolor{keywordflow}{if} (zh(i) * i\_hbbl(i) < 2.0*cs%harm\_BL\_val) &
969               z2\_wt = max(0.0, min(1.0, zh(i) * i\_hbbl(i) * i\_valbl - 1.0))
970             z2 = z2\_wt * (max(zh(i), max(zcol1(i),zcol2(i)) + dmin(i)) * i\_hbbl(i))
971             botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
972             hvel(i,k) = (1.0-botfn)*h\_arith(i,k) + botfn*h\_harm(i,k)
973 \textcolor{keywordflow}{          endif}
974 \textcolor{keywordflow}{       endif}
975 
976 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i & k loops}
977 \textcolor{keywordflow}{    endif}
978 
979     \textcolor{keyword}{call }find\_coupling\_coef(a\_cpl, hvel, do\_i, h\_harm, bbl\_thick, kv\_bbl, z\_i, h\_ml, &
980                             dt, j, g, gv, us, cs, visc, forces, work\_on\_u=.false., obc=obc)
981     \textcolor{keywordflow}{if} ( \textcolor{keyword}{allocated}(hml\_v)) \textcolor{keywordflow}{then}
982        \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then} ; hml\_v(i,j) = h\_ml(i) ;\textcolor{keywordflow}{ endif} ;\textcolor{keywordflow}{ enddo}
983 \textcolor{keywordflow}{    endif}
984     do\_any\_shelf = .false.
985     \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(forces%frac\_shelf\_v)) \textcolor{keywordflow}{then}
986       \textcolor{keywordflow}{do} i=is,ie
987         cs%a1\_shelf\_v(i,j) = 0.0
988         do\_i\_shelf(i) = (do\_i(i) .and. forces%frac\_shelf\_v(i,j) > 0.0)
989         \textcolor{keywordflow}{if} (do\_i\_shelf(i)) do\_any\_shelf = .true.
990 \textcolor{keywordflow}{      enddo}
991       \textcolor{keywordflow}{if} (do\_any\_shelf) \textcolor{keywordflow}{then}
992         \textcolor{keywordflow}{if} (cs%harmonic\_visc) \textcolor{keywordflow}{then}
993           \textcolor{keyword}{call }find\_coupling\_coef(a\_shelf, hvel, do\_i\_shelf, h\_harm, bbl\_thick, &
994                                   kv\_bbl, z\_i, h\_ml, dt, j, g, gv, us, cs, visc, &
995                                   forces, work\_on\_u=.false., obc=obc, shelf=.true.)
996         \textcolor{keywordflow}{else}  \textcolor{comment}{! Find upwind-biased thickness near the surface.}
997           \textcolor{comment}{! Perhaps this needs to be done more carefully, via find\_eta.}
998           \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
999             zh(i) = 0.0 ; ztop\_min(i) = min(zcol1(i), zcol2(i))
1000             i\_htbl(i) = 1.0 / (visc%tbl\_thick\_shelf\_v(i,j)*gv%Z\_to\_H + h\_neglect)
1001 \textcolor{keywordflow}{          endif} ;\textcolor{keywordflow}{ enddo}
1002           \textcolor{keywordflow}{do} k=1,nz
1003             \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
1004               zcol1(i) = zcol1(i) - h(i,j,k) ; zcol2(i) = zcol2(i) - h(i,j+1,k)
1005               zh(i) = zh(i) + h\_harm(i,k)
1006 
1007               hvel\_shelf(i,k) = hvel(i,k)
1008               \textcolor{keywordflow}{if} (v(i,j,k) * h\_delta(i,k) > 0) \textcolor{keywordflow}{then}
1009                 \textcolor{keywordflow}{if} (zh(i) * i\_htbl(i) < cs%harm\_BL\_val) \textcolor{keywordflow}{then}
1010                   hvel\_shelf(i,k) = min(hvel(i,k), h\_harm(i,k))
1011                 \textcolor{keywordflow}{else}
1012                   z2\_wt = 1.0  ; \textcolor{keywordflow}{if} (zh(i) * i\_htbl(i) < 2.0*cs%harm\_BL\_val) &
1013                     z2\_wt = max(0.0, min(1.0, zh(i) * i\_htbl(i) * i\_valbl - 1.0))
1014                   z2 = z2\_wt * (max(zh(i), ztop\_min(i) - min(zcol1(i),zcol2(i))) * i\_htbl(i))
1015                   topfn = 1.0 / (1.0 + 0.09*z2**6)
1016                   hvel\_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h\_arith(i,k) + topfn*h\_harm(i,k))
1017 \textcolor{keywordflow}{                endif}
1018 \textcolor{keywordflow}{             endif}
1019 \textcolor{keywordflow}{            endif} ;\textcolor{keywordflow}{ enddo}
1020 \textcolor{keywordflow}{          enddo}
1021           \textcolor{keyword}{call }find\_coupling\_coef(a\_shelf, hvel\_shelf, do\_i\_shelf, h\_harm, &
1022                                   bbl\_thick, kv\_bbl, z\_i, h\_ml, dt, j, g, gv, us, cs, &
1023                                   visc, forces, work\_on\_u=.false., obc=obc, shelf=.true.)
1024 \textcolor{keywordflow}{        endif}
1025         \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) cs%a1\_shelf\_v(i,j) = a\_shelf(i,1) ;\textcolor{keywordflow}{ enddo}
1026 \textcolor{keywordflow}{      endif}
1027 \textcolor{keywordflow}{    endif}
1028 
1029     \textcolor{keywordflow}{if} (do\_any\_shelf) \textcolor{keywordflow}{then}
1030       \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
1031         cs%a\_v(i,j,k) = min(a\_cpl\_max, forces%frac\_shelf\_v(i,j)  * a\_shelf(i,k) + &
1032                                        (1.0-forces%frac\_shelf\_v(i,j)) * a\_cpl(i,k))
1033 \textcolor{comment}{! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH}
1034 \textcolor{comment}{!        CS%a\_v(i,J,K) = min(a\_cpl\_max, forces%frac\_shelf\_v(i,J)  * max(a\_shelf(i,K), a\_cpl(i,K)) + &}
1035                     \textcolor{comment}{!                   (1.0-forces%frac\_shelf\_v(i,J)) * a\_cpl(i,K))}
1036       \textcolor{keywordflow}{elseif} (do\_i(i)) \textcolor{keywordflow}{then}
1037         cs%a\_v(i,j,k) = min(a\_cpl\_max, a\_cpl(i,k))
1038 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1039       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i\_shelf(i)) \textcolor{keywordflow}{then}
1040         \textcolor{comment}{! Should we instead take the inverse of the average of the inverses?}
1041         cs%h\_v(i,j,k) = forces%frac\_shelf\_v(i,j)  * hvel\_shelf(i,k) + &
1042                    (1.0-forces%frac\_shelf\_v(i,j)) * hvel(i,k) + h\_neglect
1043       \textcolor{keywordflow}{elseif} (do\_i(i)) \textcolor{keywordflow}{then}
1044         cs%h\_v(i,j,k) = hvel(i,k) + h\_neglect
1045 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1046     \textcolor{keywordflow}{else}
1047       \textcolor{keywordflow}{do} k=1,nz+1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) cs%a\_v(i,j,k) = min(a\_cpl\_max, a\_cpl(i,k)) ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1048       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) cs%h\_v(i,j,k) = hvel(i,k) + h\_neglect ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1049 \textcolor{keywordflow}{    endif}
1050 
1051     \textcolor{comment}{! Diagnose total Kv at v-points}
1052     \textcolor{keywordflow}{if} (cs%id\_Kv\_v > 0) \textcolor{keywordflow}{then}
1053       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
1054         \textcolor{keywordflow}{if} (do\_i(i)) kv\_v(i,j,k) = 0.5 * gv%H\_to\_Z*(cs%a\_v(i,j,k)+cs%a\_v(i,j,k+1)) * cs%h\_v(i,j,k)
1055 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo}
1056 \textcolor{keywordflow}{    endif}
1057 
1058 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end of v-point j loop}
1059 
1060   \textcolor{keywordflow}{if} (cs%debug) \textcolor{keywordflow}{then}
1061     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"vertvisc\_coef h\_[uv]"}, cs%h\_u, cs%h\_v, g%HI, haloshift=0, &
1062                   scale=gv%H\_to\_m, scalar\_pair=.true.)
1063     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"vertvisc\_coef a\_[uv]"}, cs%a\_u, cs%a\_v, g%HI, haloshift=0, &
1064                   scale=us%Z\_to\_m*us%s\_to\_T, scalar\_pair=.true.)
1065     \textcolor{keywordflow}{if} (\textcolor{keyword}{allocated}(hml\_u) .and. \textcolor{keyword}{allocated}(hml\_v)) &
1066       \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"vertvisc\_coef hML\_[uv]"}, hml\_u, hml\_v, g%HI, &
1067                     haloshift=0, scale=gv%H\_to\_m, scalar\_pair=.true.)
1068 \textcolor{keywordflow}{  endif}
1069 
1070 \textcolor{comment}{! Offer diagnostic fields for averaging.}
1071   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(visc%Kv\_slow) .and. (cs%id\_Kv\_slow > 0)) &
1072       \textcolor{keyword}{call }post\_data(cs%id\_Kv\_slow, visc%Kv\_slow, cs%diag)
1073   \textcolor{keywordflow}{if} (cs%id\_Kv\_u > 0) \textcolor{keyword}{call }post\_data(cs%id\_Kv\_u, kv\_u, cs%diag)
1074   \textcolor{keywordflow}{if} (cs%id\_Kv\_v > 0) \textcolor{keyword}{call }post\_data(cs%id\_Kv\_v, kv\_v, cs%diag)
1075   \textcolor{keywordflow}{if} (cs%id\_au\_vv > 0) \textcolor{keyword}{call }post\_data(cs%id\_au\_vv, cs%a\_u, cs%diag)
1076   \textcolor{keywordflow}{if} (cs%id\_av\_vv > 0) \textcolor{keyword}{call }post\_data(cs%id\_av\_vv, cs%a\_v, cs%diag)
1077   \textcolor{keywordflow}{if} (cs%id\_h\_u > 0) \textcolor{keyword}{call }post\_data(cs%id\_h\_u, cs%h\_u, cs%diag)
1078   \textcolor{keywordflow}{if} (cs%id\_h\_v > 0) \textcolor{keyword}{call }post\_data(cs%id\_h\_v, cs%h\_v, cs%diag)
1079   \textcolor{keywordflow}{if} (cs%id\_hML\_u > 0) \textcolor{keyword}{call }post\_data(cs%id\_hML\_u, hml\_u, cs%diag)
1080   \textcolor{keywordflow}{if} (cs%id\_hML\_v > 0) \textcolor{keyword}{call }post\_data(cs%id\_hML\_v, hml\_v, cs%diag)
1081 
1082   \textcolor{keywordflow}{if} (\textcolor{keyword}{allocated}(hml\_u)) \textcolor{keyword}{deallocate}(hml\_u)
1083   \textcolor{keywordflow}{if} (\textcolor{keyword}{allocated}(hml\_v)) \textcolor{keyword}{deallocate}(hml\_v)
1084 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__vert__friction_a158d29cf5c1d79ca7c5f327706855972}\label{namespacemom__vert__friction_a158d29cf5c1d79ca7c5f327706855972}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!vertvisc\+\_\+end@{vertvisc\+\_\+end}}
\index{vertvisc\+\_\+end@{vertvisc\+\_\+end}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{vertvisc\+\_\+end()}{vertvisc\_end()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+vert\+\_\+friction\+::vertvisc\+\_\+end (\begin{DoxyParamCaption}\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS }\end{DoxyParamCaption})}



Clean up and deallocate the vertical friction module. 


\begin{DoxyParams}{Parameters}
{\em cs} & Vertical viscosity control structure that will be deallocated in this subroutine. \\
\hline
\end{DoxyParams}


Definition at line 1890 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
1890   \textcolor{keywordtype}{type}(vertvisc\_CS), \textcolor{keywordtype}{pointer} :: CS\textcolor{comment}{ !< Vertical viscosity control structure that}
1891 \textcolor{comment}{                                   !! will be deallocated in this subroutine.}
1892 
1893   dealloc\_(cs%a\_u) ; dealloc\_(cs%h\_u)
1894   dealloc\_(cs%a\_v) ; dealloc\_(cs%h\_v)
1895   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs%a1\_shelf\_u)) \textcolor{keyword}{deallocate}(cs%a1\_shelf\_u)
1896   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs%a1\_shelf\_v)) \textcolor{keyword}{deallocate}(cs%a1\_shelf\_v)
1897   \textcolor{keyword}{deallocate}(cs)
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__vert__friction_a4bd5c8772584e41890bff55ccd52507c}\label{namespacemom__vert__friction_a4bd5c8772584e41890bff55ccd52507c}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!vertvisc\+\_\+init@{vertvisc\+\_\+init}}
\index{vertvisc\+\_\+init@{vertvisc\+\_\+init}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{vertvisc\+\_\+init()}{vertvisc\_init()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+vert\+\_\+friction\+::vertvisc\+\_\+init (\begin{DoxyParamCaption}\item[{type(ocean\+\_\+internal\+\_\+state), intent(in), target}]{M\+IS,  }\item[{type(time\+\_\+type), intent(in), target}]{Time,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(param\+\_\+file\+\_\+type), intent(in)}]{param\+\_\+file,  }\item[{type(diag\+\_\+ctrl), intent(inout), target}]{diag,  }\item[{type(accel\+\_\+diag\+\_\+ptrs), intent(inout)}]{A\+Dp,  }\item[{type(directories), intent(in)}]{dirs,  }\item[{integer, intent(inout), target}]{ntrunc,  }\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS }\end{DoxyParamCaption})}



Initialize the vertical friction module. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em mis} & The \char`\"{}\+M\+O\+M Internal State\char`\"{}, a set of pointers\\
\hline
\mbox{\tt in}  & {\em time} & Current model time\\
\hline
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in}  & {\em param\+\_\+file} & File to parse for parameters\\
\hline
\mbox{\tt in,out}  & {\em diag} & Diagnostic control structure\\
\hline
\mbox{\tt in,out}  & {\em adp} & Acceleration diagnostic pointers\\
\hline
\mbox{\tt in}  & {\em dirs} & Relevant directory paths\\
\hline
\mbox{\tt in,out}  & {\em ntrunc} & Number of velocity truncations\\
\hline
 & {\em cs} & Vertical viscosity control structure \\
\hline
\end{DoxyParams}


Definition at line 1580 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
1580   \textcolor{keywordtype}{type}(ocean\_internal\_state), &
1581                    \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)}    :: MIS\textcolor{comment}{    !< The "MOM Internal State", a set of pointers}
1582 \textcolor{comment}{                                                   !! to the fields and accelerations that make}
1583 \textcolor{comment}{                                                   !! up the ocean's physical state}
1584   \textcolor{keywordtype}{type}(time\_type), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(in)}    :: Time\textcolor{comment}{   !< Current model time}
1585   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: G\textcolor{comment}{      !< Ocean grid structure}
1586   \textcolor{keywordtype}{type}(verticalGrid\_type), \textcolor{keywordtype}{intent(in)}    :: GV\textcolor{comment}{     !< Ocean vertical grid structure}
1587   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: US\textcolor{comment}{     !< A dimensional unit scaling type}
1588   \textcolor{keywordtype}{type}(param\_file\_type),   \textcolor{keywordtype}{intent(in)}    :: param\_file\textcolor{comment}{ !< File to parse for parameters}
1589   \textcolor{keywordtype}{type}(diag\_ctrl), \textcolor{keywordtype}{target}, \textcolor{keywordtype}{intent(inout)} :: diag\textcolor{comment}{   !< Diagnostic control structure}
1590   \textcolor{keywordtype}{type}(accel\_diag\_ptrs),   \textcolor{keywordtype}{intent(inout)} :: ADp\textcolor{comment}{    !< Acceleration diagnostic pointers}
1591   \textcolor{keywordtype}{type}(directories),       \textcolor{keywordtype}{intent(in)}    :: dirs\textcolor{comment}{   !< Relevant directory paths}
1592   \textcolor{keywordtype}{integer}, \textcolor{keywordtype}{target},         \textcolor{keywordtype}{intent(inout)} :: ntrunc\textcolor{comment}{ !< Number of velocity truncations}
1593   \textcolor{keywordtype}{type}(vertvisc\_CS),       \textcolor{keywordtype}{pointer}       :: CS\textcolor{comment}{     !< Vertical viscosity control structure}
1594 
1595   \textcolor{comment}{! Local variables}
1596 
1597   \textcolor{keywordtype}{real} :: hmix\_str\_dflt
1598   \textcolor{keywordtype}{real} :: Kv\_dflt \textcolor{comment}{! A default viscosity [m2 s-1].}
1599   \textcolor{keywordtype}{real} :: Hmix\_m  \textcolor{comment}{! A boundary layer thickness [m].}
1600   \textcolor{keywordtype}{logical} :: default\_2018\_answers
1601   \textcolor{keywordtype}{integer} :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
1602 \textcolor{comment}{! This include declares and sets the variable "version".}
1603 \textcolor{preprocessor}{#include "version\_variable.h"}
1604 \textcolor{preprocessor}{}  \textcolor{keywordtype}{character(len=40)}  :: mdl = \textcolor{stringliteral}{"MOM\_vert\_friction"} \textcolor{comment}{! This module's name.}
1605   \textcolor{keywordtype}{character(len=40)}  :: thickness\_units
1606 
1607   \textcolor{keywordflow}{if} (\textcolor{keyword}{associated}(cs)) \textcolor{keywordflow}{then}
1608     \textcolor{keyword}{call }mom\_error(warning, \textcolor{stringliteral}{"vertvisc\_init called with an associated "}// &
1609                             \textcolor{stringliteral}{"control structure."})
1610     \textcolor{keywordflow}{return}
1611 \textcolor{keywordflow}{  endif}
1612   \textcolor{keyword}{allocate}(cs)
1613 
1614   \textcolor{keywordflow}{if} (gv%Boussinesq) then; thickness\_units = \textcolor{stringliteral}{"m"}
1615   else; thickness\_units = \textcolor{stringliteral}{"kg m-2"};\textcolor{keywordflow}{ endif}
1616 
1617   isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1618   isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1619 
1620   cs%diag => diag ; cs%ntrunc => ntrunc ; ntrunc = 0
1621 
1622 \textcolor{comment}{! Default, read and log parameters}
1623   \textcolor{keyword}{call }log\_version(param\_file, mdl, version, \textcolor{stringliteral}{""}, log\_to\_all=.true., debugging=.true.)
1624   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DEFAULT\_2018\_ANSWERS"}, default\_2018\_answers, &
1625                  \textcolor{stringliteral}{"This sets the default value for the various \_2018\_ANSWERS parameters."}, &
1626                  default=.false.)
1627   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"VERT\_FRICTION\_2018\_ANSWERS"}, cs%answers\_2018, &
1628                  \textcolor{stringliteral}{"If true, use the order of arithmetic and expressions that recover the answers "}//&
1629                  \textcolor{stringliteral}{"from the end of 2018.  Otherwise, use expressions that do not use an arbitrary "}//&
1630                  \textcolor{stringliteral}{"hard-coded maximum viscous coupling coefficient between layers."}, &
1631                  default=default\_2018\_answers)
1632   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"BOTTOMDRAGLAW"}, cs%bottomdraglaw, &
1633                  \textcolor{stringliteral}{"If true, the bottom stress is calculated with a drag "}//&
1634                  \textcolor{stringliteral}{"law of the form c\_drag*|u|*u. The velocity magnitude "}//&
1635                  \textcolor{stringliteral}{"may be an assumed value or it may be based on the "}//&
1636                  \textcolor{stringliteral}{"actual velocity in the bottommost HBBL, depending on "}//&
1637                  \textcolor{stringliteral}{"LINEAR\_DRAG."}, default=.true.)
1638   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CHANNEL\_DRAG"}, cs%Channel\_drag, &
1639                  \textcolor{stringliteral}{"If true, the bottom drag is exerted directly on each "}//&
1640                  \textcolor{stringliteral}{"layer proportional to the fraction of the bottom it "}//&
1641                  \textcolor{stringliteral}{"overlies."}, default=.false.)
1642   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DIRECT\_STRESS"}, cs%direct\_stress, &
1643                  \textcolor{stringliteral}{"If true, the wind stress is distributed over the "}//&
1644                  \textcolor{stringliteral}{"topmost HMIX\_STRESS of fluid (like in HYCOM), and KVML "}//&
1645                  \textcolor{stringliteral}{"may be set to a very small value."}, default=.false.)
1646   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DYNAMIC\_VISCOUS\_ML"}, cs%dynamic\_viscous\_ML, &
1647                  \textcolor{stringliteral}{"If true, use a bulk Richardson number criterion to "}//&
1648                  \textcolor{stringliteral}{"determine the mixed layer thickness for viscosity."}, &
1649                  default=.false.)
1650   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"U\_TRUNC\_FILE"}, cs%u\_trunc\_file, &
1651                  \textcolor{stringliteral}{"The absolute path to a file into which the accelerations "}//&
1652                  \textcolor{stringliteral}{"leading to zonal velocity truncations are written. "}//&
1653                  \textcolor{stringliteral}{"Undefine this for efficiency if this diagnostic is not "}//&
1654                  \textcolor{stringliteral}{"needed."}, default=\textcolor{stringliteral}{" "}, debuggingparam=.true.)
1655   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"V\_TRUNC\_FILE"}, cs%v\_trunc\_file, &
1656                  \textcolor{stringliteral}{"The absolute path to a file into which the accelerations "}//&
1657                  \textcolor{stringliteral}{"leading to meridional velocity truncations are written. "}//&
1658                  \textcolor{stringliteral}{"Undefine this for efficiency if this diagnostic is not "}//&
1659                  \textcolor{stringliteral}{"needed."}, default=\textcolor{stringliteral}{" "}, debuggingparam=.true.)
1660   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HARMONIC\_VISC"}, cs%harmonic\_visc, &
1661                  \textcolor{stringliteral}{"If true, use the harmonic mean thicknesses for "}//&
1662                  \textcolor{stringliteral}{"calculating the vertical viscosity."}, default=.false.)
1663   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HARMONIC\_BL\_SCALE"}, cs%harm\_BL\_val, &
1664                  \textcolor{stringliteral}{"A scale to determine when water is in the boundary "}//&
1665                  \textcolor{stringliteral}{"layers based solely on harmonic mean thicknesses for "}//&
1666                  \textcolor{stringliteral}{"the purpose of determining the extent to which the "}//&
1667                  \textcolor{stringliteral}{"thicknesses used in the viscosities are upwinded."}, &
1668                  default=0.0, units=\textcolor{stringliteral}{"nondim"})
1669   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"DEBUG"}, cs%debug, default=.false.)
1670 
1671   \textcolor{keywordflow}{if} (gv%nkml < 1) &
1672     \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HMIX\_FIXED"}, cs%Hmix, &
1673                  \textcolor{stringliteral}{"The prescribed depth over which the near-surface "}//&
1674                  \textcolor{stringliteral}{"viscosity and diffusivity are elevated when the bulk "}//&
1675                  \textcolor{stringliteral}{"mixed layer is not used."}, units=\textcolor{stringliteral}{"m"}, scale=gv%m\_to\_H, &
1676                  unscaled=hmix\_m, fail\_if\_missing=.true.)
1677   \textcolor{keywordflow}{if} (cs%direct\_stress) \textcolor{keywordflow}{then}
1678     \textcolor{keywordflow}{if} (gv%nkml < 1) \textcolor{keywordflow}{then}
1679       \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HMIX\_STRESS"}, cs%Hmix\_stress, &
1680                  \textcolor{stringliteral}{"The depth over which the wind stress is applied if "}//&
1681                  \textcolor{stringliteral}{"DIRECT\_STRESS is true."}, units=\textcolor{stringliteral}{"m"}, default=hmix\_m, scale=gv%m\_to\_H)
1682     \textcolor{keywordflow}{else}
1683       \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HMIX\_STRESS"}, cs%Hmix\_stress, &
1684                  \textcolor{stringliteral}{"The depth over which the wind stress is applied if "}//&
1685                  \textcolor{stringliteral}{"DIRECT\_STRESS is true."}, units=\textcolor{stringliteral}{"m"}, fail\_if\_missing=.true., scale=gv%m\_to\_H)
1686 \textcolor{keywordflow}{    endif}
1687     \textcolor{keywordflow}{if} (cs%Hmix\_stress <= 0.0) \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"vertvisc\_init: "} // &
1688        \textcolor{stringliteral}{"HMIX\_STRESS must be set to a positive value if DIRECT\_STRESS is true."})
1689 \textcolor{keywordflow}{  endif}
1690   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KV"}, cs%Kv, &
1691                  \textcolor{stringliteral}{"The background kinematic viscosity in the interior. "}//&
1692                  \textcolor{stringliteral}{"The molecular value, ~1e-6 m2 s-1, may be used."}, &
1693                  units=\textcolor{stringliteral}{"m2 s-1"}, fail\_if\_missing=.true., scale=us%m2\_s\_to\_Z2\_T, unscaled=kv\_dflt)
1694 
1695   \textcolor{keywordflow}{if} (gv%nkml < 1) \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KVML"}, cs%Kvml, &
1696                  \textcolor{stringliteral}{"The kinematic viscosity in the mixed layer.  A typical "}//&
1697                  \textcolor{stringliteral}{"value is ~1e-2 m2 s-1. KVML is not used if "}//&
1698                  \textcolor{stringliteral}{"BULKMIXEDLAYER is true.  The default is set by KV."}, &
1699                  units=\textcolor{stringliteral}{"m2 s-1"}, default=kv\_dflt, scale=us%m2\_s\_to\_Z2\_T)
1700   \textcolor{keywordflow}{if} (.not.cs%bottomdraglaw) \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"KVBBL"}, cs%Kvbbl, &
1701                  \textcolor{stringliteral}{"The kinematic viscosity in the benthic boundary layer. "}//&
1702                  \textcolor{stringliteral}{"A typical value is ~1e-2 m2 s-1. KVBBL is not used if "}//&
1703                  \textcolor{stringliteral}{"BOTTOMDRAGLAW is true.  The default is set by KV."}, &
1704                  units=\textcolor{stringliteral}{"m2 s-1"}, default=kv\_dflt, scale=us%m2\_s\_to\_Z2\_T)
1705   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"HBBL"}, cs%Hbbl, &
1706                  \textcolor{stringliteral}{"The thickness of a bottom boundary layer with a "}//&
1707                  \textcolor{stringliteral}{"viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "}//&
1708                  \textcolor{stringliteral}{"the thickness over which near-bottom velocities are "}//&
1709                  \textcolor{stringliteral}{"averaged for the drag law if BOTTOMDRAGLAW is defined "}//&
1710                  \textcolor{stringliteral}{"but LINEAR\_DRAG is not."}, units=\textcolor{stringliteral}{"m"}, fail\_if\_missing=.true., scale=gv%m\_to\_H)
1711   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"MAXVEL"}, cs%maxvel, &
1712                  \textcolor{stringliteral}{"The maximum velocity allowed before the velocity "}//&
1713                  \textcolor{stringliteral}{"components are truncated."}, units=\textcolor{stringliteral}{"m s-1"}, default=3.0e8, scale=us%m\_s\_to\_L\_T)
1714   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CFL\_BASED\_TRUNCATIONS"}, cs%CFL\_based\_trunc, &
1715                  \textcolor{stringliteral}{"If true, base truncations on the CFL number, and not an "}//&
1716                  \textcolor{stringliteral}{"absolute speed."}, default=.true.)
1717   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CFL\_TRUNCATE"}, cs%CFL\_trunc, &
1718                  \textcolor{stringliteral}{"The value of the CFL number that will cause velocity "}//&
1719                  \textcolor{stringliteral}{"components to be truncated; instability can occur past 0.5."}, &
1720                  units=\textcolor{stringliteral}{"nondim"}, default=0.5)
1721   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CFL\_REPORT"}, cs%CFL\_report, &
1722                  \textcolor{stringliteral}{"The value of the CFL number that causes accelerations "}//&
1723                  \textcolor{stringliteral}{"to be reported; the default is CFL\_TRUNCATE."}, &
1724                  units=\textcolor{stringliteral}{"nondim"}, default=cs%CFL\_trunc)
1725   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CFL\_TRUNCATE\_RAMP\_TIME"}, cs%truncRampTime, &
1726                  \textcolor{stringliteral}{"The time over which the CFL truncation value is ramped "}//&
1727                  \textcolor{stringliteral}{"up at the beginning of the run."}, &
1728                  units=\textcolor{stringliteral}{"s"}, default=0.)
1729   cs%CFL\_truncE = cs%CFL\_trunc
1730   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"CFL\_TRUNCATE\_START"}, cs%CFL\_truncS, &
1731                  \textcolor{stringliteral}{"The start value of the truncation CFL number used when "}//&
1732                  \textcolor{stringliteral}{"ramping up CFL\_TRUNC."}, &
1733                  units=\textcolor{stringliteral}{"nondim"}, default=0.)
1734   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"STOKES\_MIXING\_COMBINED"}, cs%StokesMixing, &
1735                  \textcolor{stringliteral}{"Flag to use Stokes drift Mixing via the Lagrangian "}//&
1736                  \textcolor{stringliteral}{" current (Eulerian plus Stokes drift). "}//&
1737                  \textcolor{stringliteral}{" Still needs work and testing, so not recommended for use."},&
1738                  default=.false.)
1739   \textcolor{comment}{!BGR 04/04/2018\{}
1740   \textcolor{comment}{! StokesMixing is required for MOM6 for some Langmuir mixing parameterization.}
1741   \textcolor{comment}{!   The code used here has not been developed for vanishing layers or in}
1742   \textcolor{comment}{!   conjunction with any bottom friction.  Therefore, the following line is}
1743   \textcolor{comment}{!   added so this functionality cannot be used without user intervention in}
1744   \textcolor{comment}{!   the code.  This will prevent general use of this functionality until proper}
1745   \textcolor{comment}{!   care is given to the previously mentioned issues.  Comment out the following}
1746   \textcolor{comment}{!   MOM\_error to use, but do so at your own risk and with these points in mind.}
1747   \textcolor{comment}{!\}}
1748   \textcolor{keywordflow}{if} (cs%StokesMixing) \textcolor{keywordflow}{then}
1749     \textcolor{keyword}{call }mom\_error(fatal, \textcolor{stringliteral}{"Stokes mixing requires user intervention in the code.\(\backslash\)n"}//&
1750                           \textcolor{stringliteral}{"  Model now exiting.  See MOM\_vert\_friction.F90 for \(\backslash\)n"}//&
1751                           \textcolor{stringliteral}{"  details (search 'BGR 04/04/2018' to locate comment)."})
1752 \textcolor{keywordflow}{  endif}
1753   \textcolor{keyword}{call }get\_param(param\_file, mdl, \textcolor{stringliteral}{"VEL\_UNDERFLOW"}, cs%vel\_underflow, &
1754                  \textcolor{stringliteral}{"A negligibly small velocity magnitude below which velocity "}//&
1755                  \textcolor{stringliteral}{"components are set to 0.  A reasonable value might be "}//&
1756                  \textcolor{stringliteral}{"1e-30 m/s, which is less than an Angstrom divided by "}//&
1757                  \textcolor{stringliteral}{"the age of the universe."}, units=\textcolor{stringliteral}{"m s-1"}, default=0.0, scale=us%m\_s\_to\_L\_T)
1758 
1759   alloc\_(cs%a\_u(isdb:iedb,jsd:jed,nz+1)) ; cs%a\_u(:,:,:) = 0.0
1760   alloc\_(cs%h\_u(isdb:iedb,jsd:jed,nz))   ; cs%h\_u(:,:,:) = 0.0
1761   alloc\_(cs%a\_v(isd:ied,jsdb:jedb,nz+1)) ; cs%a\_v(:,:,:) = 0.0
1762   alloc\_(cs%h\_v(isd:ied,jsdb:jedb,nz))   ; cs%h\_v(:,:,:) = 0.0
1763 
1764   cs%id\_Kv\_slow = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Kv\_slow'}, diag%axesTi, time, &
1765      \textcolor{stringliteral}{'Slow varying vertical viscosity'}, \textcolor{stringliteral}{'m2 s-1'}, conversion=us%Z2\_T\_to\_m2\_s)
1766 
1767   cs%id\_Kv\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Kv\_u'}, diag%axesCuL, time, &
1768      \textcolor{stringliteral}{'Total vertical viscosity at u-points'}, \textcolor{stringliteral}{'m2 s-1'}, conversion=us%Z2\_T\_to\_m2\_s)
1769 
1770   cs%id\_Kv\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Kv\_v'}, diag%axesCvL, time, &
1771      \textcolor{stringliteral}{'Total vertical viscosity at v-points'}, \textcolor{stringliteral}{'m2 s-1'}, conversion=us%Z2\_T\_to\_m2\_s)
1772 
1773   cs%id\_au\_vv = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'au\_visc'}, diag%axesCui, time, &
1774      \textcolor{stringliteral}{'Zonal Viscous Vertical Coupling Coefficient'}, \textcolor{stringliteral}{'m s-1'}, conversion=us%Z\_to\_m*us%s\_to\_T)
1775 
1776   cs%id\_av\_vv = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'av\_visc'}, diag%axesCvi, time, &
1777      \textcolor{stringliteral}{'Meridional Viscous Vertical Coupling Coefficient'}, \textcolor{stringliteral}{'m s-1'}, conversion=us%Z\_to\_m*us%s\_to\_T)
1778 
1779   cs%id\_h\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Hu\_visc'}, diag%axesCuL, time, &
1780      \textcolor{stringliteral}{'Thickness at Zonal Velocity Points for Viscosity'}, thickness\_units, &
1781      conversion=gv%H\_to\_m)
1782 
1783   cs%id\_h\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'Hv\_visc'}, diag%axesCvL, time, &
1784      \textcolor{stringliteral}{'Thickness at Meridional Velocity Points for Viscosity'}, thickness\_units, &
1785      conversion=gv%H\_to\_m)
1786 
1787   cs%id\_hML\_u = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'HMLu\_visc'}, diag%axesCu1, time, &
1788      \textcolor{stringliteral}{'Mixed Layer Thickness at Zonal Velocity Points for Viscosity'}, thickness\_units, &
1789      conversion=gv%H\_to\_m)
1790 
1791   cs%id\_hML\_v = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'HMLv\_visc'}, diag%axesCv1, time, &
1792      \textcolor{stringliteral}{'Mixed Layer Thickness at Meridional Velocity Points for Viscosity'}, thickness\_units, &
1793      conversion=gv%H\_to\_m)
1794 
1795   cs%id\_du\_dt\_visc = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'du\_dt\_visc'}, diag%axesCuL, &
1796      time, \textcolor{stringliteral}{'Zonal Acceleration from Vertical Viscosity'}, \textcolor{stringliteral}{'m s-2'}, conversion=us%L\_T2\_to\_m\_s2)
1797   \textcolor{keywordflow}{if} (cs%id\_du\_dt\_visc > 0) \textcolor{keyword}{call }safe\_alloc\_ptr(adp%du\_dt\_visc,isdb,iedb,jsd,jed,nz)
1798   cs%id\_dv\_dt\_visc = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'dv\_dt\_visc'}, diag%axesCvL, &
1799      time, \textcolor{stringliteral}{'Meridional Acceleration from Vertical Viscosity'}, \textcolor{stringliteral}{'m s-2'}, conversion=us%L\_T2\_to\_m\_s2)
1800   \textcolor{keywordflow}{if} (cs%id\_dv\_dt\_visc > 0) \textcolor{keyword}{call }safe\_alloc\_ptr(adp%dv\_dt\_visc,isd,ied,jsdb,jedb,nz)
1801 
1802   cs%id\_taux\_bot = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'taux\_bot'}, diag%axesCu1, &
1803      time, \textcolor{stringliteral}{'Zonal Bottom Stress from Ocean to Earth'}, \textcolor{stringliteral}{'Pa'}, &
1804      conversion=us%RZ\_to\_kg\_m2*us%L\_T2\_to\_m\_s2)
1805   cs%id\_tauy\_bot = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'tauy\_bot'}, diag%axesCv1, &
1806      time, \textcolor{stringliteral}{'Meridional Bottom Stress from Ocean to Earth'}, \textcolor{stringliteral}{'Pa'}, &
1807      conversion=us%RZ\_to\_kg\_m2*us%L\_T2\_to\_m\_s2)
1808 
1809   \textcolor{comment}{!CS%id\_hf\_du\_dt\_visc = register\_diag\_field('ocean\_model', 'hf\_du\_dt\_visc', diag%axesCuL, Time, &}
1810   \textcolor{comment}{!    'Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', 'm s-2', &}
1811   \textcolor{comment}{!    v\_extensive=.true., conversion=US%L\_T2\_to\_m\_s2)}
1812   \textcolor{comment}{!if (CS%id\_hf\_du\_dt\_visc > 0) then}
1813   \textcolor{comment}{!  call safe\_alloc\_ptr(CS%hf\_du\_dt\_visc,IsdB,IedB,jsd,jed,nz)}
1814   \textcolor{comment}{!  call safe\_alloc\_ptr(ADp%du\_dt\_visc,IsdB,IedB,jsd,jed,nz)}
1815   \textcolor{comment}{!  call safe\_alloc\_ptr(ADp%diag\_hfrac\_u,IsdB,IedB,jsd,jed,nz)}
1816   \textcolor{comment}{!endif}
1817 
1818   \textcolor{comment}{!CS%id\_hf\_dv\_dt\_visc = register\_diag\_field('ocean\_model', 'hf\_dv\_dt\_visc', diag%axesCvL, Time, &}
1819   \textcolor{comment}{!    'Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', 'm s-2', &}
1820   \textcolor{comment}{!    v\_extensive=.true., conversion=US%L\_T2\_to\_m\_s2)}
1821   \textcolor{comment}{!if (CS%id\_hf\_dv\_dt\_visc > 0) then}
1822   \textcolor{comment}{!  call safe\_alloc\_ptr(CS%hf\_dv\_dt\_visc,isd,ied,JsdB,JedB,nz)}
1823   \textcolor{comment}{!  call safe\_alloc\_ptr(ADp%dv\_dt\_visc,isd,ied,JsdB,JedB,nz)}
1824   \textcolor{comment}{!  call safe\_alloc\_ptr(ADp%diag\_hfrac\_v,isd,ied,Jsd,JedB,nz)}
1825   \textcolor{comment}{!endif}
1826 
1827   cs%id\_hf\_du\_dt\_visc\_2d = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'hf\_du\_dt\_visc\_2d'}, diag%axesCu1, time, &
1828       \textcolor{stringliteral}{'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity'}, \textcolor{stringliteral}{'m s-2'}, &
1829       conversion=us%L\_T2\_to\_m\_s2)
1830   \textcolor{keywordflow}{if} (cs%id\_hf\_du\_dt\_visc\_2d > 0) \textcolor{keywordflow}{then}
1831     \textcolor{keyword}{call }safe\_alloc\_ptr(adp%du\_dt\_visc,isdb,iedb,jsd,jed,nz)
1832     \textcolor{keyword}{call }safe\_alloc\_ptr(adp%diag\_hfrac\_u,isdb,iedb,jsd,jed,nz)
1833 \textcolor{keywordflow}{  endif}
1834 
1835   cs%id\_hf\_dv\_dt\_visc\_2d = register\_diag\_field(\textcolor{stringliteral}{'ocean\_model'}, \textcolor{stringliteral}{'hf\_dv\_dt\_visc\_2d'}, diag%axesCv1, time, &
1836       \textcolor{stringliteral}{'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity'}, \textcolor{stringliteral}{'m s-2'}, &
1837       conversion=us%L\_T2\_to\_m\_s2)
1838   \textcolor{keywordflow}{if} (cs%id\_hf\_dv\_dt\_visc\_2d > 0) \textcolor{keywordflow}{then}
1839     \textcolor{keyword}{call }safe\_alloc\_ptr(adp%dv\_dt\_visc,isd,ied,jsdb,jedb,nz)
1840     \textcolor{keyword}{call }safe\_alloc\_ptr(adp%diag\_hfrac\_v,isd,ied,jsd,jedb,nz)
1841 \textcolor{keywordflow}{  endif}
1842 
1843   \textcolor{keywordflow}{if} ((len\_trim(cs%u\_trunc\_file) > 0) .or. (len\_trim(cs%v\_trunc\_file) > 0)) &
1844     \textcolor{keyword}{call }pointaccel\_init(mis, time, g, param\_file, diag, dirs, cs%PointAccel\_CSp)
1845 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__vert__friction_ac14bc439f3b7ae3000fa9883c95ebfe2}\label{namespacemom__vert__friction_ac14bc439f3b7ae3000fa9883c95ebfe2}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!vertvisc\+\_\+limit\+\_\+vel@{vertvisc\+\_\+limit\+\_\+vel}}
\index{vertvisc\+\_\+limit\+\_\+vel@{vertvisc\+\_\+limit\+\_\+vel}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{vertvisc\+\_\+limit\+\_\+vel()}{vertvisc\_limit\_vel()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+vert\+\_\+friction\+::vertvisc\+\_\+limit\+\_\+vel (\begin{DoxyParamCaption}\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(gv)), intent(inout)}]{u,  }\item[{real, dimension(szi\+\_\+(g),szjb\+\_\+(g),szk\+\_\+(gv)), intent(inout)}]{v,  }\item[{real, dimension(szi\+\_\+(g),szj\+\_\+(g),szk\+\_\+(gv)), intent(in)}]{h,  }\item[{type(accel\+\_\+diag\+\_\+ptrs), intent(in)}]{A\+Dp,  }\item[{type(cont\+\_\+diag\+\_\+ptrs), intent(in)}]{C\+Dp,  }\item[{type(mech\+\_\+forcing), intent(in)}]{forces,  }\item[{type(vertvisc\+\_\+type), intent(in)}]{visc,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS }\end{DoxyParamCaption})}



Velocity components which exceed a threshold for physically reasonable values are truncated. Optionally, any column with excessive velocities may be sent to a diagnostic reporting subroutine. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
\mbox{\tt in,out}  & {\em u} & Zonal velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em v} & Meridional velocity \mbox{[}L T-\/1 $\sim$$>$ m s-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h} & Layer thickness \mbox{[}H $\sim$$>$ m or kg m-\/2\mbox{]}\\
\hline
\mbox{\tt in}  & {\em adp} & Acceleration diagnostic pointers\\
\hline
\mbox{\tt in}  & {\em cdp} & Continuity diagnostic pointers\\
\hline
\mbox{\tt in}  & {\em forces} & A structure with the driving mechanical forces\\
\hline
\mbox{\tt in}  & {\em visc} & Viscosities and bottom drag\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
 & {\em cs} & Vertical viscosity control structure \\
\hline
\end{DoxyParams}


Definition at line 1369 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
1369   \textcolor{keywordtype}{type}(ocean\_grid\_type),   \textcolor{keywordtype}{intent(in)}    :: G\textcolor{comment}{      !< Ocean grid structure}
1370   \textcolor{keywordtype}{type}(verticalGrid\_type), \textcolor{keywordtype}{intent(in)}    :: GV\textcolor{comment}{     !< Ocean vertical grid structure}
1371   \textcolor{keywordtype}{type}(unit\_scale\_type),   \textcolor{keywordtype}{intent(in)}    :: US\textcolor{comment}{     !< A dimensional unit scaling type}
1372   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(GV))}, &
1373                            \textcolor{keywordtype}{intent(inout)} :: u\textcolor{comment}{      !< Zonal velocity [L T-1 ~> m s-1]}
1374   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(GV))}, &
1375                            \textcolor{keywordtype}{intent(inout)} :: v\textcolor{comment}{      !< Meridional velocity [L T-1 ~> m s-1]}
1376   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJ\_(G),SZK\_(GV))}, &
1377                            \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{      !< Layer thickness [H ~> m or kg m-2]}
1378   \textcolor{keywordtype}{type}(accel\_diag\_ptrs),   \textcolor{keywordtype}{intent(in)}    :: ADp\textcolor{comment}{    !< Acceleration diagnostic pointers}
1379   \textcolor{keywordtype}{type}(cont\_diag\_ptrs),    \textcolor{keywordtype}{intent(in)}    :: CDp\textcolor{comment}{    !< Continuity diagnostic pointers}
1380   \textcolor{keywordtype}{type}(mech\_forcing),      \textcolor{keywordtype}{intent(in)}    :: forces\textcolor{comment}{ !< A structure with the driving mechanical forces}
1381   \textcolor{keywordtype}{type}(vertvisc\_type),     \textcolor{keywordtype}{intent(in)}    :: visc\textcolor{comment}{   !< Viscosities and bottom drag}
1382   \textcolor{keywordtype}{real},                    \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{     !< Time increment [T ~> s]}
1383   \textcolor{keywordtype}{type}(vertvisc\_CS),       \textcolor{keywordtype}{pointer}       :: CS\textcolor{comment}{     !< Vertical viscosity control structure}
1384 
1385   \textcolor{comment}{! Local variables}
1386 
1387   \textcolor{keywordtype}{real} :: maxvel           \textcolor{comment}{! Velocities components greater than maxvel}
1388   \textcolor{keywordtype}{real} :: truncvel         \textcolor{comment}{! are truncated to truncvel, both [L T-1 ~> m s-1].}
1389   \textcolor{keywordtype}{real} :: CFL              \textcolor{comment}{! The local CFL number.}
1390   \textcolor{keywordtype}{real} :: H\_report         \textcolor{comment}{! A thickness below which not to report truncations.}
1391   \textcolor{keywordtype}{real} :: dt\_Rho0          \textcolor{comment}{! The timestep divided by the Boussinesq density [m2 T2 s-1 L-1 Z-1 R-1 ~> s m3
       kg-1].}
1392   \textcolor{keywordtype}{real} :: vel\_report(SZIB\_(G),SZJB\_(G))   \textcolor{comment}{! The velocity to report [L T-1 ~> m s-1]}
1393   \textcolor{keywordtype}{real} :: u\_old(SZIB\_(G),SZJ\_(G),SZK\_(G)) \textcolor{comment}{! The previous u-velocity [L T-1 ~> m s-1]}
1394   \textcolor{keywordtype}{real} :: v\_old(SZI\_(G),SZJB\_(G),SZK\_(G)) \textcolor{comment}{! The previous v-velocity [L T-1 ~> m s-1]}
1395   \textcolor{keywordtype}{logical} :: trunc\_any, dowrite(SZIB\_(G),SZJB\_(G))
1396   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1397   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1398   isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1399 
1400   maxvel = cs%maxvel
1401   truncvel = 0.9*maxvel
1402   h\_report = 6.0 * gv%Angstrom\_H
1403   dt\_rho0 = (us%L\_T\_to\_m\_s*us%Z\_to\_m) * dt / gv%Rho0
1404 
1405   \textcolor{keywordflow}{if} (len\_trim(cs%u\_trunc\_file) > 0) \textcolor{keywordflow}{then}
1406     \textcolor{comment}{!$OMP parallel do default(shared) private(trunc\_any,CFL)}
1407     \textcolor{keywordflow}{do} j=js,je
1408       trunc\_any = .false.
1409       \textcolor{keywordflow}{do} i=isq,ieq ; dowrite(i,j) = .false. ;\textcolor{keywordflow}{ enddo}
1410       \textcolor{keywordflow}{if} (cs%CFL\_based\_trunc) \textcolor{keywordflow}{then}
1411         \textcolor{keywordflow}{do} i=isq,ieq ; vel\_report(i,j) = 3.0e8*us%m\_s\_to\_L\_T ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! Speed of light default.}
1412         \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
1413           \textcolor{keywordflow}{if} (abs(u(i,j,k)) < cs%vel\_underflow) u(i,j,k) = 0.0
1414           \textcolor{keywordflow}{if} (u(i,j,k) < 0.0) \textcolor{keywordflow}{then}
1415             cfl = (-u(i,j,k) * dt) * (g%dy\_Cu(i,j) * g%IareaT(i+1,j))
1416           \textcolor{keywordflow}{else}
1417             cfl = (u(i,j,k) * dt) * (g%dy\_Cu(i,j) * g%IareaT(i,j))
1418 \textcolor{keywordflow}{          endif}
1419           \textcolor{keywordflow}{if} (cfl > cs%CFL\_trunc) trunc\_any = .true.
1420           \textcolor{keywordflow}{if} (cfl > cs%CFL\_report) \textcolor{keywordflow}{then}
1421             dowrite(i,j) = .true.
1422             vel\_report(i,j) = min(vel\_report(i,j), abs(u(i,j,k)))
1423 \textcolor{keywordflow}{          endif}
1424 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1425       \textcolor{keywordflow}{else}
1426         \textcolor{keywordflow}{do} i=isq,ieq; vel\_report(i,j) = maxvel;\textcolor{keywordflow}{ enddo}
1427         \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
1428           \textcolor{keywordflow}{if} (abs(u(i,j,k)) < cs%vel\_underflow) \textcolor{keywordflow}{then} ; u(i,j,k) = 0.0
1429           \textcolor{keywordflow}{elseif} (abs(u(i,j,k)) > maxvel) \textcolor{keywordflow}{then}
1430             dowrite(i,j) = .true. ; trunc\_any = .true.
1431 \textcolor{keywordflow}{          endif}
1432 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1433 \textcolor{keywordflow}{      endif}
1434 
1435       \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (dowrite(i,j)) \textcolor{keywordflow}{then}
1436         u\_old(i,j,:) = u(i,j,:)
1437 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
1438 
1439       \textcolor{keywordflow}{if} (trunc\_any) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (cs%CFL\_based\_trunc) \textcolor{keywordflow}{then}
1440         \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
1441           \textcolor{keywordflow}{if} ((u(i,j,k) * (dt * g%dy\_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL\_trunc) \textcolor{keywordflow}{then}
1442             u(i,j,k) = (-0.9*cs%CFL\_trunc) * (g%areaT(i+1,j) / (dt * g%dy\_Cu(i,j)))
1443             \textcolor{keywordflow}{if} (h(i,j,k) + h(i+1,j,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1444           \textcolor{keywordflow}{elseif} ((u(i,j,k) * (dt * g%dy\_Cu(i,j))) * g%IareaT(i,j) > cs%CFL\_trunc) \textcolor{keywordflow}{then}
1445             u(i,j,k) = (0.9*cs%CFL\_trunc) * (g%areaT(i,j) / (dt * g%dy\_Cu(i,j)))
1446             \textcolor{keywordflow}{if} (h(i,j,k) + h(i+1,j,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1447 \textcolor{keywordflow}{          endif}
1448 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1449       \textcolor{keywordflow}{else}
1450         \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (abs(u(i,j,k)) > maxvel) \textcolor{keywordflow}{then}
1451           u(i,j,k) = sign(truncvel,u(i,j,k))
1452           \textcolor{keywordflow}{if} (h(i,j,k) + h(i+1,j,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1453 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1454 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ endif}
1455 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! j-loop}
1456   \textcolor{keywordflow}{else}  \textcolor{comment}{! Do not report accelerations leading to large velocities.}
1457     \textcolor{keywordflow}{if} (cs%CFL\_based\_trunc) \textcolor{keywordflow}{then}
1458 \textcolor{comment}{!$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,dt,G,CS,h,H\_report)}
1459       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=isq,ieq
1460         \textcolor{keywordflow}{if} (abs(u(i,j,k)) < cs%vel\_underflow) \textcolor{keywordflow}{then} ; u(i,j,k) = 0.0
1461         \textcolor{keywordflow}{elseif} ((u(i,j,k) * (dt * g%dy\_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL\_trunc) \textcolor{keywordflow}{then}
1462           u(i,j,k) = (-0.9*cs%CFL\_trunc) * (g%areaT(i+1,j) / (dt * g%dy\_Cu(i,j)))
1463           \textcolor{keywordflow}{if} (h(i,j,k) + h(i+1,j,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1464         \textcolor{keywordflow}{elseif} ((u(i,j,k) * (dt * g%dy\_Cu(i,j))) * g%IareaT(i,j) > cs%CFL\_trunc) \textcolor{keywordflow}{then}
1465           u(i,j,k) = (0.9*cs%CFL\_trunc) * (g%areaT(i,j) / (dt * g%dy\_Cu(i,j)))
1466           \textcolor{keywordflow}{if} (h(i,j,k) + h(i+1,j,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1467 \textcolor{keywordflow}{        endif}
1468 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1469     \textcolor{keywordflow}{else}
1470 \textcolor{comment}{!$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,G,CS,truncvel,maxvel,h,H\_report)}
1471       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=js,je ; \textcolor{keywordflow}{do} i=isq,ieq
1472         \textcolor{keywordflow}{if} (abs(u(i,j,k)) < cs%vel\_underflow) \textcolor{keywordflow}{then} ; u(i,j,k) = 0.0
1473         \textcolor{keywordflow}{elseif} (abs(u(i,j,k)) > maxvel) \textcolor{keywordflow}{then}
1474           u(i,j,k) = sign(truncvel, u(i,j,k))
1475           \textcolor{keywordflow}{if} (h(i,j,k) + h(i+1,j,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1476 \textcolor{keywordflow}{        endif}
1477 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1478 \textcolor{keywordflow}{    endif}
1479 \textcolor{keywordflow}{  endif}
1480 
1481   \textcolor{keywordflow}{if} (len\_trim(cs%u\_trunc\_file) > 0) \textcolor{keywordflow}{then}
1482     \textcolor{keywordflow}{do} j=js,je; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (dowrite(i,j)) \textcolor{keywordflow}{then}
1483 \textcolor{comment}{!   Here the diagnostic reporting subroutines are called if}
1484 \textcolor{comment}{! unphysically large values were found.}
1485       \textcolor{keyword}{call }write\_u\_accel(i, j, u\_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel\_CSp, &
1486                vel\_report(i,j), forces%taux(i,j)*dt\_rho0, a=cs%a\_u, hv=cs%h\_u)
1487 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1488 \textcolor{keywordflow}{  endif}
1489 
1490   \textcolor{keywordflow}{if} (len\_trim(cs%v\_trunc\_file) > 0) \textcolor{keywordflow}{then}
1491     \textcolor{comment}{!$OMP parallel do default(shared) private(trunc\_any,CFL)}
1492     \textcolor{keywordflow}{do} j=jsq,jeq
1493       trunc\_any = .false.
1494       \textcolor{keywordflow}{do} i=is,ie ; dowrite(i,j) = .false. ;\textcolor{keywordflow}{ enddo}
1495       \textcolor{keywordflow}{if} (cs%CFL\_based\_trunc) \textcolor{keywordflow}{then}
1496         \textcolor{keywordflow}{do} i=is,ie ; vel\_report(i,j) = 3.0e8*us%m\_s\_to\_L\_T ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! Speed of light default.}
1497         \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
1498           \textcolor{keywordflow}{if} (abs(v(i,j,k)) < cs%vel\_underflow) v(i,j,k) = 0.0
1499           \textcolor{keywordflow}{if} (v(i,j,k) < 0.0) \textcolor{keywordflow}{then}
1500             cfl = (-v(i,j,k) * dt) * (g%dx\_Cv(i,j) * g%IareaT(i,j+1))
1501           \textcolor{keywordflow}{else}
1502             cfl = (v(i,j,k) * dt) * (g%dx\_Cv(i,j) * g%IareaT(i,j))
1503 \textcolor{keywordflow}{          endif}
1504           \textcolor{keywordflow}{if} (cfl > cs%CFL\_trunc) trunc\_any = .true.
1505           \textcolor{keywordflow}{if} (cfl > cs%CFL\_report) \textcolor{keywordflow}{then}
1506             dowrite(i,j) = .true.
1507             vel\_report(i,j) = min(vel\_report(i,j), abs(v(i,j,k)))
1508 \textcolor{keywordflow}{          endif}
1509 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1510       \textcolor{keywordflow}{else}
1511         \textcolor{keywordflow}{do} i=is,ie ; vel\_report(i,j) = maxvel ;\textcolor{keywordflow}{ enddo}
1512         \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
1513           \textcolor{keywordflow}{if} (abs(v(i,j,k)) < cs%vel\_underflow) \textcolor{keywordflow}{then} ; v(i,j,k) = 0.0
1514           \textcolor{keywordflow}{elseif} (abs(v(i,j,k)) > maxvel) \textcolor{keywordflow}{then}
1515             dowrite(i,j) = .true. ; trunc\_any = .true.
1516 \textcolor{keywordflow}{          endif}
1517 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1518 \textcolor{keywordflow}{      endif}
1519 
1520       \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (dowrite(i,j)) \textcolor{keywordflow}{then}
1521         v\_old(i,j,:) = v(i,j,:)
1522 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ enddo}
1523 
1524       \textcolor{keywordflow}{if} (trunc\_any) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{if} (cs%CFL\_based\_trunc) \textcolor{keywordflow}{then}
1525         \textcolor{keywordflow}{do} k=1,nz; \textcolor{keywordflow}{do} i=is,ie
1526           \textcolor{keywordflow}{if} ((v(i,j,k) * (dt * g%dx\_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL\_trunc) \textcolor{keywordflow}{then}
1527             v(i,j,k) = (-0.9*cs%CFL\_trunc) * (g%areaT(i,j+1) / (dt * g%dx\_Cv(i,j)))
1528             \textcolor{keywordflow}{if} (h(i,j,k) + h(i,j+1,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1529           \textcolor{keywordflow}{elseif} ((v(i,j,k) * (dt * g%dx\_Cv(i,j))) * g%IareaT(i,j) > cs%CFL\_trunc) \textcolor{keywordflow}{then}
1530             v(i,j,k) = (0.9*cs%CFL\_trunc) * (g%areaT(i,j) / (dt * g%dx\_Cv(i,j)))
1531             \textcolor{keywordflow}{if} (h(i,j,k) + h(i,j+1,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1532 \textcolor{keywordflow}{          endif}
1533 \textcolor{keywordflow}{        enddo} ;\textcolor{keywordflow}{ enddo}
1534       \textcolor{keywordflow}{else}
1535         \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (abs(v(i,j,k)) > maxvel) \textcolor{keywordflow}{then}
1536           v(i,j,k) = sign(truncvel,v(i,j,k))
1537           \textcolor{keywordflow}{if} (h(i,j,k) + h(i,j+1,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1538 \textcolor{keywordflow}{        endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1539 \textcolor{keywordflow}{      endif} ;\textcolor{keywordflow}{ endif}
1540 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! J-loop}
1541   \textcolor{keywordflow}{else}  \textcolor{comment}{! Do not report accelerations leading to large velocities.}
1542     \textcolor{keywordflow}{if} (cs%CFL\_based\_trunc) \textcolor{keywordflow}{then}
1543       \textcolor{comment}{!$OMP parallel do default(shared)}
1544       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} i=is,ie
1545         \textcolor{keywordflow}{if} (abs(v(i,j,k)) < cs%vel\_underflow) \textcolor{keywordflow}{then} ; v(i,j,k) = 0.0
1546         \textcolor{keywordflow}{elseif} ((v(i,j,k) * (dt * g%dx\_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL\_trunc) \textcolor{keywordflow}{then}
1547           v(i,j,k) = (-0.9*cs%CFL\_trunc) * (g%areaT(i,j+1) / (dt * g%dx\_Cv(i,j)))
1548           \textcolor{keywordflow}{if} (h(i,j,k) + h(i,j+1,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1549         \textcolor{keywordflow}{elseif} ((v(i,j,k) * (dt * g%dx\_Cv(i,j))) * g%IareaT(i,j) > cs%CFL\_trunc) \textcolor{keywordflow}{then}
1550           v(i,j,k) = (0.9*cs%CFL\_trunc) * (g%areaT(i,j) / (dt * g%dx\_Cv(i,j)))
1551           \textcolor{keywordflow}{if} (h(i,j,k) + h(i,j+1,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1552 \textcolor{keywordflow}{        endif}
1553 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1554     \textcolor{keywordflow}{else}
1555       \textcolor{comment}{!$OMP parallel do default(shared)}
1556       \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} j=jsq,jeq ; \textcolor{keywordflow}{do} i=is,ie
1557         \textcolor{keywordflow}{if} (abs(v(i,j,k)) < cs%vel\_underflow) \textcolor{keywordflow}{then} ; v(i,j,k) = 0.0
1558         \textcolor{keywordflow}{elseif} (abs(v(i,j,k)) > maxvel) \textcolor{keywordflow}{then}
1559           v(i,j,k) = sign(truncvel, v(i,j,k))
1560           \textcolor{keywordflow}{if} (h(i,j,k) + h(i,j+1,k) > h\_report) cs%ntrunc = cs%ntrunc + 1
1561 \textcolor{keywordflow}{        endif}
1562 \textcolor{keywordflow}{      enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1563 \textcolor{keywordflow}{    endif}
1564 \textcolor{keywordflow}{  endif}
1565 
1566   \textcolor{keywordflow}{if} (len\_trim(cs%v\_trunc\_file) > 0) \textcolor{keywordflow}{then}
1567     \textcolor{keywordflow}{do} j=jsq,jeq; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (dowrite(i,j)) \textcolor{keywordflow}{then}
1568 \textcolor{comment}{!   Here the diagnostic reporting subroutines are called if}
1569 \textcolor{comment}{! unphysically large values were found.}
1570       \textcolor{keyword}{call }write\_v\_accel(i, j, v\_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel\_CSp, &
1571                vel\_report(i,j), forces%tauy(i,j)*dt\_rho0, a=cs%a\_v, hv=cs%h\_v)
1572 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
1573 \textcolor{keywordflow}{  endif}
1574 
\end{DoxyCode}
\mbox{\Hypertarget{namespacemom__vert__friction_ad0b7bdc6695ee0c4207797bc94672863}\label{namespacemom__vert__friction_ad0b7bdc6695ee0c4207797bc94672863}} 
\index{mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}!vertvisc\+\_\+remnant@{vertvisc\+\_\+remnant}}
\index{vertvisc\+\_\+remnant@{vertvisc\+\_\+remnant}!mom\+\_\+vert\+\_\+friction@{mom\+\_\+vert\+\_\+friction}}
\subsubsection{\texorpdfstring{vertvisc\+\_\+remnant()}{vertvisc\_remnant()}}
{\footnotesize\ttfamily subroutine, public mom\+\_\+vert\+\_\+friction\+::vertvisc\+\_\+remnant (\begin{DoxyParamCaption}\item[{type(vertvisc\+\_\+type), intent(in)}]{visc,  }\item[{real, dimension(szib\+\_\+(g),szj\+\_\+(g),szk\+\_\+(gv)), intent(inout)}]{visc\+\_\+rem\+\_\+u,  }\item[{real, dimension( g \%isd\+: g \%ied, g \%jsdb\+: g \%jedb, gv \%ke), intent(inout)}]{visc\+\_\+rem\+\_\+v,  }\item[{real, intent(in)}]{dt,  }\item[{type(ocean\+\_\+grid\+\_\+type), intent(in)}]{G,  }\item[{type(verticalgrid\+\_\+type), intent(in)}]{GV,  }\item[{type(unit\+\_\+scale\+\_\+type), intent(in)}]{US,  }\item[{type(\mbox{\hyperlink{structmom__vert__friction_1_1vertvisc__cs}{vertvisc\+\_\+cs}}), pointer}]{CS }\end{DoxyParamCaption})}



Calculate the fraction of momentum originally in a layer that remains after a time-\/step of viscosity, and the fraction of a time-\/step\textquotesingle{}s worth of barotropic acceleration that a layer experiences after viscosity is applied. 


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em g} & Ocean grid structure\\
\hline
\mbox{\tt in}  & {\em gv} & Ocean vertical grid structure\\
\hline
\mbox{\tt in}  & {\em visc} & Viscosities and bottom drag\\
\hline
\mbox{\tt in,out}  & {\em visc\+\_\+rem\+\_\+u} & Fraction of a time-\/step\textquotesingle{}s worth of a\\
\hline
\mbox{\tt in,out}  & {\em visc\+\_\+rem\+\_\+v} & Fraction of a time-\/step\textquotesingle{}s worth of a\\
\hline
\mbox{\tt in}  & {\em dt} & Time increment \mbox{[}T $\sim$$>$ s\mbox{]}\\
\hline
\mbox{\tt in}  & {\em us} & A dimensional unit scaling type\\
\hline
 & {\em cs} & Vertical viscosity control structure \\
\hline
\end{DoxyParams}


Definition at line 509 of file M\+O\+M\+\_\+vert\+\_\+friction.\+F90.


\begin{DoxyCode}
509   \textcolor{keywordtype}{type}(ocean\_grid\_type), \textcolor{keywordtype}{intent(in)}   :: G\textcolor{comment}{    !< Ocean grid structure}
510   \textcolor{keywordtype}{type}(verticalGrid\_type), \textcolor{keywordtype}{intent(in)} :: GV\textcolor{comment}{   !< Ocean vertical grid structure}
511   \textcolor{keywordtype}{type}(vertvisc\_type),   \textcolor{keywordtype}{intent(in)}   :: visc\textcolor{comment}{ !< Viscosities and bottom drag}
512   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZIB\_(G),SZJ\_(G),SZK\_(GV))}, &
513                          \textcolor{keywordtype}{intent(inout)} :: visc\_rem\_u\textcolor{comment}{ !< Fraction of a time-step's worth of a}
514 \textcolor{comment}{                                              !! barotopic acceleration that a layer experiences after}
515 \textcolor{comment}{                                              !! viscosity is applied in the zonal direction [nondim]}
516   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(SZI\_(G),SZJB\_(G),SZK\_(GV))}, &
517                          \textcolor{keywordtype}{intent(inout)} :: visc\_rem\_v\textcolor{comment}{ !< Fraction of a time-step's worth of a}
518 \textcolor{comment}{                                              !! barotopic acceleration that a layer experiences after}
519 \textcolor{comment}{                                              !! viscosity is applied in the meridional direction [nondim]}
520   \textcolor{keywordtype}{real},                  \textcolor{keywordtype}{intent(in)}    :: dt\textcolor{comment}{  !< Time increment [T ~> s]}
521   \textcolor{keywordtype}{type}(unit\_scale\_type), \textcolor{keywordtype}{intent(in)}    :: US\textcolor{comment}{  !< A dimensional unit scaling type}
522   \textcolor{keywordtype}{type}(vertvisc\_CS),     \textcolor{keywordtype}{pointer}       :: CS\textcolor{comment}{  !< Vertical viscosity control structure}
523 
524   \textcolor{comment}{! Local variables}
525 
526   \textcolor{keywordtype}{real} :: b1(SZIB\_(G))          \textcolor{comment}{! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].}
527   \textcolor{keywordtype}{real} :: c1(SZIB\_(G),SZK\_(G))  \textcolor{comment}{! A variable used by the tridiagonal solver [nondim].}
528   \textcolor{keywordtype}{real} :: d1(SZIB\_(G))          \textcolor{comment}{! d1=1-c1 is used by the tridiagonal solver [nondim].}
529   \textcolor{keywordtype}{real} :: Ray(SZIB\_(G),SZK\_(G)) \textcolor{comment}{! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].}
530   \textcolor{keywordtype}{real} :: b\_denom\_1   \textcolor{comment}{! The first term in the denominator of b1 [H ~> m or kg m-2].}
531   \textcolor{keywordtype}{real} :: dt\_Z\_to\_H        \textcolor{comment}{! The time step times the conversion from Z to the}
532                            \textcolor{comment}{! units of thickness [T H Z-1 ~> s or s kg m-3].}
533   \textcolor{keywordtype}{logical} :: do\_i(SZIB\_(G))
534 
535   \textcolor{keywordtype}{integer} :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz
536   is = g%isc ; ie = g%iec
537   isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
538 
539   \textcolor{keywordflow}{if} (.not.\textcolor{keyword}{associated}(cs)) \textcolor{keyword}{call }mom\_error(fatal,\textcolor{stringliteral}{"MOM\_vert\_friction(visc): "}// &
540          \textcolor{stringliteral}{"Module must be initialized before it is used."})
541 
542   dt\_z\_to\_h = dt*gv%Z\_to\_H
543 
544   \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; ray(i,k) = 0.0 ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
545 
546   \textcolor{comment}{! Find the zonal viscous using a modification of a standard tridagonal solver.}
547 \textcolor{comment}{!$OMP parallel do default(none) shared(G,Isq,Ieq,CS,nz,visc,dt\_Z\_to\_H,visc\_rem\_u) &}
548 \textcolor{comment}{!$OMP                     firstprivate(Ray)                                       &}
549 \textcolor{comment}{!$OMP                          private(do\_i,b\_denom\_1,b1,d1,c1)}
550   \textcolor{keywordflow}{do} j=g%jsc,g%jec
551     \textcolor{keywordflow}{do} i=isq,ieq ; do\_i(i) = (g%mask2dCu(i,j) > 0) ;\textcolor{keywordflow}{ enddo}
552 
553     \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=isq,ieq
554       ray(i,k) = visc%Ray\_u(i,j,k)
555 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
556 
557     \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
558       b\_denom\_1 = cs%h\_u(i,j,1) + dt\_z\_to\_h * (ray(i,1) + cs%a\_u(i,j,1))
559       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h*cs%a\_u(i,j,2))
560       d1(i) = b\_denom\_1 * b1(i)
561       visc\_rem\_u(i,j,1) = b1(i) * cs%h\_u(i,j,1)
562 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo}
563     \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
564       c1(i,k) = dt\_z\_to\_h * cs%a\_u(i,j,k)*b1(i)
565       b\_denom\_1 = cs%h\_u(i,j,k) + dt\_z\_to\_h * (ray(i,k) + cs%a\_u(i,j,k)*d1(i))
566       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h * cs%a\_u(i,j,k+1))
567       d1(i) = b\_denom\_1 * b1(i)
568       visc\_rem\_u(i,j,k) = (cs%h\_u(i,j,k) + dt\_z\_to\_h * cs%a\_u(i,j,k) * visc\_rem\_u(i,j,k-1)) * b1(i)
569 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
570     \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=isq,ieq ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
571       visc\_rem\_u(i,j,k) = visc\_rem\_u(i,j,k) + c1(i,k+1)*visc\_rem\_u(i,j,k+1)
572 
573 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i and k loops}
574 
575 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end u-component j loop}
576 
577   \textcolor{comment}{! Now find the meridional viscous using a modification.}
578 \textcolor{comment}{!$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,G,CS,visc,dt\_Z\_to\_H,visc\_rem\_v,nz) &}
579 \textcolor{comment}{!$OMP                     firstprivate(Ray)                                             &}
580 \textcolor{comment}{!$OMP                          private(do\_i,b\_denom\_1,b1,d1,c1)}
581   \textcolor{keywordflow}{do} j=jsq,jeq
582     \textcolor{keywordflow}{do} i=is,ie ; do\_i(i) = (g%mask2dCv(i,j) > 0) ;\textcolor{keywordflow}{ enddo}
583 
584     \textcolor{keywordflow}{if} (cs%Channel\_drag) \textcolor{keywordflow}{then} ; \textcolor{keywordflow}{do} k=1,nz ; \textcolor{keywordflow}{do} i=is,ie
585       ray(i,k) = visc%Ray\_v(i,j,k)
586 \textcolor{keywordflow}{    enddo} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ endif}
587 
588     \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
589       b\_denom\_1 = cs%h\_v(i,j,1) + dt\_z\_to\_h * (ray(i,1) + cs%a\_v(i,j,1))
590       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h*cs%a\_v(i,j,2))
591       d1(i) = b\_denom\_1 * b1(i)
592       visc\_rem\_v(i,j,1) = b1(i) * cs%h\_v(i,j,1)
593 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo}
594     \textcolor{keywordflow}{do} k=2,nz ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
595       c1(i,k) = dt\_z\_to\_h * cs%a\_v(i,j,k)*b1(i)
596       b\_denom\_1 = cs%h\_v(i,j,k) + dt\_z\_to\_h * (ray(i,k) + cs%a\_v(i,j,k)*d1(i))
597       b1(i) = 1.0 / (b\_denom\_1 + dt\_z\_to\_h * cs%a\_v(i,j,k+1))
598       d1(i) = b\_denom\_1 * b1(i)
599       visc\_rem\_v(i,j,k) = (cs%h\_v(i,j,k) + dt\_z\_to\_h * cs%a\_v(i,j,k) * visc\_rem\_v(i,j,k-1)) * b1(i)
600 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo}
601     \textcolor{keywordflow}{do} k=nz-1,1,-1 ; \textcolor{keywordflow}{do} i=is,ie ; \textcolor{keywordflow}{if} (do\_i(i)) \textcolor{keywordflow}{then}
602       visc\_rem\_v(i,j,k) = visc\_rem\_v(i,j,k) + c1(i,k+1)*visc\_rem\_v(i,j,k+1)
603 \textcolor{keywordflow}{    endif} ;\textcolor{keywordflow}{ enddo} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! i and k loops}
604 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end of v-component J loop}
605 
606   \textcolor{keywordflow}{if} (cs%debug) \textcolor{keywordflow}{then}
607     \textcolor{keyword}{call }uvchksum(\textcolor{stringliteral}{"visc\_rem\_[uv]"}, visc\_rem\_u, visc\_rem\_v, g%HI, haloshift=0, &
608                   scalar\_pair=.true.)
609 \textcolor{keywordflow}{  endif}
610 
\end{DoxyCode}
