--- manual/s_phys_pkgs/text/seaice.tex	2011/02/28 16:05:07	1.14
+++ manual/s_phys_pkgs/text/seaice.tex	2015/01/21 18:18:22	1.22
@@ -1,4 +1,4 @@
-% $Header: /home/ubuntu/mnt/e9_copy/manual/s_phys_pkgs/text/seaice.tex,v 1.14 2011/02/28 16:05:07 mlosch Exp $
+% $Header: /home/ubuntu/mnt/e9_copy/manual/s_phys_pkgs/text/seaice.tex,v 1.22 2015/01/21 18:18:22 mlosch Exp $
 % $Name:  $
 
 %%EH3  Copied from "MITgcm/pkg/seaice/seaice_description.tex"
@@ -61,9 +61,9 @@
 (see Section \ref{sec:buildingCode}).
 
 Parts of the SEAICE code can be enabled or disabled at compile time
-via CPP preprocessor flags. These options are set in either
-\code{SEAICE\_OPTIONS.h} or in \code{ECCO\_CPPOPTIONS.h}.
-Table \ref{tab:pkg:seaice:cpp} summarizes these options.
+via CPP preprocessor flags. These options are set in 
+\code{SEAICE\_OPTIONS.h}.
+Table \ref{tab:pkg:seaice:cpp} summarizes the most important ones.
 
 \begin{table}[!ht]
 \centering
@@ -80,22 +80,22 @@
         \code{SEAICE\_CGRID} & 
           LSR solver on C-grid (rather than original B-grid) \\
         \code{SEAICE\_ALLOW\_EVP} & 
-          use EVP rather than LSR rheology solver \\
+          enable use of EVP rheology solver \\
+        \code{SEAICE\_ALLOW\_JFNK} & 
+          enable use of JFNK rheology solver \\
         \code{SEAICE\_EXTERNAL\_FLUXES} & 
           use EXF-computed fluxes as starting point \\
-        \code{SEAICE\_MULTICATEGORY} & 
-          enable 8-category thermodynamics (by default undefined)\\
+        \code{SEAICE\_ZETA\_SMOOTHREG} & 
+          use differentialable regularization for viscosities \\
         \code{SEAICE\_VARIABLE\_FREEZING\_POINT} & 
           enable linear dependence of the freezing point on salinity
           (by default undefined)\\
         \code{ALLOW\_SEAICE\_FLOODING} & 
           enable snow to ice conversion for submerged sea-ice \\
-        \code{SEAICE\_SALINITY} & 
-          enable "salty" sea-ice (by default undefined) \\
-        \code{SEAICE\_AGE} & 
-          enable "age tracer" sea-ice (by default undefined) \\
-        \code{SEAICE\_CAP\_HEFF} & 
-          enable capping of sea-ice thickness to MAX\_HEFF \\ \hline
+        \code{SEAICE\_VARIABLE\_SALINITY} & 
+          enable sea-ice with variable salinity (by default undefined) \\
+        \code{SEAICE\_SITRACER} & 
+          enable sea-ice tracer package (by default undefined) \\
         \code{SEAICE\_BICE\_STRESS} &
           B-grid only for backward compatiblity: turn on ice-stress on
           ocean\\
@@ -105,7 +105,8 @@
       \hline
     \end{tabular}
   }
-  \caption{~}
+  \caption{Some of the most relevant CPP preprocessor flags in the
+    \code{seaice}-package.} 
 \end{table}
 
 %----------------------------------------------------------------------
@@ -113,9 +114,9 @@
 \subsubsection{Run-time parameters
 \label{sec:pkg:seaice:runtime}}
 
-Run-time parameters are set in files 
-\code{data.pkg} (read in \code{packages\_readparms.F}),
-and \code{data.seaice} (read in \code{seaice\_readparms.F}).
+Run-time parameters (see Table~\ref{tab:pkg:seaice:runtimeparms}) are set
+in files \code{data.pkg} (read in \code{packages\_readparms.F}), and
+\code{data.seaice} (read in \code{seaice\_readparms.F}).
 
 \paragraph{Enabling the package}
 ~ \\
@@ -139,8 +140,6 @@
   over grid cell in meters; initializes variable \code{HSNOW};
 \item[\code{HsaltFile}:] Initial salinity of sea ice averaged over grid
   cell in g/m$^2$; initializes variable \code{HSALT};
-\item[\code{IceAgeFile}:] Initial ice age of sea ice averaged over grid
-  cell in seconds; initializes variable \code{ICEAGE};
 \end{description}
 
 %----------------------------------------------------------------------
@@ -182,14 +181,14 @@
 first introduced by \citet{hib79, hib80}. In order to adapt this model
 to the requirements of coupled ice-ocean state estimation, many
 important aspects of the original code have been modified and
-improved:
+improved \citep{losch10:_mitsim}:
 \begin{itemize}
 \item the code has been rewritten for an Arakawa C-grid, both B- and
   C-grid variants are available; the C-grid code allows for no-slip
   and free-slip lateral boundary conditions;
-\item two different solution methods for solving the nonlinear
-  momentum equations have been adopted: LSOR \citep{zhang97}, and EVP
-  \citep{hun97};
+\item three different solution methods for solving the nonlinear
+  momentum equations have been adopted: LSOR \citep{zhang97}, EVP
+  \citep{hun97}, JFNK \citep{lemieux10,losch14:_jfnk};
 \item ice-ocean stress can be formulated as in \citet{hibler87} or as in
   \citet{cam08}; 
 \item ice variables are advected by sophisticated, conservative
@@ -214,26 +213,39 @@
 to use the VP model as the default dynamic component of our ice
 model. To do this we extended the line successive over relaxation
 (LSOR) method of \citet{zhang97} for use in a parallel
-configuration.
+configuration. An EVP model and a free-drift implemtation can be
+selected with runtime flags.
 
-Note, that by default the seaice-package includes the orginial
+\paragraph{Compatibility with ice-thermodynamics package \code{thsice}\label{sec:pkg:seaice:thsice}}~\\
+%
+Note, that by default the \code{seaice}-package includes the orginial
 so-called zero-layer thermodynamics following \citet{hib80} with a
 snow cover as in \citet{zha98a}. The zero-layer thermodynamic model
 assumes that ice does not store heat and, therefore, tends to
 exaggerate the seasonal variability in ice thickness.  This
 exaggeration can be significantly reduced by using
-\citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic model
-that permits heat storage in ice.  Recently, the three-layer
-thermodynamic model has been reformulated by \citet{win00}.  The
-reformulation improves model physics by representing the brine content
-of the upper ice with a variable heat capacity.  It also improves
-model numerics and consumes less computer time and memory.  The Winton
-sea-ice thermodynamics have been ported to the MIT GCM; they currently
-reside under pkg/thsice. The package pkg/thsice is fully compatible
-with pkg/seaice and with pkg/exf. When turned on together with
-pkg/seaice, the zero-layer thermodynamics are replaced by the Winton
-thermodynamics.
+\citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic
+model that permits heat storage in ice. Recently, the three-layer thermodynamic model has been reformulated by
+\citet{win00}.  The reformulation improves model physics by
+representing the brine content of the upper ice with a variable heat
+capacity.  It also improves model numerics and consumes less computer
+time and memory. 
+
+The Winton sea-ice thermodynamics have been ported to the MIT GCM;
+they currently reside under \code{pkg/thsice}.  The package
+\code{thsice} is described in section~\ref{sec:pkg:thsice}; it is
+fully compatible with the packages \code{seaice} and \code{exf}. When
+turned on together with \code{seaice}, the zero-layer thermodynamics
+are replaced by the Winton thermodynamics. In order to use the
+\code{seaice}-package with the thermodynamics of \code{thsice},
+compile both packages and turn both package on in \code{data.pkg}; see
+an example in \code{global\_ocean.cs32x15/input.icedyn}. Note, that
+once \code{thsice} is turned on, the variables and diagnostics
+associated to the default thermodynamics are meaningless, and the
+diagnostics of \code{thsice} have to be used instead.
 
+\paragraph{Surface forcing\label{sec:pkg:seaice:surfaceforcing}}~\\
+%
 The sea ice model requires the following input fields: 10-m winds, 2-m
 air temperature and specific humidity, downward longwave and shortwave
 radiations, precipitation, evaporation, and river and glacier runoff.
@@ -244,8 +256,8 @@
 global: in ice-free regions bulk formulae are used to estimate oceanic
 forcing from the atmospheric fields.
 
-\paragraph{Dynamics\label{sec:pkg:seaice:dynamics}}
-
+\paragraph{Dynamics\label{sec:pkg:seaice:dynamics}}~\\
+%
 \newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}}
 \newcommand{\vtau}{\vek{\mathbf{\tau}}}
 The momentum equation of the sea-ice model is
@@ -283,6 +295,8 @@
 densities; and $R_{air/ocean}$ are rotation matrices that act on the
 wind/current vectors.
 
+\paragraph{Viscous-Plastic (VP) Rheology\label{sec:pkg:seaice:VPrheology}}~\\
+%
 For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
 be related to the ice strain rate and strength by a nonlinear
 viscous-plastic (VP) constitutive law \citep{hib79, zhang97}:
@@ -302,7 +316,7 @@
 The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on
 both thickness $h$ and compactness (concentration) $c$:
 \begin{equation}
-  P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},
+  P_{\max} = P^{*}c\,h\,\exp\{-C^{*}\cdot(1-c)\},
 \label{eq:icestrength}
 \end{equation}
 with the constants $P^{*}$ (run-time parameter \code{SEAICE\_strength}) and
@@ -334,28 +348,187 @@
 = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state
 always lies on the elliptic yield curve by definition.
 
-In the so-called truncated ellipse method the shear viscosity $\eta$
-is capped to suppress any tensile stress \citep{hibler97, geiger98}:
+Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in
+\code{SEAICE\_OPTIONS.h} before compiling replaces the method for
+bounding $\zeta$ by a smooth (differentiable) expression:
 \begin{equation}
-  \label{eq:etatem}
-  \eta = \min\left(\frac{\zeta}{e^2},
-  \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
-  {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
-      +4\dot{\epsilon}_{12}^2}}\right).
+  \label{eq:zetaregsmooth}
+  \begin{split}
+  \zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min})
+      \,\zeta_{\max}}\right)\\
+  &= \frac{P}{2\Delta^*}
+  \tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right) 
+  \end{split}
 \end{equation}
-To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in
-\code{SEAICE\_OPTIONS.h} and turn it on with
-\code{SEAICEuseTEM=.TRUE.} in \code{data.seaice}. 
+where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions
+by zero.
+
+\paragraph{LSR and  JFNK solver \label{sec:pkg:seaice:LSRJFNK}}~\\
+%
+% By default, the VP-model is integrated by a Pickwith the
+% semi-implicit line successive over relaxation (LSOR)-solver of
+% \citet{zhang97}, which allows for long time steps that, in our case,
+% are limited by the explicit treatment of the Coriolis term. The
+% explicit treatment of the Coriolis term does not represent a severe
+% limitation because it restricts the time step to approximately the
+% same length as in the ocean model where the Coriolis term is also
+% treated explicitly.
 
-In the current implementation, the VP-model is integrated with the
-semi-implicit line successive over relaxation (LSOR)-solver of
-\citet{zhang97}, which allows for long time steps that, in our case,
-are limited by the explicit treatment of the Coriolis term. The
-explicit treatment of the Coriolis term does not represent a severe
-limitation because it restricts the time step to approximately the
-same length as in the ocean model where the Coriolis term is also
-treated explicitly.
+\newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}} 
+%
+In the matrix notation, the discretized momentum equations can be
+written as
+\begin{equation}
+  \label{eq:matrixmom}
+  \mat{A}(\vek{x})\,\vek{x} = \vek{b}(\vek{x}).
+\end{equation}
+The solution vector $\vek{x}$ consists of the two velocity components
+$u$ and $v$ that contain the velocity variables at all grid points and
+at one time level. The standard (and default) method for solving
+Eq.\,(\ref{eq:matrixmom}) in the sea ice component of the
+\mbox{MITgcm}, as in many sea ice models, is an iterative Picard
+solver: in the $k$-th iteration a linearized form
+$\mat{A}(\vek{x}^{k-1})\,\vek{x}^{k} = \vek{b}(\vek{x}^{k-1})$ is
+solved (in the case of the MITgcm it is a Line Successive (over)
+Relaxation (LSR) algorithm \citep{zhang97}).  Picard solvers converge
+slowly, but generally the iteration is terminated after only a few
+non-linear steps \citep{zhang97, lemieux09} and the calculation
+continues with the next time level. This method is the default method
+in the MITgcm. The number of non-linear iteration steps or pseudo-time
+steps can be controlled by the runtime parameter
+\code{NPSEUDOTIMESTEPS} (default is 2).
+
+In order to overcome the poor convergence of the Picard-solver,
+\citet{lemieux10} introduced a Jacobian-free Newton-Krylov solver for
+the sea ice momentum equations. This solver is also implemented in the
+MITgcm \citep{losch14:_jfnk}. The Newton method transforms minimizing
+the residual $\vek{F}(\vek{x}) = \mat{A}(\vek{x})\,\vek{x} -
+\vek{b}(\vek{x})$ to finding the roots of a multivariate Taylor
+expansion of the residual \vek{F} around the previous ($k-1$) estimate
+$\vek{x}^{k-1}$:
+\begin{equation}
+  \label{eq:jfnktaylor}
+  \vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k}) = 
+  \vek{F}(\vek{x}^{k-1}) + \vek{F}'(\vek{x}^{k-1})\,\delta\vek{x}^{k}
+\end{equation}
+with the Jacobian $\mat{J}\equiv\vek{F}'$. The root
+$\vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k})=0$ is found by solving
+\begin{equation}
+  \label{eq:jfnklin}
+  \mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k} = -\vek{F}(\vek{x}^{k-1})
+\end{equation}
+for $\delta\vek{x}^{k}$. The next ($k$-th) estimate is given by
+$\vek{x}^{k}=\vek{x}^{k-1}+a\,\delta\vek{x}^{k}$. In order to avoid
+overshoots the factor $a$ is iteratively reduced in a line search
+($a=1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \ldots$) until
+$\|\vek{F}(\vek{x}^k)\| < \|\vek{F}(\vek{x}^{k-1})\|$, where
+$\|\cdot\|=\int\cdot\,dx^2$ is the $L_2$-norm. In practice, the line
+search is stopped at $a=\frac{1}{8}$. The line search starts after
+$\code{SEAICE\_JFNK\_lsIter}$ non-linear Newton iterations (off by
+default).
+
+
+Forming the Jacobian $\mat{J}$ explicitly is often avoided as ``too
+error prone and time consuming'' \citep{knoll04:_jfnk}. Instead,
+Krylov methods only require the action of \mat{J} on an arbitrary
+vector \vek{w} and hence allow a matrix free algorithm for solving
+Eq.\,(\ref{eq:jfnklin}) \citep{knoll04:_jfnk}. The action of \mat{J}
+can be approximated by a first-order Taylor series expansion:
+\begin{equation}
+  \label{eq:jfnkjacvecfd}
+  \mat{J}(\vek{x}^{k-1})\,\vek{w} \approx
+  \frac{\vek{F}(\vek{x}^{k-1}+\epsilon\vek{w}) - \vek{F}(\vek{x}^{k-1})}
+  {\epsilon} 
+\end{equation}
+or computed exactly with the help of automatic differentiation (AD)
+tools. \code{SEAICE\_JFNKepsilon} sets the step size
+$\epsilon$. 
+
+We use the Flexible Generalized Minimum RESidual method
+\citep[FGMRES,][]{saad93:_fgmres} with right-hand side preconditioning
+to solve Eq.\,(\ref{eq:jfnklin}) iteratively starting from a first
+guess of $\delta\vek{x}^{k}_{0} = 0$. For the preconditioning matrix
+\mat{P} we choose a simplified form of the system matrix
+$\mat{A}(\vek{x}^{k-1})$ \citep{lemieux10} where $\vek{x}^{k-1}$ is
+the estimate of the previous Newton step $k-1$. The transformed
+equation\,(\ref{eq:jfnklin}) becomes
+\begin{equation}
+  \label{eq:jfnklinpc}
+  \mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z} =
+  -\vek{F}(\vek{x}^{k-1}), 
+  \quad\text{with}\quad \delta\vek{z}=\mat{P}\delta\vek{x}^{k}.
+\end{equation}
+The Krylov method iteratively improves the approximate solution
+to~(\ref{eq:jfnklinpc}) in subspace ($\vek{r}_0$,
+$\mat{J}\mat{P}^{-1}\vek{r}_0$, $(\mat{J}\mat{P}^{-1})^2\vek{r}_0$,
+\ldots, $(\mat{J}\mat{P}^{-1})^m\vek{r}_0$) with increasing $m$;
+$\vek{r}_0 = -\vek{F}(\vek{x}^{k-1})
+-\mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k}_{0}$
+%-\vek{F}(\vek{x}^{k-1})
+%-\mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z}$ 
+is the initial residual of
+(\ref{eq:jfnklin}); $\vek{r}_0=-\vek{F}(\vek{x}^{k-1})$ with the first
+guess $\delta\vek{x}^{k}_{0}=0$. We allow a Krylov-subspace of
+dimension~$m=50$ and we do not use restarts. The preconditioning operation
+involves applying $\mat{P}^{-1}$ to the basis vectors $\vek{v}_0,
+\vek{v}_1, \vek{v}_2, \ldots, \vek{v}_m$ of the Krylov subspace. This
+operation is approximated by solving the linear system
+$\mat{P}\,\vek{w}=\vek{v}_i$. Because $\mat{P} \approx
+\mat{A}(\vek{x}^{k-1})$, we can use the LSR-algorithm \citep{zhang97}
+already implemented in the Picard solver. Each preconditioning
+operation uses a fixed number of 10~LSR-iterations avoiding any
+termination criterion. More details and results can be found in
+\citet{lemieux10, losch14:_jfnk}.
+
+To use the JFNK-solver set \code{SEAICEuseJFNK = .TRUE.} in the
+namelist file \code{data.seaice}; \code{SEAICE\_ALLOW\_JFNK} needs to
+be defined in \code{SEAICE\_OPTIONS.h} and we recommend using a smooth
+regularization of $\zeta$ by defining \code{SEAICE\_ZETA\_SMOOTHREG}
+(see above) for better convergence.  The non-linear Newton iteration
+is terminated when the $L_2$-norm of the residual is reduced by
+$\gamma_{\mathrm{nl}}$ (runtime parameter \code{JFNKgamma\_nonlin =
+  1.e-4} will already lead to expensive simulations) with respect to
+the initial norm: $\|\vek{F}(\vek{x}^k)\| <
+\gamma_{\mathrm{nl}}\|\vek{F}(\vek{x}^0)\|$.  Within a non-linear
+iteration, the linear FGMRES solver is terminated when the residual is
+smaller than $\gamma_k\|\vek{F}(\vek{x}^{k-1})\|$ where $\gamma_k$ is
+determined by
+\begin{equation}
+  \label{eq:jfnkgammalin}
+  \gamma_k = 
+  \begin{cases} 
+    \gamma_0 &\text{for $\|\vek{F}(\vek{x}^{k-1})\| \geq r$},  \\ 
+    \max\left(\gamma_{\min},
+    \frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)  
+%    \phi\left(\frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)^\alpha\right)  
+    &\text{for $\|\vek{F}(\vek{x}^{k-1})\| < r$,}
+  \end{cases}
+\end{equation}
+so that the linear tolerance parameter $\gamma_k$ decreases with the
+non-linear Newton step as the non-linear solution is approached. This
+inexact Newton method is generally more robust and computationally
+more efficient than exact methods \citep[e.g.,][]{knoll04:_jfnk}.
+% \footnote{The general idea behind
+%   inexact Newton methods is this: The Krylov solver is ``only''
+%   used to find an intermediate solution of the linear
+%   equation\,(\ref{eq:jfnklin}) that is used to improve the approximation of
+%   the actual equation\,(\ref{eq:matrixmom}). With the choice of a
+%   relatively weak lower limit for FGMRES convergence 
+%   $\gamma_{\min}$ we make sure that the time spent in the FGMRES
+%   solver is reduced at the cost of more Newton iterations. Newton
+%   iterations are cheaper than Krylov iterations so that this choice
+%   improves the overall efficiency.} 
+Typical parameter choices are
+$\gamma_0=\code{JFNKgamma\_lin\_max}=0.99$,
+$\gamma_{\min}=\code{JFNKgamma\_lin\_min}=0.1$, and $r = 
+\code{JFNKres\_tFac}\times\|\vek{F}(\vek{x}^{0})\|$ with
+$\code{JFNKres\_tFac} = \frac{1}{2}$. We recommend a maximum number of
+non-linear iterations $\code{SEAICEnewtonIterMax} = 100$ and a maximum
+number of Krylov iterations $\code{SEAICEkrylovIterMax} = 50$, because
+the Krylov subspace has a fixed dimension of 50.
 
+\paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\
+%
 \citet{hun97}'s introduced an elastic contribution to the strain
 rate in order to regularize Eq.~\ref{eq:vpequation} in such a way that
 the resulting elastic-viscous-plastic (EVP) and VP models are
@@ -373,11 +546,12 @@
 %used and compared the present sea-ice model study.
 The EVP-model uses an explicit time stepping scheme with a short
 timestep. According to the recommendation of \citet{hun97}, the
-EVP-model is stepped forward in time 120 times within the physical
-ocean model time step (although this parameter is under debate), to
-allow for elastic waves to disappear.  Because the scheme does not
-require a matrix inversion it is fast in spite of the small internal
-timestep and simple to implement on parallel computers
+EVP-model should be stepped forward in time 120 times
+($\code{SEAICE\_deltaTevp} = \code{SEAICIE\_deltaTdyn}/120$) within
+the physical ocean model time step (although this parameter is under
+debate), to allow for elastic waves to disappear.  Because the scheme
+does not require a matrix inversion it is fast in spite of the small
+internal timestep and simple to implement on parallel computers
 \citep{hun97}. For completeness, we repeat the equations for the
 components of the stress tensor $\sigma_{1} =
 \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
@@ -398,11 +572,12 @@
   \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
   &= \frac{P}{4T\Delta} D_S 
 \end{align}
-Here, the elastic parameter $E$ is redefined in terms of a damping timescale
-$T$ for elastic waves \[E=\frac{\zeta}{T}.\]
-$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and
-the external (long) timestep $\Delta{t}$. \citet{hun97} recommend
-$E_{0} = \frac{1}{3}$ (which is the default value in the code).
+Here, the elastic parameter $E$ is redefined in terms of a damping
+timescale $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
+$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and the external
+(long) timestep $\Delta{t}$.  $E_{0} = \frac{1}{3}$ is the default
+value in the code and close to what \citet{hun97} and
+\citet{hun01} recommend.
 
 To use the EVP solver, make sure that both \code{SEAICE\_CGRID} and
 \code{SEAICE\_ALLOW\_EVP} are defined in \code{SEAICE\_OPTIONS.h}
@@ -417,6 +592,77 @@
 $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly
 \code{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.
 
+\paragraph{More stable variant of Elastic-Viscous-Plastic Dynamics:  EVP*\label{sec:pkg:seaice:EVPstar}}~\\ 
+%
+The genuine EVP schemes appears to give noisy solutions \citep{hun01,
+  lemieux12, bouillon13}. This has lead to a modified EVP or EVP*
+\citep{lemieux12, bouillon13, kimmritz15}; here, refer to these
+variants by EVP*. The main idea is to modify the ``natural''
+time-discretization of the momentum equations:
+\begin{equation}
+  \label{eq:evpstar}
+  m\frac{D\vec{u}}{Dt} \approx m\frac{u^{p+1}-u^{n}}{\Delta{t}}
+  + \beta^{*}\frac{u^{p+1}-u^{p}}{\Delta{t}_{\mathrm{EVP}}}
+\end{equation}
+where $n$ is the previous time step index, and $p$ is the previous
+sub-cycling index. The term allows the definition of a residual
+$|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to
+$0$ and a re-interpretation of EVP as a pure iterative solver where
+the sub-cycling has lost all time-relation \citep{bouillon13,
+  kimmritz15}. Using the terminology of \citet{kimmritz15}, the
+evolution equations of stress $\sigma_{ij}$ and momentum $\vec{u}$ can
+be written as:
+\begin{align}
+  \label{eq:evpstarsigma}
+  \sigma_{ij}^{p+1}&=\sigma_{ij}^p+\frac{1}{\alpha}
+  \Big(\sigma_{ij}(\vec{u}^p)-\sigma_{ij}^p\Big),
+  \phantom{\int}\\
+  \label{eq:evpstarmom}
+  \vec{u}^{p+1}&=\vec{u}^p+\frac{1}{\beta}
+  \Big(\frac{\Delta t}{m}\nabla \cdot{\bf \sigma}^{p+1}+
+  \frac{\Delta t}{m}\vec{R}^{p+1/2}+\vec{u}_n-\vec{u}^p\Big).
+\end{align}
+$\vec{R}$ contains all terms in the momentum equations except for the
+rheology terms and the time derivative, $\alpha$ and $\beta$ are free
+parameters (\code{SEAICE\_evpAlpha}, \code{SEAICE\_evpBeta}) that
+replace the time stepping parameters \code{SEAICE\_deltaTevp}
+($\Delta{T}_{\mathrm{EVP}}$), \code{SEAICE\_elasticParm} ($E_{0}$), or
+\code{SEAICE\_evpTauRelax} ($T$). $\alpha$ and $\beta$ determine the
+speed of convergence and the stability. Usually, it makes sense to use
+$\alpha = \beta$, and \code{SEAICEnEVPstarSteps} $>> \alpha = \beta$
+\citep{kimmritz15}.
+
+In order to use EVP* in the MITgcm, set \code{SEAICEuseEVPstar =
+  .TRUE.,} in \code{data.seaice}. \code{SEAICEuseEVPrev =.TRUE.,} uses
+the actual form of equations (\ref{eq:evpstarsigma}) and
+(\ref{eq:evpstarmom}) with fewer implicit terms and the factor of
+$e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2})
+and (\ref{eq:evpstresstensor12}).  This turns out to improve
+convergence \citep{bouillon13}.
+
+Note, that for historical reasons, \code{SEAICE\_deltaTevp} needs to
+be set to some value in order to use also EVP*. Also note, that
+probably because of the C-grid staggering of velocities and stresses,
+EVP* does not converge as successfully as in \citet{kimmritz15}.
+
+
+\paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\
+%
+In the so-called truncated ellipse method the shear viscosity $\eta$
+is capped to suppress any tensile stress \citep{hibler97, geiger98}:
+\begin{equation}
+  \label{eq:etatem}
+  \eta = \min\left(\frac{\zeta}{e^2},
+  \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
+  {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
+      +4\dot{\epsilon}_{12}^2}}\right).
+\end{equation}
+To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in
+\code{SEAICE\_OPTIONS.h} and turn it on with
+\code{SEAICEuseTEM} in \code{data.seaice}. 
+
+\paragraph{Ice-Ocean stress \label{sec:pkg:seaice:iceoceanstress}}~\\
+%
 Moving sea ice exerts a stress on the ocean which is the opposite of
 the stress $\vtau_{ocean}$ in Eq.~\ref{eq:momseaice}. This stess is
 applied directly to the surface layer of the ocean model. An
@@ -446,7 +692,8 @@
 % $P$ at vorticity points.
 
 \paragraph{Finite-volume discretization of the stress tensor
-  divergence\label{sec:pkg:seaice:discretization}}
+  divergence\label{sec:pkg:seaice:discretization}}~\\
+%
 On an Arakawa C~grid, ice thickness and concentration and thus ice
 strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are
 naturally defined a C-points in the center of the grid
@@ -508,7 +755,8 @@
 [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2
 ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence
 $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is
-discretized in finite volumes. This conveniently avoids dealing with
+discretized in finite volumes \citep[see
+also][]{losch10:_mitsim}. This conveniently avoids dealing with
 further metric terms, as these are ``hidden'' in the differential cell
 widths. For the $u$-equation ($\alpha=1$) we have:
 \begin{align}
@@ -607,8 +855,8 @@
 analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set
 $\sigma_{21}^{Z}=0$, or equivalently $\eta_{i,j}^{Z}=0$, on boundaries.
 
-\paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}
-
+\paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\
+%
 In its original formulation the sea ice model \citep{menemenlis05}
 uses simple thermodynamics following the appendix of
 \citet{sem76}. This formulation does not allow storage of heat,
@@ -662,6 +910,8 @@
 \code{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter
 \code{SEAICEuseFlooding=.true.}.
 
+\paragraph{Advection of thermodynamic variables\label{sec:pkg:seaice:advection}}~\\
+%
 Effective ice thickness (ice volume per unit area,
 $c\cdot{h}$), concentration $c$ and effective snow thickness
 ($c\cdot{h}_{s}$) are advected by ice velocities:
@@ -681,34 +931,35 @@
 (negative thickness or concentration). These schemes conserve volume
 and horizontal area and are unconditionally stable, so that we can set
 $D_{X}=0$. Run-timeflags: \code{SEAICEadvScheme} (default=2, is the
-historic 2nd-order, centered difference scheme), \code{DIFF1}
+historic 2nd-order, centered difference scheme), \code{DIFF1} =
+$D_{X}/\Delta{x}$
 (default=0.004).
 
-There is considerable doubt about the reliability of a ``zero-layer''
-thermodynamic model --- \citet{semtner84} found significant errors in
-phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in
-such models --- so that today many sea ice models employ more complex
-thermodynamics. The MITgcm sea ice model provides the option to use
-the thermodynamics model of \citet{win00}, which in turn is based
-on the 3-layer model of \citet{sem76} and which treats brine
-content by means of enthalpy conservation. This scheme requires
-additional state variables, namely the enthalpy of the two ice layers
-(instead of effective ice salinity), to be advected by ice velocities.
+The MITgcm sea ice model provides the option to use
+the thermodynamics model of \citet{win00}, which in turn is based on
+the 3-layer model of \citet{sem76} and which treats brine content by
+means of enthalpy conservation; the corresponding package
+\code{thsice} is described in section~\ref{sec:pkg:thsice}. This
+scheme requires additional state variables, namely the enthalpy of the
+two ice layers (instead of effective ice salinity), to be advected by
+ice velocities.
 %
 The internal sea ice temperature is inferred from ice enthalpy.  To
 avoid unphysical (negative) values for ice thickness and
 concentration, a positive 2nd-order advection scheme with a SuperBee
-flux limiter \citep{roe:85} is used in this study to advect all
-sea-ice-related quantities of the \citet{win00} thermodynamic
-model.  Because of the non-linearity of the advection scheme, care
-must be taken in advecting these quantities: when simply using ice
-velocity to advect enthalpy, the total energy (i.e., the volume
-integral of enthalpy) is not conserved. Alternatively, one can advect
-the energy content (i.e., product of ice-volume and enthalpy) but then
-false enthalpy extrema can occur, which then leads to unrealistic ice
-temperature.  In the currently implemented solution, the sea-ice mass
-flux is used to advect the enthalpy in order to ensure conservation of
-enthalpy and to prevent false enthalpy extrema.
+flux limiter \citep{roe:85} should be used to advect all
+sea-ice-related quantities of the \citet{win00} thermodynamic model
+(runtime flag \code{thSIceAdvScheme=77} and
+\code{thSIce\_diffK}=$D_{X}$=0 in \code{data.ice}, defaults are 0).  Because of the
+non-linearity of the advection scheme, care must be taken in advecting
+these quantities: when simply using ice velocity to advect enthalpy,
+the total energy (i.e., the volume integral of enthalpy) is not
+conserved. Alternatively, one can advect the energy content (i.e.,
+product of ice-volume and enthalpy) but then false enthalpy extrema
+can occur, which then leads to unrealistic ice temperature.  In the
+currently implemented solution, the sea-ice mass flux is used to
+advect the enthalpy in order to ensure conservation of enthalpy and to
+prevent false enthalpy extrema. %
 
 %----------------------------------------------------------------------
 
@@ -774,58 +1025,7 @@
 Available output fields are summarized in 
 Table \ref{tab:pkg:seaice:diagnostics}.
 
-\begin{table}[!ht]
-\centering
-\label{tab:pkg:seaice:diagnostics}
-{\footnotesize
-\begin{verbatim}
----------+----+----+----------------+-----------------
- <-Name->|Levs|grid|<--  Units   -->|<- Tile (max=80c)
----------+----+----+----------------+-----------------
- SIarea  |  1 |SM  |m^2/m^2         |SEAICE fractional ice-covered area [0 to 1]
- SIheff  |  1 |SM  |m               |SEAICE effective ice thickness
- SIuice  |  1 |UU  |m/s             |SEAICE zonal ice velocity, >0 from West to East
- SIvice  |  1 |VV  |m/s             |SEAICE merid. ice velocity, >0 from South to North
- SIhsnow |  1 |SM  |m               |SEAICE snow thickness
- SIhsalt |  1 |SM  |g/m^2           |SEAICE effective salinity
- SIatmFW |  1 |SM  |kg/m^2/s        |Net freshwater flux from the atmosphere (+=down)
- SIuwind |  1 |SM  |m/s             |SEAICE zonal 10-m wind speed, >0 increases uVel
- SIvwind |  1 |SM  |m/s             |SEAICE meridional 10-m wind speed, >0 increases uVel
- SIfu    |  1 |UU  |N/m^2           |SEAICE zonal surface wind stress, >0 increases uVel
- SIfv    |  1 |VV  |N/m^2           |SEAICE merid. surface wind stress, >0 increases vVel
- SIempmr |  1 |SM  |kg/m^2/s        |SEAICE upward freshwater flux, > 0 increases salt
- SIqnet  |  1 |SM  |W/m^2           |SEAICE upward heatflux, turb+rad, >0 decreases theta
- SIqsw   |  1 |SM  |W/m^2           |SEAICE upward shortwave radiat., >0 decreases theta
- SIpress |  1 |SM  |m^2/s^2         |SEAICE strength (with upper and lower limit)
- SIzeta  |  1 |SM  |m^2/s           |SEAICE nonlinear bulk viscosity
- SIeta   |  1 |SM  |m^2/s           |SEAICE nonlinear shear viscosity
- SIsigI  |  1 |SM  |no units        |SEAICE normalized principle stress, component one
- SIsigII |  1 |SM  |no units        |SEAICE normalized principle stress, component two
- SIthdgrh|  1 |SM  |m/s             |SEAICE thermodynamic growth rate of effective ice thickness
- SIsnwice|  1 |SM  |m/s             |SEAICE ice formation rate due to flooding
- SIuheff |  1 |UU  |m^2/s           |Zonal Transport of effective ice thickness
- SIvheff |  1 |VV  |m^2/s           |Meridional Transport of effective ice thickness
- ADVxHEFF|  1 |UU  |m.m^2/s         |Zonal      Advective Flux of eff ice thickn
- ADVyHEFF|  1 |VV  |m.m^2/s         |Meridional Advective Flux of eff ice thickn
- DFxEHEFF|  1 |UU  |m.m^2/s         |Zonal      Diffusive Flux of eff ice thickn
- DFyEHEFF|  1 |VV  |m.m^2/s         |Meridional Diffusive Flux of eff ice thickn
- ADVxAREA|  1 |UU  |m^2/m^2.m^2/s   |Zonal      Advective Flux of fract area
- ADVyAREA|  1 |VV  |m^2/m^2.m^2/s   |Meridional Advective Flux of fract area
- DFxEAREA|  1 |UU  |m^2/m^2.m^2/s   |Zonal      Diffusive Flux of fract area
- DFyEAREA|  1 |VV  |m^2/m^2.m^2/s   |Meridional Diffusive Flux of fract area
- ADVxSNOW|  1 |UU  |m.m^2/s         |Zonal      Advective Flux of eff snow thickn
- ADVySNOW|  1 |VV  |m.m^2/s         |Meridional Advective Flux of eff snow thickn
- DFxESNOW|  1 |UU  |m.m^2/s         |Zonal      Diffusive Flux of eff snow thickn
- DFyESNOW|  1 |VV  |m.m^2/s         |Meridional Diffusive Flux of eff snow thickn
- ADVxSSLT|  1 |UU  |psu.m^2/s       |Zonal      Advective Flux of seaice salinity
- ADVySSLT|  1 |VV  |psu.m^2/s       |Meridional Advective Flux of seaice salinity
- DFxESSLT|  1 |UU  |psu.m^2/s       |Zonal      Diffusive Flux of seaice salinity
- DFyESSLT|  1 |VV  |psu.m^2/s       |Meridional Diffusive Flux of seaice salinity
-\end{verbatim}
-}
-\caption{Available diagnostics of the seaice-package}
-\end{table}
-
+\input{s_phys_pkgs/text/seaice_diags.tex}
 
 %\subsubsection{Package Reference}
 
@@ -833,7 +1033,13 @@
 \label{sec:pkg:seaice:experiments}
 
 \begin{itemize}
-\item{Labrador Sea experiment in lab\_sea verification directory. }
+\item{Labrador Sea experiment in \code{lab\_sea} verification directory. }
+\item \code{seaice\_obcs}, based on \code{lab\_sea}
+\item \code{offline\_exf\_seaice/input.seaicetd}, based on \code{lab\_sea}
+\item \code{global\_ocean.cs32x15/input.icedyn} and
+  \code{global\_ocean.cs32x15/input.seaice}, global
+  cubed-sphere-experiment with combinations of \code{seaice} and
+  \code{thsice}
 \end{itemize}
 
 
