/[MITgcm]/manual/s_phys_pkgs/text/seaice.tex

Diff of /manual/s_phys_pkgs/text/seaice.tex

Parent Directory | Revision Log | View Revision Graph Revision Graph | View Patch Patch

-revision 1.18 by heimbach,
Sat Aug 13 15:50:05 2011 UTC
+revision 1.25 by mlosch,
Tue Mar 29 14:50:54 2016 UTC
 Line 61 
 no additional CPP options are required.
  (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
+ via CPP preprocessor flags. These options are set in
- \code{SEAICE\_OPTIONS.h} or in \code{ECCO\_CPPOPTIONS.h}.
+ \code{SEAICE\_OPTIONS.h}.
- Table \ref{tab:pkg:seaice:cpp} summarizes these options.
+ Table \ref{tab:pkg:seaice:cpp} summarizes the most important ones. For
+ more options see the default \code{pkg/seaice/SEAICE\_OPTIONS.h}.
  \begin{table}[!ht]
  \centering
-Line 80 
 Table \ref{tab:pkg:seaice:cpp} summarize
+Line 81 
 Table \ref{tab:pkg:seaice:cpp} summarize
          \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} &
+         \code{SEAICE\_ZETA\_SMOOTHREG} &
-           enable 8-category thermodynamics (by default undefined)\\
+           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} &
+         \code{SEAICE\_VARIABLE\_SALINITY} &
-           enable "salty" sea-ice (by default undefined) \\
+           enable sea-ice with variable salinity (by default undefined) \\
-         \code{SEAICE\_AGE} &
+         \code{SEAICE\_SITRACER} &
-           enable "age tracer" sea-ice (by default undefined) \\
+           enable sea-ice tracer package (by default undefined) \\
-         \code{SEAICE\_CAP\_HEFF} &
-           enable capping of sea-ice thickness to MAX\_HEFF \\ \hline
          \code{SEAICE\_BICE\_STRESS} &
            B-grid only for backward compatiblity: turn on ice-stress on
            ocean\\
-Line 105 
 Table \ref{tab:pkg:seaice:cpp} summarize
+Line 106 
 Table \ref{tab:pkg:seaice:cpp} summarize
        \hline
      \end{tabular}
    }
-   \caption{~}
+   \caption{Some of the most relevant CPP preprocessor flags in the
+     \code{seaice}-package.}
  \end{table}
  %----------------------------------------------------------------------
-Line 139 
 Table~\ref{tab:pkg:seaice:runtimeparms}
+Line 141 
 Table~\ref{tab:pkg:seaice:runtimeparms}
    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}
  %----------------------------------------------------------------------
 Line 182 
 viscous-plastic (VP) dynamic-thermodynam
  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
+ \item three different solution methods for solving the nonlinear
-   momentum equations have been adopted: LSOR \citep{zhang97}, and EVP
+   momentum equations have been adopted: LSOR \citep{zhang97}, EVP
-   \citep{hun97};
+   \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
 Line 296 
 air and ocean drag coefficients; $\rho_{
  densities; and $R_{air/ocean}$ are rotation matrices that act on the
  wind/current vectors.
- \paragraph{Viscous-Plastic (VP) Rheology and LSOR solver \label{sec:pkg:seaice:VPdynamics}}~\\
+ \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
 Line 317 
 The ice strain rate is given by
  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
 Line 349 
 recommended). For stress tensor computat
  = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state
  always lies on the elliptic yield curve by definition.
- In the current implementation, the VP-model is integrated with the
+ Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in
- semi-implicit line successive over relaxation (LSOR)-solver of
+ \code{SEAICE\_OPTIONS.h} before compiling replaces the method for
- \citet{zhang97}, which allows for long time steps that, in our case,
+ bounding $\zeta$ by a smooth (differentiable) expression:
- are limited by the explicit treatment of the Coriolis term. The
+ \begin{equation}
- explicit treatment of the Coriolis term does not represent a severe
+   \label{eq:zetaregsmooth}
- limitation because it restricts the time step to approximately the
+   \begin{split}
- same length as in the ocean model where the Coriolis term is also
+   \zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min})
- treated explicitly.
+       \,\zeta_{\max}}\right)\\
+   &= \frac{P}{2\Delta^*}
+   \tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right)
+   \end{split}
+ \end{equation}
+ 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.
+ \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 =
+.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.
+ Setting \code{SEAICEuseStrImpCpl = .TRUE.,} turns on ``strength
+ implicit coupling'' \citep{hutchings04} in the LSR-solver and in the
+ LSR-preconditioner for the JFNK-solver. In this mode, the different
+ contributions of the stress divergence terms are re-ordered in order
+ to increase the diagonal dominance of the system
+ matrix. Unfortunately, the convergence rate of the LSR solver is
+ increased only slightly, while the JFNK-convergence appears to be
+ unaffected.
  \paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\
  %
-Line 423 
 the damping time scale $T$ accordingly,
+Line 602 
 the damping time scale $T$ accordingly,
  $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly
  \code{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.
+ \paragraph{More stable variants of Elastic-Viscous-Plastic Dynamics:
+   EVP* , mEVP, and aEVP \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, we refer to these
+ variants by modified EVP (mEVP) and adaptive EVP (aEVP)
+ \citep{kimmritz16}. 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 extra ``intertial'' term
+ $m\,(u^{p+1}-u^{n})/\Delta{t})$ allows the definition of a residual
+ $|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to
+ $0$. In this way EVP can be re-interpreted as a pure iterative solver
+ where the sub-cycling has no association with time-relation (through
+ $\Delta{t}_{\mathrm{EVP}}$) \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}+\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} $\gg
+ (\alpha,\,\beta)$ \citep{kimmritz15}. Currently, there is no
+ termination criterion and the number of mEVP iterations is fixed to
+ \code{SEAICEnEVPstarSteps}.
+ In order to use mEVP in the MITgcm, set \code{SEAICEuseEVPstar =
+   .TRUE.,} in \code{data.seaice}. If \code{SEAICEuseEVPrev =.TRUE.,}
+ the actual form of equations (\ref{eq:evpstarsigma}) and
+ (\ref{eq:evpstarmom}) is used with fewer implicit terms and the factor
+ of $e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2})
+ and (\ref{eq:evpstresstensor12}). Although this modifies the original
+ EVP-equations, it turns out to improve convergence \citep{bouillon13}.
+ Another variant is the aEVP scheme \citep{kimmritz16}, where the value
+ of $\alpha$ is set dynamically based on the stability criterion
+ \begin{equation}
+   \label{eq:aevpalpha}
+   \alpha = \beta = \max\left( \tilde{c}\pi\sqrt{c \frac{\zeta}{A_{c}}
+     \frac{\Delta{t}}{\max(m,10^{-4}\text{\,kg})}},\alpha_{\min} \right)
+ \end{equation}
+ with the grid cell area $A_c$ and the ice and snow mass $m$.  This
+ choice sacrifices speed of convergence for stability with the result
+ that aEVP converges quickly to VP where $\alpha$ can be small and more
+ slowly in areas where the equations are stiff. In practice, aEVP leads
+ to an overall better convergence than mEVP \citep{kimmritz16}.
+ %
+ To use aEVP in the MITgcm set \code{SEAICEaEVPcoeff} $= \tilde{c}$;
+ this also sets the default values of \code{SEAICEaEVPcStar} ($c=4$)
+ and \code{SEAICEaEVPalphaMin} ($\alpha_{\min}=5$). Good convergence
+ has been obtained with setting these values \citep{kimmritz16}:
+ \code{SEAICEaEVPcoeff = 0.5, SEAICEnEVPstarSteps = 500,
+   SEAICEuseEVPstar = .TRUE., SEAICEuseEVPrev = .TRUE.}
+ Note, that probably because of the C-grid staggering of velocities and
+ stresses, mEVP may not converge as successfully as in
+ \citet{kimmritz15}, and that convergence at very high resolution
+ (order 5\,km) has not been studied yet.
  \paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\
  %
  In the so-called truncated ellipse method the shear viscosity $\eta$
-Line 431 
 is capped to suppress any tensile stress
+Line 688 
 is capped to suppress any tensile stress
    \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
+   {\sqrt{\max(\Delta_{\min}^{2},(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2
-       +4\dot{\epsilon}_{12}^2}}\right).
+       +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
-Line 532 
 relation $\sigma_{\alpha\beta} = 2\eta\d
+Line 789 
 relation $\sigma_{\alpha\beta} = 2\eta\d
  [(\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}
-Line 633 
 $\sigma_{21}^{Z}=0$, or equivalently $\e
+Line 891 
 $\sigma_{21}^{Z}=0$, or equivalently $\e
  \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\
  %
+ \noindent\textbf{NOTE: THIS SECTION IS TERRIBLY OUT OF DATE}\\
  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,
-Line 648 
 way to that of \citet{parkinson79} and \
+Line 907 
 way to that of \citet{parkinson79} and \
  The conductive heat flux depends strongly on the ice thickness $h$.
  However, the ice thickness in the model represents a mean over a
  potentially very heterogeneous thickness distribution.  In order to
- parameterize a sub-grid scale distribution for heat flux
+ parameterize a sub-grid scale distribution for heat flux computations,
- computations, the mean ice thickness $h$ is split into seven thickness
+ the mean ice thickness $h$ is split into $N$ thickness categories
- categories $H_{n}$ that are equally distributed between $2h$ and a
+ $H_{n}$ that are equally distributed between $2h$ and a minimum
- minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
+ imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$
- \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
+ for $n\in[1,N]$. The heat fluxes computed for each thickness category
- thickness category is area-averaged to give the total heat flux
+ is area-averaged to give the total heat flux \citep{hibler84}. To use
- \citep{hibler84}. To use this thickness category parameterization set
+ this thickness category parameterization set \code{SEAICE\_multDim} to
- \code{\#define SEAICE\_MULTICATEGORY}; note that this requires
+ the number of desired categories (7 is a good guess, for anything
- different restart files and switching this flag on in the middle of an
+ larger than 7 modify \code{SEAICE\_SIZE.h}) in
- integration is not possible.
+ \code{data.seaice}; note that this requires different restart files
+ and switching this flag on in the middle of an integration is not
+ advised. In order to include the same distribution for snow, set
+ \code{SEAICE\_useMultDimSnow = .TRUE.}; only then, the
+ parameterization of always having a fraction of thin ice is efficient
+ and generally thicker ice is produced \citep{castro-morales14}.
  The atmospheric heat flux is balanced by an oceanic heat flux from
  below.  The oceanic flux is proportional to

 Legend:



Removed from v.1.18
 


changed lines


 
Added in v.1.25
 Legend:



Removed from v.1.18
 


changed lines


 
Added in v.1.25
-Removed from v.1.18
+Added in v.1.25

	ViewVC Help
Powered by ViewVC 1.1.22