--- manual/s_phys_pkgs/text/kpp.tex	2005/07/14 21:04:52	1.4
+++ manual/s_phys_pkgs/text/kpp.tex	2005/08/02 22:26:58	1.7
@@ -1,7 +1,522 @@
-\section{KPP Package: Ocean vertical mixing -- 
-the nonlocal K-profile parameterization scheme}
+\subsection{KPP: Nonlocal K-Profile Parameterization for 
+Vertical Mixing}
 
 \label{sec:pkg:kpp}
 \begin{rawhtml}
 <!-- CMIREDIR:package_kpp: -->
 \end{rawhtml}
+
+Authors: Dimitris Menemenlis and Patrick Heimbach
+
+\subsubsection{Introduction
+\label{sec:pkg:kpp:intro}}
+
+The nonlocal K-Profile Parameterization (KPP) scheme
+of \cite{lar-eta:94} unifies the treatment of a variety of 
+unresolved processes involved in vertical mixing. 
+To consider it as one mixing scheme is, in the view of the authors,
+somewhat misleading since it consists of several entities
+to deal with distinct mixing processes in the ocean's surface 
+boundary layer, and the interior:
+%
+\begin{enumerate}
+%
+\item
+mixing in the interior is goverened by
+shear instability (modeled as function of the local gradient
+Richardson number), internal wave activity (assumed constant),
+and double-diffusion (not implemented here).
+%
+\item
+a boundary layer depth $h$ or \texttt{hbl} is determined
+at each grid point, based on a critical value of turbulent
+processes parameterized by a bulk Richardson number;
+%
+\item
+mixing is strongly enhanced in the boundary layer under the
+stabilizing or destabilizing influence of surface forcing 
+(buoyancy and momentum) enabling boundary layer properties
+to penetrate well into the thermocline; 
+mixing is represented through a polynomial profile whose 
+coefficients are determined subject to several contraints;
+%
+\item
+the boundary-layer profile is made to agree with similarity
+theory of turbulence and is matched, in the asymptotic sense
+(function and derivative agree at the boundary),
+to the interior thus fixing the polynomial coefficients; 
+matching allows for some fraction of the boundary layer mixing
+to affect the interior, and vice versa;
+%
+\item
+a ``non-local'' term $\hat{\gamma}$ or \texttt{ghat}
+which is independent of the vertical property gradient further 
+enhances mixing where the water column is unstable
+%
+\end{enumerate}
+%
+The scheme has been extensively compared to observations
+(see e.g. \cite{lar-eta:97}) and is now coomon in many
+ocean models.
+
+The following sections will describe the KPP package
+configuration and compiling (\ref{sec:pkg:kpp:comp}),
+the settings and choices of runtime parameters
+(\ref{sec:pkg:kpp:runtime}),
+more detailed description of equations to which these
+parameters relate (\ref{sec:pkg:kpp:equations}),
+and key subroutines where they are used (\ref{sec:pkg:kpp:subroutines}),
+and diagnostics output of KPP-derived diffusivities, viscosities
+and boundary-layer/mixed-layer depths.
+
+%----------------------------------------------------------------------
+
+\subsubsection{KPP configuration and compiling
+\label{sec:pkg:kpp:comp}}
+
+As with all MITgcm packages, KPP can be turned on or off at compile time
+%
+\begin{itemize}
+%
+\item
+using the \texttt{packages.conf} file by adding \texttt{kpp} to it,
+%
+\item
+or using \texttt{genmake2} adding
+\texttt{-enable=kpp} or \texttt{-disable=kpp} switches
+%
+\end{itemize}
+(see Section \ref{sect:buildingCode}).
+
+Parts of the KPP code can be enabled or disabled at compile time
+via CPP preprocessor flags. These options are set in
+\texttt{KPP\_OPTIONS.h}. Table \ref{tab:pkg:kpp:cpp} summarizes them.
+
+\begin{table}[h!]
+\centering
+  \label{tab:pkg:kpp:cpp}
+  {\footnotesize
+    \begin{tabular}{|l|l|}
+      \hline 
+      \textbf{CPP option}  &  \textbf{Description}  \\
+      \hline \hline
+        \texttt{\_KPP\_RL} & 
+          ~ \\
+        \texttt{FRUGAL\_KPP} & 
+          ~ \\
+        \texttt{KPP\_SMOOTH\_SHSQ} & 
+          ~ \\
+        \texttt{KPP\_SMOOTH\_DVSQ} & 
+          ~ \\
+        \texttt{KPP\_SMOOTH\_DENS} & 
+          ~ \\
+        \texttt{KPP\_SMOOTH\_VISC} & 
+          ~ \\
+        \texttt{KPP\_SMOOTH\_DIFF} & 
+          ~ \\
+        \texttt{KPP\_ESTIMATE\_UREF} & 
+          ~ \\
+        \texttt{INCLUDE\_DIAGNOSTICS\_INTERFACE\_CODE} & 
+          ~ \\
+        \texttt{KPP\_GHAT} & 
+          ~ \\
+        \texttt{EXCLUDE\_KPP\_SHEAR\_MIX} & 
+          ~ \\
+      \hline
+    \end{tabular}
+  }
+  \caption{~}
+\end{table}
+
+
+%----------------------------------------------------------------------
+
+\subsubsection{Run-time parameters
+\label{sec:pkg:kpp:runtime}}
+
+Run-time parameters are set in files 
+\texttt{data.pkg} and \texttt{data.kpp}
+which are read in \texttt{kpp\_readparms.F}.
+Run-time parameters may be broken into 3 categories:
+(i) switching on/off the package at runtime,
+(ii) required MITgcm flags,
+(iii) package flags and parameters.
+
+\paragraph{Enabling the package}
+~ \\
+%
+The KPP package is switched on at runtime by setting
+\texttt{useKPP = .TRUE.} in \texttt{data.pkg}.
+
+\paragraph{Required MITgcm flags}
+~ \\
+%
+The following flags/parameters of the MITgcm dynamical
+kernel need to be set in conjunction with KPP:
+
+\begin{tabular}{ll}
+\texttt{implicitViscosity = .TRUE.} & enable implicit vertical viscosity \\
+\texttt{implicitDiffusion = .TRUE.} & enable implicit vertical diffusion \\
+\end{tabular}
+
+
+\paragraph{Package flags and parameters}
+~ \\
+%
+Table \ref{tab:pkg:kpp:runtime_flags} summarizes the
+runtime flags that are set in \texttt{data.pkg}, and
+their default values.
+
+\begin{table}[h!]
+\centering
+  \label{tab:pkg:kpp:runtime_flags}
+  {\footnotesize
+    \begin{tabular}{|l|c|l|}
+      \hline 
+      \textbf{Flag/parameter} & \textbf{default} &  \textbf{Description}  \\
+      \hline \hline
+         \multicolumn{3}{|c|}{\textit{I/O related parameters} } \\
+         \hline
+        kpp\_freq & \texttt{deltaTClock} & 
+           Recomputation frequency for KPP fields \\
+        kpp\_dumpFreq & \texttt{dumpFreq} & 
+           Dump frequency of KPP field snapshots \\
+        kpp\_taveFreq & \texttt{taveFreq} & 
+           Averaging and dump frequency of KPP fields \\
+        KPPmixingMaps & \texttt{.FALSE.} & 
+           include KPP diagnostic maps in STDOUT \\
+        KPPwriteState & \texttt{.FALSE.} & 
+           write KPP state to file \\
+        KPP\_ghatUseTotalDiffus & \texttt{.FALSE.} & 
+           if \texttt{.T.} compute non-local term using total vertical diffusivity \\
+        ~ & ~ &
+           if \texttt{.F.} use KPP vertical diffusivity \\
+      \hline
+      \multicolumn{3}{|c|}{\textit{Genral KPP parameters} } \\
+      \hline
+        minKPPhbl & \texttt{delRc(1)} & 
+           Minimum boundary layer depth \\
+        epsilon & 0.1 & 
+           nondimensional extent of the surface layer \\
+        vonk & 0.4 & 
+           von Karman constant \\
+        dB\_dz & 5.2E-5 1/s$^2$ & 
+           maximum dB/dz in mixed layer hMix \\
+        concs & 98.96 &
+           ~ \\
+        concv & 1.8 &
+           ~ \\
+      \hline
+      \multicolumn{3}{|c|}{\textit{Boundary layer parameters (S/R \texttt{bldepth})} } \\
+      \hline
+        Ricr & 0.3 & 
+           critical bulk Richardson number \\
+        cekman & 0.7 & 
+           coefficient for Ekman depth \\
+        cmonob & 1.0 & 
+           coefficient for Monin-Obukhov depth \\
+        concv & 1.8 & 
+           ratio of interior to  entrainment depth buoyancy frequency \\
+        hbf & 1.0 & 
+           fraction of depth to which absorbed solar radiation contributes \\
+        ~ & ~ &
+           to surface buoyancy forcing \\
+        Vtc & \texttt{~} & 
+           non-dim. coeff. for velocity scale of turbulant velocity shear \\
+        ~ & ~ &
+           ( = function of concv,concs,epsilon,vonk,Ricr) \\
+      \hline
+      \multicolumn{3}{|c|}{\textit{Boundary layer mixing parameters (S/R \texttt{blmix})} } \\
+      \hline
+        cstar & 10. & 
+           proportionality coefficient for nonlocal transport \\
+        cg & ~ & 
+           non-dimensional coefficient for counter-gradient term \\
+        ~ & ~ &
+           ( = function of cstar,vonk,concs,epsilon) \\
+      \hline
+      \multicolumn{3}{|c|}{\textit{Interior mixing parameters (S/R \texttt{Ri\_iwmix})} } \\
+      \hline
+        Riinfty & 0.7 & 
+           gradient Richardson number limit for shear instability \\
+        BVDQcon & -0.2E-4 1/s$^2$ & 
+           Brunt-V\"ai\"sal\"a squared \\
+        difm0 & 0.005 m$^2$/s & 
+           viscosity max. due to shear instability \\
+        difs0 & 0.005 m$^2$/s & 
+           tracer diffusivity max. due to shear instability \\
+        dift0 & 0.005 m$^2$/s & 
+           heat diffusivity max. due to shear instability \\
+        difmcon & 0.1 & 
+           viscosity due to convective instability \\
+        difscon & 0.1 & 
+           tracer diffusivity due to convective instability \\
+        diftcon & 0.1 & 
+           heat diffusivity due to convective instability \\
+      \hline
+      \multicolumn{3}{|c|}{\textit{Double-diffusive mixing parameters (S/R \texttt{ddmix})} } \\
+      \hline
+        Rrho0 & not used & 
+           limit for double diffusive density ratio \\
+        dsfmax & not used & 
+           maximum diffusivity in case of salt fingering \\
+         \hline
+      \hline
+    \end{tabular}
+  }
+  \caption{~}
+\end{table}
+
+
+
+%----------------------------------------------------------------------
+
+\subsubsection{Equations
+\label{sec:pkg:kpp:equations}}
+
+We restrict ourselves to writing out only the essential equations
+that relate to main processes and parameters mentioned above.
+We closely follow the notation of \cite{lar-eta:94}.
+
+\paragraph{Mixing in the boundary layer} ~ \\
+%
+~
+
+The vertical fluxes $\overline{wx}$
+of momentum and tracer properties $X$
+is composed of a gradient-flux term (proportional to
+the vertical property divergence $\partial_z X$), and
+a ``nonlocal'' term $\gamma_x$ that enhances the
+gradient-flux mixing coefficient $K_x$
+%
+\begin{equation}
+\overline{wx}(d) \, = \, -K_x \left(
+\frac{\partial X}{\partial z} \, - \, \gamma_x \right)
+\end{equation}
+
+\begin{itemize}
+%
+\item
+\textit{Boundary layer mixing profile} \\
+%
+It is expressed as the product of the boundary layer depth $h$,
+a depth-dependent turbulent velocity scale $w_x(\sigma)$ and a
+non-dimensional shape function $G(\sigma)$
+%
+\begin{equation}
+K_x(\sigma) \, = \, h \, w_x(\sigma) \, G(\sigma)
+\end{equation}
+%
+with dimensionless vertical coordinate $\sigma = d/h$.
+For details of $ w_x(\sigma)$ and $G(\sigma)$ we refer to
+\cite{lar-eta:94}.
+
+%
+\item
+\textit{Nonlocal mixing term} \\
+%
+The nonlocal transport term $\gamma$ is nonzero only for
+tracers in unstable (convective) forcing conditions.
+Thus, depending on the  stability parameter $\zeta = d/L$ 
+(with depth $d$, Monin-Obukhov length scale $L$)
+it has the following form:
+%
+\begin{eqnarray}
+\begin{array}{cl}
+\gamma_x \, = \, 0 & \zeta \, \ge \, 0 \\
+~ & ~ \\
+\left.
+\begin{array}{c}
+\gamma_m \, = \, 0 \\
+ ~ \\
+\gamma_s \, = \, C_s 
+\frac{\overline{w s_0}}{w_s(\sigma) h} \\
+ ~ \\
+\gamma_{\theta} \, = \, C_s
+\frac{\overline{w \theta_0}+\overline{w \theta_R}}{w_s(\sigma) h} \\
+\end{array}
+\right\} 
+&
+\zeta \, < \, 0 \\
+\end{array}
+\end{eqnarray}
+
+\end{itemize}
+
+
+\paragraph{Mixing in the interior} ~ \\
+%
+~
+
+\paragraph{Implicit time integration} ~ \\
+%
+~
+
+%----------------------------------------------------------------------
+
+\subsubsection{Key subroutines
+\label{sec:pkg:kpp:subroutines}}
+
+\paragraph{kpp\_calc:} Top-level routine. \\
+~
+
+\paragraph{kpp\_mix:} Intermediate-level routine \\
+~
+
+\paragraph{ri\_iwmix:} ~ \\
+%
+Compute interior viscosity and diffusivity coefficients due to
+%
+\begin{itemize}
+%
+\item
+shear instability (dependent on a local gradient Richardson number),
+%
+\item
+to background internal wave activity, and
+%
+\item
+to static instability (local Richardson number < 0).
+%
+\end{itemize}
+
+
+\paragraph{bldepth:} ~ \\
+%
+The oceanic planetary boundary layer depth, \texttt{hbl}, is determined as
+the shallowest depth where the bulk Richardson number is
+equal to the critical value, \texttt{Ricr}.
+
+Bulk Richardson numbers are evaluated by computing velocity and
+buoyancy differences between values at zgrid(kl) < 0 and surface
+reference values.
+In this configuration, the reference values are equal to the
+values in the surface layer.
+When using a very fine vertical grid, these values should be
+computed as the vertical average of velocity and buoyancy from
+the surface down to epsilon*zgrid(kl).
+
+When the bulk Richardson number at k exceeds Ricr, hbl is
+linearly interpolated between grid levels zgrid(k) and zgrid(k-1).
+
+The water column and the surface forcing are diagnosed for
+stable/ustable forcing conditions, and where hbl is relative
+to grid points (caseA), so that conditional branches can be
+avoided in later subroutines.
+
+\paragraph{blmix:} ~ \\
+%
+Compute boundary layer mixing coefficients.
+Mixing coefficients within boundary layer depend on surface
+forcing and the magnitude and gradient of interior mixing below
+the boundary layer ("matching").
+%
+\begin{enumerate}
+%
+\item
+compute velocity scales at hbl
+%
+\item
+find the interior viscosities and derivatives at hbl
+%
+\item
+compute turbulent velocity scales on the interfaces
+%
+\item
+compute the dimensionless shape functions at the interfaces
+%
+\item
+compute boundary layer diffusivities at the interfaces
+%
+\item
+compute nonlocal transport term
+%
+\item
+find diffusivities at kbl-1 grid level
+%
+\end{enumerate}
+
+\paragraph{kpp\_calc\_diff\_t/\_s, kpp\_calc\_visc:} ~  \\
+%
+Add contribution to net diffusivity/viscosity from 
+KPP diffusivity/viscosity.
+
+\paragraph{kpp\_transport\_t/\_s/\_ptr:} ~ \\
+%
+Add non local KPP transport term (ghat) to diffusive
+temperature/salinity/passive tracer flux.
+The nonlocal transport term is nonzero only for scalars
+in unstable (convective) forcing conditions. 
+
+\paragraph{Flow chart:} ~ \\
+%
+{\footnotesize
+\begin{verbatim}
+
+C     !CALLING SEQUENCE:
+c ...
+c  kpp_calc (TOP LEVEL ROUTINE)
+c  |
+c  |-- statekpp: o compute all EOS/density-related arrays
+c  |             o uses S/R FIND_ALPHA, FIND_BETA, FIND_RHO
+c  |
+c  |-- kppmix
+c  |   |--- ri_iwmix (compute interior mixing coefficients due to constant
+c  |   |              internal wave activity, static instability, 
+c  |   |              and local shear instability).
+c  |   |
+c  |   |--- bldepth (diagnose boundary layer depth)
+c  |   |
+c  |   |--- blmix (compute boundary layer diffusivities)
+c  |   |
+c  |   |--- enhance (enhance diffusivity at interface kbl - 1)
+c  |   o
+c  |
+c  |-- swfrac
+c  o
+
+\end{verbatim}
+}
+
+%----------------------------------------------------------------------
+
+\subsubsection{KPP diagnostics
+\label{sec:pkg:kpp:diagnostics}}
+
+Diagnostics output is available via the diagnostics package
+(see Section \ref{sec:pkg:diagnostics}).
+Available output fields are summarized in 
+Table \ref{tab:pkg:kpp:diagnostics}.
+
+\begin{table}[h!]
+\centering
+\label{tab:pkg:kpp:diagnostics}
+{\footnotesize
+\begin{verbatim}
+------------------------------------------------------
+ <-Name->|Levs|grid|<--  Units   -->|<- Tile (max=80c)
+------------------------------------------------------
+ KPPviscA| 23 |SM  |m^2/s           |KPP vertical eddy viscosity coefficient
+ KPPdiffS| 23 |SM  |m^2/s           |Vertical diffusion coefficient for salt & tracers
+ KPPdiffT| 23 |SM  |m^2/s           |Vertical diffusion coefficient for heat
+ KPPghat | 23 |SM  |s/m^2           |Nonlocal transport coefficient
+ KPPhbl  |  1 |SM  |m               |KPP boundary layer depth, bulk Ri criterion
+ KPPmld  |  1 |SM  |m               |Mixed layer depth, dT=.8degC density criterion
+ KPPfrac |  1 |SM  |                |Short-wave flux fraction penetrating mixing layer
+\end{verbatim}
+}
+\caption{~}
+\end{table}
+
+%----------------------------------------------------------------------
+
+\subsubsection{Reference experiments}
+
+lab\_sea:
+
+natl\_box:
+
+%----------------------------------------------------------------------
+
+\subsubsection{References}
+
