/[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.8 by mlosch, Wed May 13 12:54:45 2009 UTC revision 1.22 by mlosch, Wed Jan 21 18:18:22 2015 UTC
# Line 16  Patrick Heimbach, Chris Hill and Jinlun Line 16  Patrick Heimbach, Chris Hill and Jinlun
16    
17  %----------------------------------------------------------------------  %----------------------------------------------------------------------
18  \subsubsection{Introduction  \subsubsection{Introduction
19  \label{sec:pkg:exf:intro}}  \label{sec:pkg:seaice:intro}}
20    
21    
22  Package ``seaice'' provides a dynamic and thermodynamic interactive  Package ``seaice'' provides a dynamic and thermodynamic interactive
# Line 25  sea-ice model. Line 25  sea-ice model.
25  CPP options enable or disable different aspects of the package  CPP options enable or disable different aspects of the package
26  (Section \ref{sec:pkg:seaice:config}).  (Section \ref{sec:pkg:seaice:config}).
27  Run-Time options, flags, filenames and field-related dates/times are  Run-Time options, flags, filenames and field-related dates/times are
28  set in \texttt{data.seaice}  set in \code{data.seaice}
29  (Section \ref{sec:pkg:seaice:runtime}).  (Section \ref{sec:pkg:seaice:runtime}).
30  A description of key subroutines is given in Section  A description of key subroutines is given in Section
31  \ref{sec:pkg:seaice:subroutines}.  \ref{sec:pkg:seaice:subroutines}.
32  Input fields, units and sign conventions are summarized in  Input fields, units and sign conventions are summarized in
33  Section \ref{sec:pkg:seaice:fields_units}, and available diagnostics  Section \ref{sec:pkg:seaice:fields_units}, and available diagnostics
34  output is listed in Section \ref{sec:pkg:seaice:fields_diagnostics}.  output is listed in Section \ref{sec:pkg:seaice:diagnostics}.
35    
36  %----------------------------------------------------------------------  %----------------------------------------------------------------------
37    
# Line 46  As with all MITgcm packages, SEAICE can Line 46  As with all MITgcm packages, SEAICE can
46  \begin{itemize}  \begin{itemize}
47  %  %
48  \item  \item
49  using the \texttt{packages.conf} file by adding \texttt{seaice} to it,  using the \code{packages.conf} file by adding \code{seaice} to it,
50  %  %
51  \item  \item
52  or using \texttt{genmake2} adding  or using \code{genmake2} adding
53  \texttt{-enable=seaice} or \texttt{-disable=seaice} switches  \code{-enable=seaice} or \code{-disable=seaice} switches
54  %  %
55  \item  \item
56  \textit{required packages and CPP options}: \\  \textit{required packages and CPP options}: \\
57  SEAICE requires the external forcing package \texttt{exf} to be enabled;  SEAICE requires the external forcing package \code{exf} to be enabled;
58  no additional CPP options are required.  no additional CPP options are required.
59  %  %
60  \end{itemize}  \end{itemize}
61  (see Section \ref{sect: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  \texttt{SEAICE\_OPTIONS.h} or in \texttt{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.
67    
68  \begin{table}[h!]  \begin{table}[!ht]
69  \centering  \centering
70    \label{tab:pkg:seaice:cpp}    \label{tab:pkg:seaice:cpp}
71    {\footnotesize    {\footnotesize
# Line 73  Table \ref{tab:pkg:seaice:cpp} summarize Line 73  Table \ref{tab:pkg:seaice:cpp} summarize
73        \hline        \hline
74        \textbf{CPP option}  &  \textbf{Description}  \\        \textbf{CPP option}  &  \textbf{Description}  \\
75        \hline \hline        \hline \hline
76          \texttt{SEAICE\_DEBUG} &          \code{SEAICE\_DEBUG} &
77            Enhance STDOUT for debugging \\            Enhance STDOUT for debugging \\
78          \texttt{SEAICE\_ALLOW\_DYNAMICS} &          \code{SEAICE\_ALLOW\_DYNAMICS} &
79            sea-ice dynamics code \\            sea-ice dynamics code \\
80          \texttt{SEAICE\_CGRID} &          \code{SEAICE\_CGRID} &
81            LSR solver on C-grid (rather than original B-grid) \\            LSR solver on C-grid (rather than original B-grid) \\
82          \texttt{SEAICE\_ALLOW\_EVP} &          \code{SEAICE\_ALLOW\_EVP} &
83            use EVP rather than LSR rheology solver \\            enable use of EVP rheology solver \\
84          \texttt{SEAICE\_EXTERNAL\_FLUXES} &          \code{SEAICE\_ALLOW\_JFNK} &
85              enable use of JFNK rheology solver \\
86            \code{SEAICE\_EXTERNAL\_FLUXES} &
87            use EXF-computed fluxes as starting point \\            use EXF-computed fluxes as starting point \\
88          \texttt{SEAICE\_MULTICATEGORY} &          \code{SEAICE\_ZETA\_SMOOTHREG} &
89            enable 8-category thermodynamics (by default undefined)\\            use differentialable regularization for viscosities \\
90          \texttt{SEAICE\_VARIABLE\_FREEZING\_POINT} &          \code{SEAICE\_VARIABLE\_FREEZING\_POINT} &
91            enable linear dependence of the freezing point on salinity            enable linear dependence of the freezing point on salinity
92            (by default undefined)\\            (by default undefined)\\
93          \texttt{ALLOW\_SEAICE\_FLOODING} &          \code{ALLOW\_SEAICE\_FLOODING} &
94            enable snow to ice conversion for submerged sea-ice \\            enable snow to ice conversion for submerged sea-ice \\
95          \texttt{SEAICE\_SALINITY} &          \code{SEAICE\_VARIABLE\_SALINITY} &
96            enable "salty" sea-ice (by default undefined) \\            enable sea-ice with variable salinity (by default undefined) \\
97          \texttt{SEAICE\_AGE} &          \code{SEAICE\_SITRACER} &
98            enable "age tracer" sea-ice (by default undefined) \\            enable sea-ice tracer package (by default undefined) \\
99          \texttt{SEAICE\_CAP\_HEFF} &          \code{SEAICE\_BICE\_STRESS} &
           enable capping of sea-ice thickness to MAX\_HEFF \\ \hline  
         \texttt{SEAICE\_BICE\_STRESS} &  
100            B-grid only for backward compatiblity: turn on ice-stress on            B-grid only for backward compatiblity: turn on ice-stress on
101            ocean\\            ocean\\
102          \texttt{EXPLICIT\_SSH\_SLOPE} &          \code{EXPLICIT\_SSH\_SLOPE} &
103            B-grid only for backward compatiblity: use ETAN for tilt            B-grid only for backward compatiblity: use ETAN for tilt
104            computations rather than geostrophic velocities \\            computations rather than geostrophic velocities \\
105        \hline        \hline
106      \end{tabular}      \end{tabular}
107    }    }
108    \caption{~}    \caption{Some of the most relevant CPP preprocessor flags in the
109        \code{seaice}-package.}
110  \end{table}  \end{table}
111    
112  %----------------------------------------------------------------------  %----------------------------------------------------------------------
# Line 113  Table \ref{tab:pkg:seaice:cpp} summarize Line 114  Table \ref{tab:pkg:seaice:cpp} summarize
114  \subsubsection{Run-time parameters  \subsubsection{Run-time parameters
115  \label{sec:pkg:seaice:runtime}}  \label{sec:pkg:seaice:runtime}}
116    
117  Run-time parameters are set in files  Run-time parameters (see Table~\ref{tab:pkg:seaice:runtimeparms}) are set
118  \texttt{data.pkg} (read in \texttt{packages\_readparms.F}),  in files \code{data.pkg} (read in \code{packages\_readparms.F}), and
119  and \texttt{data.seaice} (read in \texttt{seaice\_readparms.F}).  \code{data.seaice} (read in \code{seaice\_readparms.F}).
120    
121  \paragraph{Enabling the package}  \paragraph{Enabling the package}
122  ~ \\  ~ \\
123  %  %
124  A package is switched on/off at run-time by setting  A package is switched on/off at run-time by setting
125  (e.g. for SEAICE) \texttt{useSEAICE = .TRUE.} in \texttt{data.pkg}.  (e.g. for SEAICE) \code{useSEAICE = .TRUE.} in \code{data.pkg}.
126    
127  \paragraph{General flags and parameters}  \paragraph{General flags and parameters}
128  ~ \\  ~ \\
129  %  %
130  Table~\ref{tab:pkg:seaice:runtimeparms} lists most run-time parameters.  Table~\ref{tab:pkg:seaice:runtimeparms} lists most run-time parameters.
131  \input{part6/seaice-parms.tex}  \input{s_phys_pkgs/text/seaice-parms.tex}
   
132    
133    \paragraph{Input fields and units\label{sec:pkg:seaice:fields_units}}
134    \begin{description}
135    \item[\code{HeffFile}:] Initial sea ice thickness averaged over grid cell
136      in meters; initializes variable \code{HEFF};
137    \item[\code{AreaFile}:] Initial fractional sea ice cover, range $[0,1]$;
138      initializes variable \code{AREA};
139    \item[\code{HsnowFile}:] Initial snow thickness on sea ice averaged
140      over grid cell in meters; initializes variable \code{HSNOW};
141    \item[\code{HsaltFile}:] Initial salinity of sea ice averaged over grid
142      cell in g/m$^2$; initializes variable \code{HSALT};
143    \end{description}
144    
145  %----------------------------------------------------------------------  %----------------------------------------------------------------------
146  \subsubsection{Description  \subsubsection{Description
# Line 170  viscous-plastic (VP) dynamic-thermodynam Line 181  viscous-plastic (VP) dynamic-thermodynam
181  first introduced by \citet{hib79, hib80}. In order to adapt this model  first introduced by \citet{hib79, hib80}. In order to adapt this model
182  to the requirements of coupled ice-ocean state estimation, many  to the requirements of coupled ice-ocean state estimation, many
183  important aspects of the original code have been modified and  important aspects of the original code have been modified and
184  improved:  improved \citep{losch10:_mitsim}:
185  \begin{itemize}  \begin{itemize}
186  \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
187    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
188    and free-slip lateral boundary conditions;    and free-slip lateral boundary conditions;
189  \item two different solution methods for solving the nonlinear  \item three different solution methods for solving the nonlinear
190    momentum equations have been adopted: LSOR \citep{zhang97}, and EVP    momentum equations have been adopted: LSOR \citep{zhang97}, EVP
191    \citep{hun97};    \citep{hun97}, JFNK \citep{lemieux10,losch14:_jfnk};
192  \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
193    \citet{cam08};    \citet{cam08};
194  \item ice variables are advected by sophisticated, conservative  \item ice variables are advected by sophisticated, conservative
# Line 202  relatively simpler formulation, compared Line 213  relatively simpler formulation, compared
213  to use the VP model as the default dynamic component of our ice  to use the VP model as the default dynamic component of our ice
214  model. To do this we extended the line successive over relaxation  model. To do this we extended the line successive over relaxation
215  (LSOR) method of \citet{zhang97} for use in a parallel  (LSOR) method of \citet{zhang97} for use in a parallel
216  configuration.  configuration. An EVP model and a free-drift implemtation can be
217    selected with runtime flags.
218    
219  Note, that by default the seaice-package includes the orginial  \paragraph{Compatibility with ice-thermodynamics package \code{thsice}\label{sec:pkg:seaice:thsice}}~\\
220    %
221    Note, that by default the \code{seaice}-package includes the orginial
222  so-called zero-layer thermodynamics following \citet{hib80} with a  so-called zero-layer thermodynamics following \citet{hib80} with a
223  snow cover as in \citet{zha98a}. The zero-layer thermodynamic model  snow cover as in \citet{zha98a}. The zero-layer thermodynamic model
224  assumes that ice does not store heat and, therefore, tends to  assumes that ice does not store heat and, therefore, tends to
225  exaggerate the seasonal variability in ice thickness.  This  exaggerate the seasonal variability in ice thickness.  This
226  exaggeration can be significantly reduced by using  exaggeration can be significantly reduced by using
227  \citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic model  \citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic
228  that permits heat storage in ice.  Recently, the three-layer  model that permits heat storage in ice. Recently, the three-layer thermodynamic model has been reformulated by
229  thermodynamic model has been reformulated by \citet{win00}.  The  \citet{win00}.  The reformulation improves model physics by
230  reformulation improves model physics by representing the brine content  representing the brine content of the upper ice with a variable heat
231  of the upper ice with a variable heat capacity.  It also improves  capacity.  It also improves model numerics and consumes less computer
232  model numerics and consumes less computer time and memory.  The Winton  time and memory.
233  sea-ice thermodynamics have been ported to the MIT GCM; they currently  
234  reside under pkg/thsice. The package pkg/thsice is fully compatible  The Winton sea-ice thermodynamics have been ported to the MIT GCM;
235  with pkg/seaice and with pkg/exf. When turned on together with  they currently reside under \code{pkg/thsice}.  The package
236  pkg/seaice, the zero-layer thermodynamics are replaced by the Winton  \code{thsice} is described in section~\ref{sec:pkg:thsice}; it is
237  thermodynamics.  fully compatible with the packages \code{seaice} and \code{exf}. When
238    turned on together with \code{seaice}, the zero-layer thermodynamics
239    are replaced by the Winton thermodynamics. In order to use the
240    \code{seaice}-package with the thermodynamics of \code{thsice},
241    compile both packages and turn both package on in \code{data.pkg}; see
242    an example in \code{global\_ocean.cs32x15/input.icedyn}. Note, that
243    once \code{thsice} is turned on, the variables and diagnostics
244    associated to the default thermodynamics are meaningless, and the
245    diagnostics of \code{thsice} have to be used instead.
246    
247    \paragraph{Surface forcing\label{sec:pkg:seaice:surfaceforcing}}~\\
248    %
249  The sea ice model requires the following input fields: 10-m winds, 2-m  The sea ice model requires the following input fields: 10-m winds, 2-m
250  air temperature and specific humidity, downward longwave and shortwave  air temperature and specific humidity, downward longwave and shortwave
251  radiations, precipitation, evaporation, and river and glacier runoff.  radiations, precipitation, evaporation, and river and glacier runoff.
# Line 232  surface heat flux, and net shortwave flu Line 256  surface heat flux, and net shortwave flu
256  global: in ice-free regions bulk formulae are used to estimate oceanic  global: in ice-free regions bulk formulae are used to estimate oceanic
257  forcing from the atmospheric fields.  forcing from the atmospheric fields.
258    
259  \subsubsection{Dynamics}  \paragraph{Dynamics\label{sec:pkg:seaice:dynamics}}~\\
260  \label{sec:pkg:seaice:dynamics}  %
   
261  \newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}}  \newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}}
262  \newcommand{\vtau}{\vek{\mathbf{\tau}}}  \newcommand{\vtau}{\vek{\mathbf{\tau}}}
263  The momentum equation of the sea-ice model is  The momentum equation of the sea-ice model is
# Line 264  terms are given by Line 287  terms are given by
287    \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|    \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|
288                     R_{air}  (\vek{U}_{air}  -\vek{u}), \\                     R_{air}  (\vek{U}_{air}  -\vek{u}), \\
289    \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}|    \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}|
290                     R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\                     R_{ocean}(\vek{U}_{ocean}-\vek{u}),
291  \end{align*}  \end{align*}
292  where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere  where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
293  and surface currents of the ocean, respectively; $C_{air/ocean}$ are  and surface currents of the ocean, respectively; $C_{air/ocean}$ are
# Line 272  air and ocean drag coefficients; $\rho_{ Line 295  air and ocean drag coefficients; $\rho_{
295  densities; and $R_{air/ocean}$ are rotation matrices that act on the  densities; and $R_{air/ocean}$ are rotation matrices that act on the
296  wind/current vectors.  wind/current vectors.
297    
298    \paragraph{Viscous-Plastic (VP) Rheology\label{sec:pkg:seaice:VPrheology}}~\\
299    %
300  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
301  be related to the ice strain rate and strength by a nonlinear  be related to the ice strain rate and strength by a nonlinear
302  viscous-plastic (VP) constitutive law \citep{hib79, zhang97}:  viscous-plastic (VP) constitutive law \citep{hib79, zhang97}:
# Line 291  The ice strain rate is given by Line 316  The ice strain rate is given by
316  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
317  both thickness $h$ and compactness (concentration) $c$:  both thickness $h$ and compactness (concentration) $c$:
318  \begin{equation}  \begin{equation}
319    P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},    P_{\max} = P^{*}c\,h\,\exp\{-C^{*}\cdot(1-c)\},
320  \label{eq:icestrength}  \label{eq:icestrength}
321  \end{equation}  \end{equation}
322  with the constants $P^{*}$ (run-time parameter \texttt{SEAICE\_strength}) and  with the constants $P^{*}$ (run-time parameter \code{SEAICE\_strength}) and
323  $C^{*}=20$. The nonlinear bulk and shear  $C^{*}=20$. The nonlinear bulk and shear
324  viscosities $\eta$ and $\zeta$ are functions of ice strain rate  viscosities $\eta$ and $\zeta$ are functions of ice strain rate
325  invariants and ice strength such that the principal components of the  invariants and ice strength such that the principal components of the
# Line 313  minor axis $e$ equal to $2$; they are gi Line 338  minor axis $e$ equal to $2$; they are gi
338  \end{align*}  \end{align*}
339  The bulk viscosities are bounded above by imposing both a minimum  The bulk viscosities are bounded above by imposing both a minimum
340  $\Delta_{\min}$ (for numerical reasons, run-time parameter  $\Delta_{\min}$ (for numerical reasons, run-time parameter
341  \texttt{SEAICE\_EPS} with a default value of  \code{SEAICE\_EPS} with a default value of
342  $10^{-10}\text{\,s}^{-1}$) and a maximum $\zeta_{\max} =  $10^{-10}\text{\,s}^{-1}$) and a maximum $\zeta_{\max} =
343  P_{\max}/\Delta^*$, where  P_{\max}/\Delta^*$, where
344  $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. (There is also  $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. (There is also
345  the option of bounding $\zeta$ from below by setting run-time  the option of bounding $\zeta$ from below by setting run-time
346  parameter \texttt{SEAICE\_zetaMin} $>0$, but this is generally not  parameter \code{SEAICE\_zetaMin} $>0$, but this is generally not
347  recommended). For stress tensor computation the replacement pressure $P  recommended). For stress tensor computation the replacement pressure $P
348  = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state  = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state
349  always lies on the elliptic yield curve by definition.  always lies on the elliptic yield curve by definition.
350    
351  In the so-called truncated ellipse method the shear viscosity $\eta$  Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in
352  is capped to suppress any tensile stress \citep{hibler97, geiger98}:  \code{SEAICE\_OPTIONS.h} before compiling replaces the method for
353    bounding $\zeta$ by a smooth (differentiable) expression:
354  \begin{equation}  \begin{equation}
355    \label{eq:etatem}    \label{eq:zetaregsmooth}
356    \eta = \min\left(\frac{\zeta}{e^2},    \begin{split}
357    \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}    \zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min})
358    {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2        \,\zeta_{\max}}\right)\\
359        +4\dot{\epsilon}_{12}^2}}\right).    &= \frac{P}{2\Delta^*}
360      \tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right)
361      \end{split}
362  \end{equation}  \end{equation}
363  To enable this method, set \texttt{\#define SEAICE\_ALLOW\_TEM} in  where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions
364  \texttt{SEAICE\_OPTIONS.h} and turn it on with  by zero.
 \texttt{SEAICEuseTEM=.TRUE.} in \texttt{data.seaice}.  
   
 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.  
365    
366    \paragraph{LSR and  JFNK solver \label{sec:pkg:seaice:LSRJFNK}}~\\
367    %
368    % By default, the VP-model is integrated by a Pickwith the
369    % semi-implicit line successive over relaxation (LSOR)-solver of
370    % \citet{zhang97}, which allows for long time steps that, in our case,
371    % are limited by the explicit treatment of the Coriolis term. The
372    % explicit treatment of the Coriolis term does not represent a severe
373    % limitation because it restricts the time step to approximately the
374    % same length as in the ocean model where the Coriolis term is also
375    % treated explicitly.
376    
377    \newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}}
378    %
379    In the matrix notation, the discretized momentum equations can be
380    written as
381    \begin{equation}
382      \label{eq:matrixmom}
383      \mat{A}(\vek{x})\,\vek{x} = \vek{b}(\vek{x}).
384    \end{equation}
385    The solution vector $\vek{x}$ consists of the two velocity components
386    $u$ and $v$ that contain the velocity variables at all grid points and
387    at one time level. The standard (and default) method for solving
388    Eq.\,(\ref{eq:matrixmom}) in the sea ice component of the
389    \mbox{MITgcm}, as in many sea ice models, is an iterative Picard
390    solver: in the $k$-th iteration a linearized form
391    $\mat{A}(\vek{x}^{k-1})\,\vek{x}^{k} = \vek{b}(\vek{x}^{k-1})$ is
392    solved (in the case of the MITgcm it is a Line Successive (over)
393    Relaxation (LSR) algorithm \citep{zhang97}).  Picard solvers converge
394    slowly, but generally the iteration is terminated after only a few
395    non-linear steps \citep{zhang97, lemieux09} and the calculation
396    continues with the next time level. This method is the default method
397    in the MITgcm. The number of non-linear iteration steps or pseudo-time
398    steps can be controlled by the runtime parameter
399    \code{NPSEUDOTIMESTEPS} (default is 2).
400    
401    In order to overcome the poor convergence of the Picard-solver,
402    \citet{lemieux10} introduced a Jacobian-free Newton-Krylov solver for
403    the sea ice momentum equations. This solver is also implemented in the
404    MITgcm \citep{losch14:_jfnk}. The Newton method transforms minimizing
405    the residual $\vek{F}(\vek{x}) = \mat{A}(\vek{x})\,\vek{x} -
406    \vek{b}(\vek{x})$ to finding the roots of a multivariate Taylor
407    expansion of the residual \vek{F} around the previous ($k-1$) estimate
408    $\vek{x}^{k-1}$:
409    \begin{equation}
410      \label{eq:jfnktaylor}
411      \vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k}) =
412      \vek{F}(\vek{x}^{k-1}) + \vek{F}'(\vek{x}^{k-1})\,\delta\vek{x}^{k}
413    \end{equation}
414    with the Jacobian $\mat{J}\equiv\vek{F}'$. The root
415    $\vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k})=0$ is found by solving
416    \begin{equation}
417      \label{eq:jfnklin}
418      \mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k} = -\vek{F}(\vek{x}^{k-1})
419    \end{equation}
420    for $\delta\vek{x}^{k}$. The next ($k$-th) estimate is given by
421    $\vek{x}^{k}=\vek{x}^{k-1}+a\,\delta\vek{x}^{k}$. In order to avoid
422    overshoots the factor $a$ is iteratively reduced in a line search
423    ($a=1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \ldots$) until
424    $\|\vek{F}(\vek{x}^k)\| < \|\vek{F}(\vek{x}^{k-1})\|$, where
425    $\|\cdot\|=\int\cdot\,dx^2$ is the $L_2$-norm. In practice, the line
426    search is stopped at $a=\frac{1}{8}$. The line search starts after
427    $\code{SEAICE\_JFNK\_lsIter}$ non-linear Newton iterations (off by
428    default).
429    
430    
431    Forming the Jacobian $\mat{J}$ explicitly is often avoided as ``too
432    error prone and time consuming'' \citep{knoll04:_jfnk}. Instead,
433    Krylov methods only require the action of \mat{J} on an arbitrary
434    vector \vek{w} and hence allow a matrix free algorithm for solving
435    Eq.\,(\ref{eq:jfnklin}) \citep{knoll04:_jfnk}. The action of \mat{J}
436    can be approximated by a first-order Taylor series expansion:
437    \begin{equation}
438      \label{eq:jfnkjacvecfd}
439      \mat{J}(\vek{x}^{k-1})\,\vek{w} \approx
440      \frac{\vek{F}(\vek{x}^{k-1}+\epsilon\vek{w}) - \vek{F}(\vek{x}^{k-1})}
441      {\epsilon}
442    \end{equation}
443    or computed exactly with the help of automatic differentiation (AD)
444    tools. \code{SEAICE\_JFNKepsilon} sets the step size
445    $\epsilon$.
446    
447    We use the Flexible Generalized Minimum RESidual method
448    \citep[FGMRES,][]{saad93:_fgmres} with right-hand side preconditioning
449    to solve Eq.\,(\ref{eq:jfnklin}) iteratively starting from a first
450    guess of $\delta\vek{x}^{k}_{0} = 0$. For the preconditioning matrix
451    \mat{P} we choose a simplified form of the system matrix
452    $\mat{A}(\vek{x}^{k-1})$ \citep{lemieux10} where $\vek{x}^{k-1}$ is
453    the estimate of the previous Newton step $k-1$. The transformed
454    equation\,(\ref{eq:jfnklin}) becomes
455    \begin{equation}
456      \label{eq:jfnklinpc}
457      \mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z} =
458      -\vek{F}(\vek{x}^{k-1}),
459      \quad\text{with}\quad \delta\vek{z}=\mat{P}\delta\vek{x}^{k}.
460    \end{equation}
461    The Krylov method iteratively improves the approximate solution
462    to~(\ref{eq:jfnklinpc}) in subspace ($\vek{r}_0$,
463    $\mat{J}\mat{P}^{-1}\vek{r}_0$, $(\mat{J}\mat{P}^{-1})^2\vek{r}_0$,
464    \ldots, $(\mat{J}\mat{P}^{-1})^m\vek{r}_0$) with increasing $m$;
465    $\vek{r}_0 = -\vek{F}(\vek{x}^{k-1})
466    -\mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k}_{0}$
467    %-\vek{F}(\vek{x}^{k-1})
468    %-\mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z}$
469    is the initial residual of
470    (\ref{eq:jfnklin}); $\vek{r}_0=-\vek{F}(\vek{x}^{k-1})$ with the first
471    guess $\delta\vek{x}^{k}_{0}=0$. We allow a Krylov-subspace of
472    dimension~$m=50$ and we do not use restarts. The preconditioning operation
473    involves applying $\mat{P}^{-1}$ to the basis vectors $\vek{v}_0,
474    \vek{v}_1, \vek{v}_2, \ldots, \vek{v}_m$ of the Krylov subspace. This
475    operation is approximated by solving the linear system
476    $\mat{P}\,\vek{w}=\vek{v}_i$. Because $\mat{P} \approx
477    \mat{A}(\vek{x}^{k-1})$, we can use the LSR-algorithm \citep{zhang97}
478    already implemented in the Picard solver. Each preconditioning
479    operation uses a fixed number of 10~LSR-iterations avoiding any
480    termination criterion. More details and results can be found in
481    \citet{lemieux10, losch14:_jfnk}.
482    
483    To use the JFNK-solver set \code{SEAICEuseJFNK = .TRUE.} in the
484    namelist file \code{data.seaice}; \code{SEAICE\_ALLOW\_JFNK} needs to
485    be defined in \code{SEAICE\_OPTIONS.h} and we recommend using a smooth
486    regularization of $\zeta$ by defining \code{SEAICE\_ZETA\_SMOOTHREG}
487    (see above) for better convergence.  The non-linear Newton iteration
488    is terminated when the $L_2$-norm of the residual is reduced by
489    $\gamma_{\mathrm{nl}}$ (runtime parameter \code{JFNKgamma\_nonlin =
490      1.e-4} will already lead to expensive simulations) with respect to
491    the initial norm: $\|\vek{F}(\vek{x}^k)\| <
492    \gamma_{\mathrm{nl}}\|\vek{F}(\vek{x}^0)\|$.  Within a non-linear
493    iteration, the linear FGMRES solver is terminated when the residual is
494    smaller than $\gamma_k\|\vek{F}(\vek{x}^{k-1})\|$ where $\gamma_k$ is
495    determined by
496    \begin{equation}
497      \label{eq:jfnkgammalin}
498      \gamma_k =
499      \begin{cases}
500        \gamma_0 &\text{for $\|\vek{F}(\vek{x}^{k-1})\| \geq r$},  \\
501        \max\left(\gamma_{\min},
502        \frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)  
503    %    \phi\left(\frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)^\alpha\right)  
504        &\text{for $\|\vek{F}(\vek{x}^{k-1})\| < r$,}
505      \end{cases}
506    \end{equation}
507    so that the linear tolerance parameter $\gamma_k$ decreases with the
508    non-linear Newton step as the non-linear solution is approached. This
509    inexact Newton method is generally more robust and computationally
510    more efficient than exact methods \citep[e.g.,][]{knoll04:_jfnk}.
511    % \footnote{The general idea behind
512    %   inexact Newton methods is this: The Krylov solver is ``only''
513    %   used to find an intermediate solution of the linear
514    %   equation\,(\ref{eq:jfnklin}) that is used to improve the approximation of
515    %   the actual equation\,(\ref{eq:matrixmom}). With the choice of a
516    %   relatively weak lower limit for FGMRES convergence
517    %   $\gamma_{\min}$ we make sure that the time spent in the FGMRES
518    %   solver is reduced at the cost of more Newton iterations. Newton
519    %   iterations are cheaper than Krylov iterations so that this choice
520    %   improves the overall efficiency.}
521    Typical parameter choices are
522    $\gamma_0=\code{JFNKgamma\_lin\_max}=0.99$,
523    $\gamma_{\min}=\code{JFNKgamma\_lin\_min}=0.1$, and $r =
524    \code{JFNKres\_tFac}\times\|\vek{F}(\vek{x}^{0})\|$ with
525    $\code{JFNKres\_tFac} = \frac{1}{2}$. We recommend a maximum number of
526    non-linear iterations $\code{SEAICEnewtonIterMax} = 100$ and a maximum
527    number of Krylov iterations $\code{SEAICEkrylovIterMax} = 50$, because
528    the Krylov subspace has a fixed dimension of 50.
529    
530    \paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\
531    %
532  \citet{hun97}'s introduced an elastic contribution to the strain  \citet{hun97}'s introduced an elastic contribution to the strain
533  rate in order to regularize Eq.~\ref{eq:vpequation} in such a way that  rate in order to regularize Eq.~\ref{eq:vpequation} in such a way that
534  the resulting elastic-viscous-plastic (EVP) and VP models are  the resulting elastic-viscous-plastic (EVP) and VP models are
# Line 362  identical at steady state, Line 546  identical at steady state,
546  %used and compared the present sea-ice model study.  %used and compared the present sea-ice model study.
547  The EVP-model uses an explicit time stepping scheme with a short  The EVP-model uses an explicit time stepping scheme with a short
548  timestep. According to the recommendation of \citet{hun97}, the  timestep. According to the recommendation of \citet{hun97}, the
549  EVP-model is stepped forward in time 120 times within the physical  EVP-model should be stepped forward in time 120 times
550  ocean model time step (although this parameter is under debate), to  ($\code{SEAICE\_deltaTevp} = \code{SEAICIE\_deltaTdyn}/120$) within
551  allow for elastic waves to disappear.  Because the scheme does not  the physical ocean model time step (although this parameter is under
552  require a matrix inversion it is fast in spite of the small internal  debate), to allow for elastic waves to disappear.  Because the scheme
553  timestep and simple to implement on parallel computers  does not require a matrix inversion it is fast in spite of the small
554    internal timestep and simple to implement on parallel computers
555  \citep{hun97}. For completeness, we repeat the equations for the  \citep{hun97}. For completeness, we repeat the equations for the
556  components of the stress tensor $\sigma_{1} =  components of the stress tensor $\sigma_{1} =
557  \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and  \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
# Line 387  abbreviations, the equations~\ref{eq:evp Line 572  abbreviations, the equations~\ref{eq:evp
572    \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}    \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
573    &= \frac{P}{4T\Delta} D_S    &= \frac{P}{4T\Delta} D_S
574  \end{align}  \end{align}
575  Here, the elastic parameter $E$ is redefined in terms of a damping timescale  Here, the elastic parameter $E$ is redefined in terms of a damping
576  $T$ for elastic waves \[E=\frac{\zeta}{T}.\]  timescale $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
577  $T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and  $T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and the external
578  the external (long) timestep $\Delta{t}$. \citet{hun97} recommend  (long) timestep $\Delta{t}$.  $E_{0} = \frac{1}{3}$ is the default
579  $E_{0} = \frac{1}{3}$ (which is the default value in the code).  value in the code and close to what \citet{hun97} and
580    \citet{hun01} recommend.
581    
582  To use the EVP solver, make sure that both \texttt{SEAICE\_CGRID} and  To use the EVP solver, make sure that both \code{SEAICE\_CGRID} and
583  \texttt{SEAICE\_ALLOW\_EVP} are defined in \texttt{SEAICE\_OPTIONS.h}  \code{SEAICE\_ALLOW\_EVP} are defined in \code{SEAICE\_OPTIONS.h}
584  (default). The solver is turned on by setting the sub-cycling time  (default). The solver is turned on by setting the sub-cycling time
585  step \texttt{SEAICE\_deltaTevp} to a value larger than zero. The  step \code{SEAICE\_deltaTevp} to a value larger than zero. The
586  choice of this time step is under debate. \citet{hun97} recommend  choice of this time step is under debate. \citet{hun97} recommend
587  order(120) time steps for the EVP solver within one model time step  order(120) time steps for the EVP solver within one model time step
588  $\Delta{t}$ (\texttt{deltaTmom}). One can also choose order(120) time  $\Delta{t}$ (\code{deltaTmom}). One can also choose order(120) time
589  steps within the forcing time scale, but then we recommend adjusting  steps within the forcing time scale, but then we recommend adjusting
590  the damping time scale $T$ accordingly, by setting either  the damping time scale $T$ accordingly, by setting either
591  \texttt{SEAICE\_elasticParm} ($E_{0}$), so that  \code{SEAICE\_elasticParm} ($E_{0}$), so that
592  $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly  $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly
593  \texttt{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.  \code{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.
594    
595    \paragraph{More stable variant of Elastic-Viscous-Plastic Dynamics:  EVP*\label{sec:pkg:seaice:EVPstar}}~\\
596    %
597    The genuine EVP schemes appears to give noisy solutions \citep{hun01,
598      lemieux12, bouillon13}. This has lead to a modified EVP or EVP*
599    \citep{lemieux12, bouillon13, kimmritz15}; here, refer to these
600    variants by EVP*. The main idea is to modify the ``natural''
601    time-discretization of the momentum equations:
602    \begin{equation}
603      \label{eq:evpstar}
604      m\frac{D\vec{u}}{Dt} \approx m\frac{u^{p+1}-u^{n}}{\Delta{t}}
605      + \beta^{*}\frac{u^{p+1}-u^{p}}{\Delta{t}_{\mathrm{EVP}}}
606    \end{equation}
607    where $n$ is the previous time step index, and $p$ is the previous
608    sub-cycling index. The term allows the definition of a residual
609    $|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to
610    $0$ and a re-interpretation of EVP as a pure iterative solver where
611    the sub-cycling has lost all time-relation \citep{bouillon13,
612      kimmritz15}. Using the terminology of \citet{kimmritz15}, the
613    evolution equations of stress $\sigma_{ij}$ and momentum $\vec{u}$ can
614    be written as:
615    \begin{align}
616      \label{eq:evpstarsigma}
617      \sigma_{ij}^{p+1}&=\sigma_{ij}^p+\frac{1}{\alpha}
618      \Big(\sigma_{ij}(\vec{u}^p)-\sigma_{ij}^p\Big),
619      \phantom{\int}\\
620      \label{eq:evpstarmom}
621      \vec{u}^{p+1}&=\vec{u}^p+\frac{1}{\beta}
622      \Big(\frac{\Delta t}{m}\nabla \cdot{\bf \sigma}^{p+1}+
623      \frac{\Delta t}{m}\vec{R}^{p+1/2}+\vec{u}_n-\vec{u}^p\Big).
624    \end{align}
625    $\vec{R}$ contains all terms in the momentum equations except for the
626    rheology terms and the time derivative, $\alpha$ and $\beta$ are free
627    parameters (\code{SEAICE\_evpAlpha}, \code{SEAICE\_evpBeta}) that
628    replace the time stepping parameters \code{SEAICE\_deltaTevp}
629    ($\Delta{T}_{\mathrm{EVP}}$), \code{SEAICE\_elasticParm} ($E_{0}$), or
630    \code{SEAICE\_evpTauRelax} ($T$). $\alpha$ and $\beta$ determine the
631    speed of convergence and the stability. Usually, it makes sense to use
632    $\alpha = \beta$, and \code{SEAICEnEVPstarSteps} $>> \alpha = \beta$
633    \citep{kimmritz15}.
634    
635    In order to use EVP* in the MITgcm, set \code{SEAICEuseEVPstar =
636      .TRUE.,} in \code{data.seaice}. \code{SEAICEuseEVPrev =.TRUE.,} uses
637    the actual form of equations (\ref{eq:evpstarsigma}) and
638    (\ref{eq:evpstarmom}) with fewer implicit terms and the factor of
639    $e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2})
640    and (\ref{eq:evpstresstensor12}).  This turns out to improve
641    convergence \citep{bouillon13}.
642    
643    Note, that for historical reasons, \code{SEAICE\_deltaTevp} needs to
644    be set to some value in order to use also EVP*. Also note, that
645    probably because of the C-grid staggering of velocities and stresses,
646    EVP* does not converge as successfully as in \citet{kimmritz15}.
647    
648    
649    \paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\
650    %
651    In the so-called truncated ellipse method the shear viscosity $\eta$
652    is capped to suppress any tensile stress \citep{hibler97, geiger98}:
653    \begin{equation}
654      \label{eq:etatem}
655      \eta = \min\left(\frac{\zeta}{e^2},
656      \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
657      {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
658          +4\dot{\epsilon}_{12}^2}}\right).
659    \end{equation}
660    To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in
661    \code{SEAICE\_OPTIONS.h} and turn it on with
662    \code{SEAICEuseTEM} in \code{data.seaice}.
663    
664    \paragraph{Ice-Ocean stress \label{sec:pkg:seaice:iceoceanstress}}~\\
665    %
666  Moving sea ice exerts a stress on the ocean which is the opposite of  Moving sea ice exerts a stress on the ocean which is the opposite of
667  the stress $\vtau_{ocean}$ in Eq.~\ref{eq:momseaice}. This stess is  the stress $\vtau_{ocean}$ in Eq.~\ref{eq:momseaice}. This stess is
668  applied directly to the surface layer of the ocean model. An  applied directly to the surface layer of the ocean model. An
# Line 422  advect tracers, is really an average ove Line 679  advect tracers, is really an average ove
679  velocity and the ice velocity leading to an inconsistency as the ice  velocity and the ice velocity leading to an inconsistency as the ice
680  temperature and salinity are different from the oceanic variables.  temperature and salinity are different from the oceanic variables.
681  To turn on the stress formulation of \citet{hibler87}, set  To turn on the stress formulation of \citet{hibler87}, set
682  \texttt{useHB87StressCoupling=.TRUE.} in \texttt{data.seaice}.  \code{useHB87StressCoupling=.TRUE.} in \code{data.seaice}.
683    
684    
685  % Our discretization differs from \citet{zhang97, zhang03} in the  % Our discretization differs from \citet{zhang97, zhang03} in the
# Line 434  To turn on the stress formulation of \ci Line 691  To turn on the stress formulation of \ci
691  % differences and averaging is only involved in computing $\Delta$ and  % differences and averaging is only involved in computing $\Delta$ and
692  % $P$ at vorticity points.  % $P$ at vorticity points.
693    
694  \subsubsection{Finite-volume discretization of the stress tensor  \paragraph{Finite-volume discretization of the stress tensor
695    divergence}    divergence\label{sec:pkg:seaice:discretization}}~\\
696  \label{sec:pkg:seaice:discretization}  %
697  On an Arakawa C~grid, ice thickness and concentration and thus ice  On an Arakawa C~grid, ice thickness and concentration and thus ice
698  strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are  strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are
699  naturally defined a C-points in the center of the grid  naturally defined a C-points in the center of the grid
# Line 498  relation $\sigma_{\alpha\beta} = 2\eta\d Line 755  relation $\sigma_{\alpha\beta} = 2\eta\d
755  [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2  [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2
756  ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence  ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence
757  $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is  $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is
758  discretized in finite volumes. This conveniently avoids dealing with  discretized in finite volumes \citep[see
759    also][]{losch10:_mitsim}. This conveniently avoids dealing with
760  further metric terms, as these are ``hidden'' in the differential cell  further metric terms, as these are ``hidden'' in the differential cell
761  widths. For the $u$-equation ($\alpha=1$) we have:  widths. For the $u$-equation ($\alpha=1$) we have:
762  \begin{align}  \begin{align}
# Line 515  widths. For the $u$-equation ($\alpha=1$ Line 773  widths. For the $u$-equation ($\alpha=1$
773    + \Delta{x}_1\sigma_{21}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}    + \Delta{x}_1\sigma_{21}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
774    \biggr\} \\ \notag    \biggr\} \\ \notag
775    =& \frac{1}{A_{i,j}^w} \biggl\{    =& \frac{1}{A_{i,j}^w} \biggl\{
776    (\Delta{x}_2\sigma_{11})_{i,j}^C - (\Delta{x}_2\sigma_{11})_{i-1,j}^C \\\notag    (\Delta{x}_2\sigma_{11})_{i,j}^C -
777      (\Delta{x}_2\sigma_{11})_{i-1,j}^C
778      \\\notag
779    \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{}    \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{}
780    + (\Delta{x}_1\sigma_{21})_{i,j+1}^Z - (\Delta{x}_1\sigma_{21})_{i,j}^Z    + (\Delta{x}_1\sigma_{21})_{i,j+1}^Z - (\Delta{x}_1\sigma_{21})_{i,j}^Z
781    \biggr\}    \biggr\}
782    \intertext{with}  \end{align}
783    with
784    \begin{align}
785    (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+}    (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+}
786    \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}    \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
787    \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag    \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
# Line 530  widths. For the $u$-equation ($\alpha=1$ Line 792  widths. For the $u$-equation ($\alpha=1$
792    \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}    \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
793    k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag    k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
794    \phantom{=}& - \Delta{y}_{i,j}^{F} \frac{P}{2} \\    \phantom{=}& - \Delta{y}_{i,j}^{F} \frac{P}{2} \\
   %  
795    (\Delta{x}_1\sigma_{21})_{i,j}^Z =& \phantom{+}    (\Delta{x}_1\sigma_{21})_{i,j}^Z =& \phantom{+}
796    \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}    \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
797    \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\ \notag    \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\ \notag
# Line 562  Similarly, we have for the $v$-equation Line 823  Similarly, we have for the $v$-equation
823    \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{}    \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{}
824    + (\Delta{x}_1\sigma_{22})_{i,j}^C - (\Delta{x}_1\sigma_{22})_{i,j-1}^C    + (\Delta{x}_1\sigma_{22})_{i,j}^C - (\Delta{x}_1\sigma_{22})_{i,j-1}^C
825    \biggr\}    \biggr\}
826    \intertext{with}  \end{align}
827    with
828    \begin{align}
829    (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+}    (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+}
830    \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}    \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
831    \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\\notag    \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}}
832    &+ \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}    \\\notag &
833    \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag    + \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
834    &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}    \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\\notag
   k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag  
835    &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}    &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
836      k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2}
837      \\\notag &
838      - \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
839    k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} \\ \notag    k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} \\ \notag
   %  
840    (\Delta{x}_2\sigma_{22})_{i,j}^C =& \phantom{+}    (\Delta{x}_2\sigma_{22})_{i,j}^C =& \phantom{+}
841    \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}    \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
842    \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag    \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
# Line 591  free slip boundary conditions the latera Line 855  free slip boundary conditions the latera
855  analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set  analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set
856  $\sigma_{21}^{Z}=0$, or equivalently $\eta_{i,j}^{Z}=0$, on boundaries.  $\sigma_{21}^{Z}=0$, or equivalently $\eta_{i,j}^{Z}=0$, on boundaries.
857    
858  \subsubsection{Thermodynamics}  \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\
859  \label{sec:pkg:seaice:thermodynamics}  %
   
860  In its original formulation the sea ice model \citep{menemenlis05}  In its original formulation the sea ice model \citep{menemenlis05}
861  uses simple thermodynamics following the appendix of  uses simple thermodynamics following the appendix of
862  \citet{sem76}. This formulation does not allow storage of heat,  \citet{sem76}. This formulation does not allow storage of heat,
# Line 616  minimum imposed ice thickness of $5\text Line 879  minimum imposed ice thickness of $5\text
879  \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each  \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
880  thickness category is area-averaged to give the total heat flux  thickness category is area-averaged to give the total heat flux
881  \citep{hibler84}. To use this thickness category parameterization set  \citep{hibler84}. To use this thickness category parameterization set
882  \texttt{\#define SEAICE\_MULTICATEGORY}; note that this requires  \code{\#define SEAICE\_MULTICATEGORY}; note that this requires
883  different restart files and switching this flag on in the middle of an  different restart files and switching this flag on in the middle of an
884  integration is not possible.  integration is not possible.
885    
# Line 626  $\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ w Line 889  $\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ w
889  the density and heat capacity of sea water and $T_{fr}$ is the local  the density and heat capacity of sea water and $T_{fr}$ is the local
890  freezing point temperature that is a function of salinity. This flux  freezing point temperature that is a function of salinity. This flux
891  is not assumed to instantaneously melt or create ice, but a time scale  is not assumed to instantaneously melt or create ice, but a time scale
892  of three days (run-time parameter \texttt{SEAICE\_gamma\_t}) is used  of three days (run-time parameter \code{SEAICE\_gamma\_t}) is used
893  to relax $T_{w}$ to the freezing point.  to relax $T_{w}$ to the freezing point.
894  %  %
895  The parameterization of lateral and vertical growth of sea ice follows  The parameterization of lateral and vertical growth of sea ice follows
896  that of \citet{hib79, hib80}; the so-called lead closing parameter  that of \citet{hib79, hib80}; the so-called lead closing parameter
897  $h_{0}$ (run-time parameter \texttt{HO}) has a default value of  $h_{0}$ (run-time parameter \code{HO}) has a default value of
898  0.5~meters.  0.5~meters.
899    
900  On top of the ice there is a layer of snow that modifies the heat flux  On top of the ice there is a layer of snow that modifies the heat flux
# Line 644  the snow is flooded, a simple mass conse Line 907  the snow is flooded, a simple mass conse
907  snowice formation (a flood-freeze algorithm following Archimedes'  snowice formation (a flood-freeze algorithm following Archimedes'
908  principle) turns snow into ice until the ice surface is back at $z=0$  principle) turns snow into ice until the ice surface is back at $z=0$
909  \citep{leppaeranta83}. The flood-freeze algorithm is enabled with the CPP-flag  \citep{leppaeranta83}. The flood-freeze algorithm is enabled with the CPP-flag
910  \texttt{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter  \code{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter
911  \texttt{SEAICEuseFlooding=.true.}.  \code{SEAICEuseFlooding=.true.}.
912    
913    \paragraph{Advection of thermodynamic variables\label{sec:pkg:seaice:advection}}~\\
914    %
915  Effective ice thickness (ice volume per unit area,  Effective ice thickness (ice volume per unit area,
916  $c\cdot{h}$), concentration $c$ and effective snow thickness  $c\cdot{h}$), concentration $c$ and effective snow thickness
917  ($c\cdot{h}_{s}$) are advected by ice velocities:  ($c\cdot{h}_{s}$) are advected by ice velocities:
# Line 659  where $\Gamma_X$ are the thermodynamic s Line 924  where $\Gamma_X$ are the thermodynamic s
924  diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.  diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.
925  %  %
926  From the various advection scheme that are available in the MITgcm, we  From the various advection scheme that are available in the MITgcm, we
927  choose flux-limited schemes \citep[multidimensional 2nd and 3rd-order  recommend flux-limited schemes \citep[multidimensional 2nd and
928  advection scheme with flux limiter][]{roe:85, hundsdorfer94} to  3rd-order advection scheme with flux limiter][]{roe:85, hundsdorfer94}
929  preserve sharp gradients and edges that are typical of sea ice  to preserve sharp gradients and edges that are typical of sea ice
930  distributions and to rule out unphysical over- and undershoots  distributions and to rule out unphysical over- and undershoots
931  (negative thickness or concentration). These scheme conserve volume  (negative thickness or concentration). These schemes conserve volume
932  and horizontal area and are unconditionally stable, so that we can set  and horizontal area and are unconditionally stable, so that we can set
933  $D_{X}=0$. Run-timeflags: \texttt{SEAICEadvScheme} (default=2),  $D_{X}=0$. Run-timeflags: \code{SEAICEadvScheme} (default=2, is the
934  \texttt{DIFF1} (default=0.004).  historic 2nd-order, centered difference scheme), \code{DIFF1} =
935    $D_{X}/\Delta{x}$
936  There is considerable doubt about the reliability of a ``zero-layer''  (default=0.004).
937  thermodynamic model --- \citet{semtner84} found significant errors in  
938  phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in  The MITgcm sea ice model provides the option to use
939  such models --- so that today many sea ice models employ more complex  the thermodynamics model of \citet{win00}, which in turn is based on
940  thermodynamics. The MITgcm sea ice model provides the option to use  the 3-layer model of \citet{sem76} and which treats brine content by
941  the thermodynamics model of \citet{win00}, which in turn is based  means of enthalpy conservation; the corresponding package
942  on the 3-layer model of \citet{sem76} and which treats brine  \code{thsice} is described in section~\ref{sec:pkg:thsice}. This
943  content by means of enthalpy conservation. This scheme requires  scheme requires additional state variables, namely the enthalpy of the
944  additional state variables, namely the enthalpy of the two ice layers  two ice layers (instead of effective ice salinity), to be advected by
945  (instead of effective ice salinity), to be advected by ice velocities.  ice velocities.
946  %  %
947  The internal sea ice temperature is inferred from ice enthalpy.  To  The internal sea ice temperature is inferred from ice enthalpy.  To
948  avoid unphysical (negative) values for ice thickness and  avoid unphysical (negative) values for ice thickness and
949  concentration, a positive 2nd-order advection scheme with a SuperBee  concentration, a positive 2nd-order advection scheme with a SuperBee
950  flux limiter \citep{roe:85} is used in this study to advect all  flux limiter \citep{roe:85} should be used to advect all
951  sea-ice-related quantities of the \citet{win00} thermodynamic  sea-ice-related quantities of the \citet{win00} thermodynamic model
952  model.  Because of the non-linearity of the advection scheme, care  (runtime flag \code{thSIceAdvScheme=77} and
953  must be taken in advecting these quantities: when simply using ice  \code{thSIce\_diffK}=$D_{X}$=0 in \code{data.ice}, defaults are 0).  Because of the
954  velocity to advect enthalpy, the total energy (i.e., the volume  non-linearity of the advection scheme, care must be taken in advecting
955  integral of enthalpy) is not conserved. Alternatively, one can advect  these quantities: when simply using ice velocity to advect enthalpy,
956  the energy content (i.e., product of ice-volume and enthalpy) but then  the total energy (i.e., the volume integral of enthalpy) is not
957  false enthalpy extrema can occur, which then leads to unrealistic ice  conserved. Alternatively, one can advect the energy content (i.e.,
958  temperature.  In the currently implemented solution, the sea-ice mass  product of ice-volume and enthalpy) but then false enthalpy extrema
959  flux is used to advect the enthalpy in order to ensure conservation of  can occur, which then leads to unrealistic ice temperature.  In the
960  enthalpy and to prevent false enthalpy extrema.  currently implemented solution, the sea-ice mass flux is used to
961    advect the enthalpy in order to ensure conservation of enthalpy and to
962    prevent false enthalpy extrema. %
963    
964  %----------------------------------------------------------------------  %----------------------------------------------------------------------
965    
966  \subsubsection{Key subroutines  \subsubsection{Key subroutines
967  \label{sec:pkg:seaice:subroutines}}  \label{sec:pkg:seaice:subroutines}}
968    
969  Top-level routine: \texttt{seaice\_model.F}  Top-level routine: \code{seaice\_model.F}
970    
971  {\footnotesize  {\footnotesize
972  \begin{verbatim}  \begin{verbatim}
# Line 758  Diagnostics output is available via the Line 1025  Diagnostics output is available via the
1025  Available output fields are summarized in  Available output fields are summarized in
1026  Table \ref{tab:pkg:seaice:diagnostics}.  Table \ref{tab:pkg:seaice:diagnostics}.
1027    
1028  \begin{table}[h!]  \input{s_phys_pkgs/text/seaice_diags.tex}
 \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}  
   
1029    
1030  %\subsubsection{Package Reference}  %\subsubsection{Package Reference}
1031    
# Line 817  Table \ref{tab:pkg:seaice:diagnostics}. Line 1033  Table \ref{tab:pkg:seaice:diagnostics}.
1033  \label{sec:pkg:seaice:experiments}  \label{sec:pkg:seaice:experiments}
1034    
1035  \begin{itemize}  \begin{itemize}
1036  \item{Labrador Sea experiment in lab\_sea verification directory. }  \item{Labrador Sea experiment in \code{lab\_sea} verification directory. }
1037    \item \code{seaice\_obcs}, based on \code{lab\_sea}
1038    \item \code{offline\_exf\_seaice/input.seaicetd}, based on \code{lab\_sea}
1039    \item \code{global\_ocean.cs32x15/input.icedyn} and
1040      \code{global\_ocean.cs32x15/input.seaice}, global
1041      cubed-sphere-experiment with combinations of \code{seaice} and
1042      \code{thsice}
1043  \end{itemize}  \end{itemize}
1044    
1045    
1046  %%% Local Variables:  %%% Local Variables:
1047  %%% mode: latex  %%% mode: latex
1048  %%% TeX-master: "../manual"  %%% TeX-master: "../../manual"
1049  %%% End:  %%% End:

Legend:
Removed from v.1.8  
changed lines
  Added in v.1.22

  ViewVC Help
Powered by ViewVC 1.1.22