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

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

Parent Directory Parent Directory | Revision Log 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. Line 61  no additional CPP options are required.
61  (see Section \ref{sec:buildingCode}).  (see Section \ref{sec:buildingCode}).
62    
63  Parts of the SEAICE code can be enabled or disabled at compile time  Parts of the SEAICE code can be enabled or disabled at compile time
64  via CPP preprocessor flags. These options are set in either  via CPP preprocessor flags. These options are set in
65  \code{SEAICE\_OPTIONS.h} or in \code{ECCO\_CPPOPTIONS.h}.  \code{SEAICE\_OPTIONS.h}.
66  Table \ref{tab:pkg:seaice:cpp} summarizes these options.  Table \ref{tab:pkg:seaice:cpp} summarizes the most important ones. For
67    more options see the default \code{pkg/seaice/SEAICE\_OPTIONS.h}.
68    
69  \begin{table}[!ht]  \begin{table}[!ht]
70  \centering  \centering
# Line 80  Table \ref{tab:pkg:seaice:cpp} summarize Line 81  Table \ref{tab:pkg:seaice:cpp} summarize
81          \code{SEAICE\_CGRID} &          \code{SEAICE\_CGRID} &
82            LSR solver on C-grid (rather than original B-grid) \\            LSR solver on C-grid (rather than original B-grid) \\
83          \code{SEAICE\_ALLOW\_EVP} &          \code{SEAICE\_ALLOW\_EVP} &
84            use EVP rather than LSR rheology solver \\            enable use of EVP rheology solver \\
85            \code{SEAICE\_ALLOW\_JFNK} &
86              enable use of JFNK rheology solver \\
87          \code{SEAICE\_EXTERNAL\_FLUXES} &          \code{SEAICE\_EXTERNAL\_FLUXES} &
88            use EXF-computed fluxes as starting point \\            use EXF-computed fluxes as starting point \\
89          \code{SEAICE\_MULTICATEGORY} &          \code{SEAICE\_ZETA\_SMOOTHREG} &
90            enable 8-category thermodynamics (by default undefined)\\            use differentialable regularization for viscosities \\
91          \code{SEAICE\_VARIABLE\_FREEZING\_POINT} &          \code{SEAICE\_VARIABLE\_FREEZING\_POINT} &
92            enable linear dependence of the freezing point on salinity            enable linear dependence of the freezing point on salinity
93            (by default undefined)\\            (by default undefined)\\
94          \code{ALLOW\_SEAICE\_FLOODING} &          \code{ALLOW\_SEAICE\_FLOODING} &
95            enable snow to ice conversion for submerged sea-ice \\            enable snow to ice conversion for submerged sea-ice \\
96          \code{SEAICE\_SALINITY} &          \code{SEAICE\_VARIABLE\_SALINITY} &
97            enable "salty" sea-ice (by default undefined) \\            enable sea-ice with variable salinity (by default undefined) \\
98          \code{SEAICE\_AGE} &          \code{SEAICE\_SITRACER} &
99            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  
100          \code{SEAICE\_BICE\_STRESS} &          \code{SEAICE\_BICE\_STRESS} &
101            B-grid only for backward compatiblity: turn on ice-stress on            B-grid only for backward compatiblity: turn on ice-stress on
102            ocean\\            ocean\\
# Line 105  Table \ref{tab:pkg:seaice:cpp} summarize Line 106  Table \ref{tab:pkg:seaice:cpp} summarize
106        \hline        \hline
107      \end{tabular}      \end{tabular}
108    }    }
109    \caption{~}    \caption{Some of the most relevant CPP preprocessor flags in the
110        \code{seaice}-package.}
111  \end{table}  \end{table}
112    
113  %----------------------------------------------------------------------  %----------------------------------------------------------------------
# Line 139  Table~\ref{tab:pkg:seaice:runtimeparms} Line 141  Table~\ref{tab:pkg:seaice:runtimeparms}
141    over grid cell in meters; initializes variable \code{HSNOW};    over grid cell in meters; initializes variable \code{HSNOW};
142  \item[\code{HsaltFile}:] Initial salinity of sea ice averaged over grid  \item[\code{HsaltFile}:] Initial salinity of sea ice averaged over grid
143    cell in g/m$^2$; initializes variable \code{HSALT};    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};  
144  \end{description}  \end{description}
145    
146  %----------------------------------------------------------------------  %----------------------------------------------------------------------
# Line 182  viscous-plastic (VP) dynamic-thermodynam Line 182  viscous-plastic (VP) dynamic-thermodynam
182  first introduced by \citet{hib79, hib80}. In order to adapt this model  first introduced by \citet{hib79, hib80}. In order to adapt this model
183  to the requirements of coupled ice-ocean state estimation, many  to the requirements of coupled ice-ocean state estimation, many
184  important aspects of the original code have been modified and  important aspects of the original code have been modified and
185  improved:  improved \citep{losch10:_mitsim}:
186  \begin{itemize}  \begin{itemize}
187  \item the code has been rewritten for an Arakawa C-grid, both B- and  \item the code has been rewritten for an Arakawa C-grid, both B- and
188    C-grid variants are available; the C-grid code allows for no-slip    C-grid variants are available; the C-grid code allows for no-slip
189    and free-slip lateral boundary conditions;    and free-slip lateral boundary conditions;
190  \item two different solution methods for solving the nonlinear  \item three different solution methods for solving the nonlinear
191    momentum equations have been adopted: LSOR \citep{zhang97}, and EVP    momentum equations have been adopted: LSOR \citep{zhang97}, EVP
192    \citep{hun97};    \citep{hun97}, JFNK \citep{lemieux10,losch14:_jfnk};
193  \item ice-ocean stress can be formulated as in \citet{hibler87} or as in  \item ice-ocean stress can be formulated as in \citet{hibler87} or as in
194    \citet{cam08};    \citet{cam08};
195  \item ice variables are advected by sophisticated, conservative  \item ice variables are advected by sophisticated, conservative
# Line 296  air and ocean drag coefficients; $\rho_{ Line 296  air and ocean drag coefficients; $\rho_{
296  densities; and $R_{air/ocean}$ are rotation matrices that act on the  densities; and $R_{air/ocean}$ are rotation matrices that act on the
297  wind/current vectors.  wind/current vectors.
298    
299  \paragraph{Viscous-Plastic (VP) Rheology and LSOR solver \label{sec:pkg:seaice:VPdynamics}}~\\  \paragraph{Viscous-Plastic (VP) Rheology\label{sec:pkg:seaice:VPrheology}}~\\
300  %  %
301  For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can  For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
302  be related to the ice strain rate and strength by a nonlinear  be related to the ice strain rate and strength by a nonlinear
# Line 317  The ice strain rate is given by Line 317  The ice strain rate is given by
317  The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on  The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on
318  both thickness $h$ and compactness (concentration) $c$:  both thickness $h$ and compactness (concentration) $c$:
319  \begin{equation}  \begin{equation}
320    P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},    P_{\max} = P^{*}c\,h\,\exp\{-C^{*}\cdot(1-c)\},
321  \label{eq:icestrength}  \label{eq:icestrength}
322  \end{equation}  \end{equation}
323  with the constants $P^{*}$ (run-time parameter \code{SEAICE\_strength}) and  with the constants $P^{*}$ (run-time parameter \code{SEAICE\_strength}) and
# Line 349  recommended). For stress tensor computat Line 349  recommended). For stress tensor computat
349  = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state  = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state
350  always lies on the elliptic yield curve by definition.  always lies on the elliptic yield curve by definition.
351    
352  In the current implementation, the VP-model is integrated with the  Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in
353  semi-implicit line successive over relaxation (LSOR)-solver of  \code{SEAICE\_OPTIONS.h} before compiling replaces the method for
354  \citet{zhang97}, which allows for long time steps that, in our case,  bounding $\zeta$ by a smooth (differentiable) expression:
355  are limited by the explicit treatment of the Coriolis term. The  \begin{equation}
356  explicit treatment of the Coriolis term does not represent a severe    \label{eq:zetaregsmooth}
357  limitation because it restricts the time step to approximately the    \begin{split}
358  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})
359  treated explicitly.        \,\zeta_{\max}}\right)\\
360      &= \frac{P}{2\Delta^*}
361      \tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right)
362      \end{split}
363    \end{equation}
364    where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions
365    by zero.
366    
367    \paragraph{LSR and  JFNK solver \label{sec:pkg:seaice:LSRJFNK}}~\\
368    %
369    % By default, the VP-model is integrated by a Pickwith the
370    % semi-implicit line successive over relaxation (LSOR)-solver of
371    % \citet{zhang97}, which allows for long time steps that, in our case,
372    % are limited by the explicit treatment of the Coriolis term. The
373    % explicit treatment of the Coriolis term does not represent a severe
374    % limitation because it restricts the time step to approximately the
375    % same length as in the ocean model where the Coriolis term is also
376    % treated explicitly.
377    
378    \newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}}
379    %
380    In the matrix notation, the discretized momentum equations can be
381    written as
382    \begin{equation}
383      \label{eq:matrixmom}
384      \mat{A}(\vek{x})\,\vek{x} = \vek{b}(\vek{x}).
385    \end{equation}
386    The solution vector $\vek{x}$ consists of the two velocity components
387    $u$ and $v$ that contain the velocity variables at all grid points and
388    at one time level. The standard (and default) method for solving
389    Eq.\,(\ref{eq:matrixmom}) in the sea ice component of the
390    \mbox{MITgcm}, as in many sea ice models, is an iterative Picard
391    solver: in the $k$-th iteration a linearized form
392    $\mat{A}(\vek{x}^{k-1})\,\vek{x}^{k} = \vek{b}(\vek{x}^{k-1})$ is
393    solved (in the case of the MITgcm it is a Line Successive (over)
394    Relaxation (LSR) algorithm \citep{zhang97}).  Picard solvers converge
395    slowly, but generally the iteration is terminated after only a few
396    non-linear steps \citep{zhang97, lemieux09} and the calculation
397    continues with the next time level. This method is the default method
398    in the MITgcm. The number of non-linear iteration steps or pseudo-time
399    steps can be controlled by the runtime parameter
400    \code{NPSEUDOTIMESTEPS} (default is 2).
401    
402    In order to overcome the poor convergence of the Picard-solver,
403    \citet{lemieux10} introduced a Jacobian-free Newton-Krylov solver for
404    the sea ice momentum equations. This solver is also implemented in the
405    MITgcm \citep{losch14:_jfnk}. The Newton method transforms minimizing
406    the residual $\vek{F}(\vek{x}) = \mat{A}(\vek{x})\,\vek{x} -
407    \vek{b}(\vek{x})$ to finding the roots of a multivariate Taylor
408    expansion of the residual \vek{F} around the previous ($k-1$) estimate
409    $\vek{x}^{k-1}$:
410    \begin{equation}
411      \label{eq:jfnktaylor}
412      \vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k}) =
413      \vek{F}(\vek{x}^{k-1}) + \vek{F}'(\vek{x}^{k-1})\,\delta\vek{x}^{k}
414    \end{equation}
415    with the Jacobian $\mat{J}\equiv\vek{F}'$. The root
416    $\vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k})=0$ is found by solving
417    \begin{equation}
418      \label{eq:jfnklin}
419      \mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k} = -\vek{F}(\vek{x}^{k-1})
420    \end{equation}
421    for $\delta\vek{x}^{k}$. The next ($k$-th) estimate is given by
422    $\vek{x}^{k}=\vek{x}^{k-1}+a\,\delta\vek{x}^{k}$. In order to avoid
423    overshoots the factor $a$ is iteratively reduced in a line search
424    ($a=1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \ldots$) until
425    $\|\vek{F}(\vek{x}^k)\| < \|\vek{F}(\vek{x}^{k-1})\|$, where
426    $\|\cdot\|=\int\cdot\,dx^2$ is the $L_2$-norm. In practice, the line
427    search is stopped at $a=\frac{1}{8}$. The line search starts after
428    $\code{SEAICE\_JFNK\_lsIter}$ non-linear Newton iterations (off by
429    default).
430    
431    
432    Forming the Jacobian $\mat{J}$ explicitly is often avoided as ``too
433    error prone and time consuming'' \citep{knoll04:_jfnk}. Instead,
434    Krylov methods only require the action of \mat{J} on an arbitrary
435    vector \vek{w} and hence allow a matrix free algorithm for solving
436    Eq.\,(\ref{eq:jfnklin}) \citep{knoll04:_jfnk}. The action of \mat{J}
437    can be approximated by a first-order Taylor series expansion:
438    \begin{equation}
439      \label{eq:jfnkjacvecfd}
440      \mat{J}(\vek{x}^{k-1})\,\vek{w} \approx
441      \frac{\vek{F}(\vek{x}^{k-1}+\epsilon\vek{w}) - \vek{F}(\vek{x}^{k-1})}
442      {\epsilon}
443    \end{equation}
444    or computed exactly with the help of automatic differentiation (AD)
445    tools. \code{SEAICE\_JFNKepsilon} sets the step size
446    $\epsilon$.
447    
448    We use the Flexible Generalized Minimum RESidual method
449    \citep[FGMRES,][]{saad93:_fgmres} with right-hand side preconditioning
450    to solve Eq.\,(\ref{eq:jfnklin}) iteratively starting from a first
451    guess of $\delta\vek{x}^{k}_{0} = 0$. For the preconditioning matrix
452    \mat{P} we choose a simplified form of the system matrix
453    $\mat{A}(\vek{x}^{k-1})$ \citep{lemieux10} where $\vek{x}^{k-1}$ is
454    the estimate of the previous Newton step $k-1$. The transformed
455    equation\,(\ref{eq:jfnklin}) becomes
456    \begin{equation}
457      \label{eq:jfnklinpc}
458      \mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z} =
459      -\vek{F}(\vek{x}^{k-1}),
460      \quad\text{with}\quad \delta\vek{z}=\mat{P}\delta\vek{x}^{k}.
461    \end{equation}
462    The Krylov method iteratively improves the approximate solution
463    to~(\ref{eq:jfnklinpc}) in subspace ($\vek{r}_0$,
464    $\mat{J}\mat{P}^{-1}\vek{r}_0$, $(\mat{J}\mat{P}^{-1})^2\vek{r}_0$,
465    \ldots, $(\mat{J}\mat{P}^{-1})^m\vek{r}_0$) with increasing $m$;
466    $\vek{r}_0 = -\vek{F}(\vek{x}^{k-1})
467    -\mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k}_{0}$
468    %-\vek{F}(\vek{x}^{k-1})
469    %-\mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z}$
470    is the initial residual of
471    (\ref{eq:jfnklin}); $\vek{r}_0=-\vek{F}(\vek{x}^{k-1})$ with the first
472    guess $\delta\vek{x}^{k}_{0}=0$. We allow a Krylov-subspace of
473    dimension~$m=50$ and we do not use restarts. The preconditioning operation
474    involves applying $\mat{P}^{-1}$ to the basis vectors $\vek{v}_0,
475    \vek{v}_1, \vek{v}_2, \ldots, \vek{v}_m$ of the Krylov subspace. This
476    operation is approximated by solving the linear system
477    $\mat{P}\,\vek{w}=\vek{v}_i$. Because $\mat{P} \approx
478    \mat{A}(\vek{x}^{k-1})$, we can use the LSR-algorithm \citep{zhang97}
479    already implemented in the Picard solver. Each preconditioning
480    operation uses a fixed number of 10~LSR-iterations avoiding any
481    termination criterion. More details and results can be found in
482    \citet{lemieux10, losch14:_jfnk}.
483    
484    To use the JFNK-solver set \code{SEAICEuseJFNK = .TRUE.} in the
485    namelist file \code{data.seaice}; \code{SEAICE\_ALLOW\_JFNK} needs to
486    be defined in \code{SEAICE\_OPTIONS.h} and we recommend using a smooth
487    regularization of $\zeta$ by defining \code{SEAICE\_ZETA\_SMOOTHREG}
488    (see above) for better convergence.  The non-linear Newton iteration
489    is terminated when the $L_2$-norm of the residual is reduced by
490    $\gamma_{\mathrm{nl}}$ (runtime parameter \code{JFNKgamma\_nonlin =
491      1.e-4} will already lead to expensive simulations) with respect to
492    the initial norm: $\|\vek{F}(\vek{x}^k)\| <
493    \gamma_{\mathrm{nl}}\|\vek{F}(\vek{x}^0)\|$.  Within a non-linear
494    iteration, the linear FGMRES solver is terminated when the residual is
495    smaller than $\gamma_k\|\vek{F}(\vek{x}^{k-1})\|$ where $\gamma_k$ is
496    determined by
497    \begin{equation}
498      \label{eq:jfnkgammalin}
499      \gamma_k =
500      \begin{cases}
501        \gamma_0 &\text{for $\|\vek{F}(\vek{x}^{k-1})\| \geq r$},  \\
502        \max\left(\gamma_{\min},
503        \frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)  
504    %    \phi\left(\frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)^\alpha\right)  
505        &\text{for $\|\vek{F}(\vek{x}^{k-1})\| < r$,}
506      \end{cases}
507    \end{equation}
508    so that the linear tolerance parameter $\gamma_k$ decreases with the
509    non-linear Newton step as the non-linear solution is approached. This
510    inexact Newton method is generally more robust and computationally
511    more efficient than exact methods \citep[e.g.,][]{knoll04:_jfnk}.
512    % \footnote{The general idea behind
513    %   inexact Newton methods is this: The Krylov solver is ``only''
514    %   used to find an intermediate solution of the linear
515    %   equation\,(\ref{eq:jfnklin}) that is used to improve the approximation of
516    %   the actual equation\,(\ref{eq:matrixmom}). With the choice of a
517    %   relatively weak lower limit for FGMRES convergence
518    %   $\gamma_{\min}$ we make sure that the time spent in the FGMRES
519    %   solver is reduced at the cost of more Newton iterations. Newton
520    %   iterations are cheaper than Krylov iterations so that this choice
521    %   improves the overall efficiency.}
522    Typical parameter choices are
523    $\gamma_0=\code{JFNKgamma\_lin\_max}=0.99$,
524    $\gamma_{\min}=\code{JFNKgamma\_lin\_min}=0.1$, and $r =
525    \code{JFNKres\_tFac}\times\|\vek{F}(\vek{x}^{0})\|$ with
526    $\code{JFNKres\_tFac} = \frac{1}{2}$. We recommend a maximum number of
527    non-linear iterations $\code{SEAICEnewtonIterMax} = 100$ and a maximum
528    number of Krylov iterations $\code{SEAICEkrylovIterMax} = 50$, because
529    the Krylov subspace has a fixed dimension of 50.
530    
531    Setting \code{SEAICEuseStrImpCpl = .TRUE.,} turns on ``strength
532    implicit coupling'' \citep{hutchings04} in the LSR-solver and in the
533    LSR-preconditioner for the JFNK-solver. In this mode, the different
534    contributions of the stress divergence terms are re-ordered in order
535    to increase the diagonal dominance of the system
536    matrix. Unfortunately, the convergence rate of the LSR solver is
537    increased only slightly, while the JFNK-convergence appears to be
538    unaffected.
539    
540  \paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\  \paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\
541  %  %
# Line 423  the damping time scale $T$ accordingly, Line 602  the damping time scale $T$ accordingly,
602  $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly  $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly
603  \code{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.  \code{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.
604    
605    \paragraph{More stable variants of Elastic-Viscous-Plastic Dynamics:
606      EVP* , mEVP, and aEVP \label{sec:pkg:seaice:EVPstar}}~\\
607    %
608    The genuine EVP schemes appears to give noisy solutions \citep{hun01,
609      lemieux12, bouillon13}. This has lead to a modified EVP or EVP*
610    \citep{lemieux12, bouillon13, kimmritz15}; here, we refer to these
611    variants by modified EVP (mEVP) and adaptive EVP (aEVP)
612    \citep{kimmritz16}. The main idea is to modify the ``natural''
613    time-discretization of the momentum equations:
614    \begin{equation}
615      \label{eq:evpstar}
616      m\frac{D\vec{u}}{Dt} \approx m\frac{u^{p+1}-u^{n}}{\Delta{t}}
617      + \beta^{*}\frac{u^{p+1}-u^{p}}{\Delta{t}_{\mathrm{EVP}}}
618    \end{equation}
619    where $n$ is the previous time step index, and $p$ is the previous
620    sub-cycling index. The extra ``intertial'' term
621    $m\,(u^{p+1}-u^{n})/\Delta{t})$ allows the definition of a residual
622    $|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to
623    $0$. In this way EVP can be re-interpreted as a pure iterative solver
624    where the sub-cycling has no association with time-relation (through
625    $\Delta{t}_{\mathrm{EVP}}$) \citep{bouillon13, kimmritz15}. Using the
626    terminology of \citet{kimmritz15}, the evolution equations of stress
627    $\sigma_{ij}$ and momentum $\vec{u}$ can be written as:
628    \begin{align}
629      \label{eq:evpstarsigma}
630      \sigma_{ij}^{p+1}&=\sigma_{ij}^p+\frac{1}{\alpha}
631      \Big(\sigma_{ij}(\vec{u}^p)-\sigma_{ij}^p\Big),
632      \phantom{\int}\\
633      \label{eq:evpstarmom}
634      \vec{u}^{p+1}&=\vec{u}^p+\frac{1}{\beta}
635      \Big(\frac{\Delta t}{m}\nabla \cdot{\bf \sigma}^{p+1}+
636      \frac{\Delta t}{m}\vec{R}^{p}+\vec{u}_n-\vec{u}^p\Big).
637    \end{align}
638    $\vec{R}$ contains all terms in the momentum equations except for the
639    rheology terms and the time derivative; $\alpha$ and $\beta$ are free
640    parameters (\code{SEAICE\_evpAlpha}, \code{SEAICE\_evpBeta}) that
641    replace the time stepping parameters \code{SEAICE\_deltaTevp}
642    ($\Delta{T}_{\mathrm{EVP}}$), \code{SEAICE\_elasticParm} ($E_{0}$), or
643    \code{SEAICE\_evpTauRelax} ($T$). $\alpha$ and $\beta$ determine the
644    speed of convergence and the stability. Usually, it makes sense to use
645    $\alpha = \beta$, and \code{SEAICEnEVPstarSteps} $\gg
646    (\alpha,\,\beta)$ \citep{kimmritz15}. Currently, there is no
647    termination criterion and the number of mEVP iterations is fixed to
648    \code{SEAICEnEVPstarSteps}.
649    
650    In order to use mEVP in the MITgcm, set \code{SEAICEuseEVPstar =
651      .TRUE.,} in \code{data.seaice}. If \code{SEAICEuseEVPrev =.TRUE.,}
652    the actual form of equations (\ref{eq:evpstarsigma}) and
653    (\ref{eq:evpstarmom}) is used with fewer implicit terms and the factor
654    of $e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2})
655    and (\ref{eq:evpstresstensor12}). Although this modifies the original
656    EVP-equations, it turns out to improve convergence \citep{bouillon13}.
657    
658    Another variant is the aEVP scheme \citep{kimmritz16}, where the value
659    of $\alpha$ is set dynamically based on the stability criterion
660    \begin{equation}
661      \label{eq:aevpalpha}
662      \alpha = \beta = \max\left( \tilde{c}\pi\sqrt{c \frac{\zeta}{A_{c}}
663        \frac{\Delta{t}}{\max(m,10^{-4}\text{\,kg})}},\alpha_{\min} \right)
664    \end{equation}
665    with the grid cell area $A_c$ and the ice and snow mass $m$.  This
666    choice sacrifices speed of convergence for stability with the result
667    that aEVP converges quickly to VP where $\alpha$ can be small and more
668    slowly in areas where the equations are stiff. In practice, aEVP leads
669    to an overall better convergence than mEVP \citep{kimmritz16}.
670    %
671    To use aEVP in the MITgcm set \code{SEAICEaEVPcoeff} $= \tilde{c}$;
672    this also sets the default values of \code{SEAICEaEVPcStar} ($c=4$)
673    and \code{SEAICEaEVPalphaMin} ($\alpha_{\min}=5$). Good convergence
674    has been obtained with setting these values \citep{kimmritz16}:
675    \code{SEAICEaEVPcoeff = 0.5, SEAICEnEVPstarSteps = 500,
676      SEAICEuseEVPstar = .TRUE., SEAICEuseEVPrev = .TRUE.}
677    
678    Note, that probably because of the C-grid staggering of velocities and
679    stresses, mEVP may not converge as successfully as in
680    \citet{kimmritz15}, and that convergence at very high resolution
681    (order 5\,km) has not been studied yet.
682    
683  \paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\  \paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\
684  %  %
685  In the so-called truncated ellipse method the shear viscosity $\eta$  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
688    \label{eq:etatem}    \label{eq:etatem}
689    \eta = \min\left(\frac{\zeta}{e^2},    \eta = \min\left(\frac{\zeta}{e^2},
690    \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}    \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
691    {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2    {\sqrt{\max(\Delta_{\min}^{2},(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2
692        +4\dot{\epsilon}_{12}^2}}\right).        +4\dot{\epsilon}_{12}^2})}\right).
693  \end{equation}  \end{equation}
694  To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in  To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in
695  \code{SEAICE\_OPTIONS.h} and turn it on with  \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
789  [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2  [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2
790  ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence  ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence
791  $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is  $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is
792  discretized in finite volumes. This conveniently avoids dealing with  discretized in finite volumes \citep[see
793    also][]{losch10:_mitsim}. This conveniently avoids dealing with
794  further metric terms, as these are ``hidden'' in the differential cell  further metric terms, as these are ``hidden'' in the differential cell
795  widths. For the $u$-equation ($\alpha=1$) we have:  widths. For the $u$-equation ($\alpha=1$) we have:
796  \begin{align}  \begin{align}
# Line 633  $\sigma_{21}^{Z}=0$, or equivalently $\e Line 891  $\sigma_{21}^{Z}=0$, or equivalently $\e
891    
892  \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\  \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\
893  %  %
894    \noindent\textbf{NOTE: THIS SECTION IS TERRIBLY OUT OF DATE}\\
895  In its original formulation the sea ice model \citep{menemenlis05}  In its original formulation the sea ice model \citep{menemenlis05}
896  uses simple thermodynamics following the appendix of  uses simple thermodynamics following the appendix of
897  \citet{sem76}. This formulation does not allow storage of heat,  \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 \
907  The conductive heat flux depends strongly on the ice thickness $h$.  The conductive heat flux depends strongly on the ice thickness $h$.
908  However, the ice thickness in the model represents a mean over a  However, the ice thickness in the model represents a mean over a
909  potentially very heterogeneous thickness distribution.  In order to  potentially very heterogeneous thickness distribution.  In order to
910  parameterize a sub-grid scale distribution for heat flux  parameterize a sub-grid scale distribution for heat flux computations,
911  computations, the mean ice thickness $h$ is split into seven thickness  the mean ice thickness $h$ is split into $N$ thickness categories
912  categories $H_{n}$ that are equally distributed between $2h$ and a  $H_{n}$ that are equally distributed between $2h$ and a minimum
913  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$
914  \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
915  thickness category is area-averaged to give the total heat flux  is area-averaged to give the total heat flux \citep{hibler84}. To use
916  \citep{hibler84}. To use this thickness category parameterization set  this thickness category parameterization set \code{SEAICE\_multDim} to
917  \code{\#define SEAICE\_MULTICATEGORY}; note that this requires  the number of desired categories (7 is a good guess, for anything
918  different restart files and switching this flag on in the middle of an  larger than 7 modify \code{SEAICE\_SIZE.h}) in
919  integration is not possible.  \code{data.seaice}; note that this requires different restart files
920    and switching this flag on in the middle of an integration is not
921    advised. In order to include the same distribution for snow, set
922    \code{SEAICE\_useMultDimSnow = .TRUE.}; only then, the
923    parameterization of always having a fraction of thin ice is efficient
924    and generally thicker ice is produced \citep{castro-morales14}.
925    
926    
927  The atmospheric heat flux is balanced by an oceanic heat flux from  The atmospheric heat flux is balanced by an oceanic heat flux from
928  below.  The oceanic flux is proportional to  below.  The oceanic flux is proportional to

Legend:
Removed from v.1.18  
changed lines
  Added in v.1.25

  ViewVC Help
Powered by ViewVC 1.1.22