\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}