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 |
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\\ |
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 |
%---------------------------------------------------------------------- |
%---------------------------------------------------------------------- |
115 |
\subsubsection{Run-time parameters |
\subsubsection{Run-time parameters |
116 |
\label{sec:pkg:seaice:runtime}} |
\label{sec:pkg:seaice:runtime}} |
117 |
|
|
118 |
Run-time parameters are set in files |
Run-time parameters (see Table~\ref{tab:pkg:seaice:runtimeparms}) are set |
119 |
\code{data.pkg} (read in \code{packages\_readparms.F}), |
in files \code{data.pkg} (read in \code{packages\_readparms.F}), and |
120 |
and \code{data.seaice} (read in \code{seaice\_readparms.F}). |
\code{data.seaice} (read in \code{seaice\_readparms.F}). |
121 |
|
|
122 |
\paragraph{Enabling the package} |
\paragraph{Enabling the package} |
123 |
~ \\ |
~ \\ |
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 |
%---------------------------------------------------------------------- |
%---------------------------------------------------------------------- |
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 |
214 |
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 |
215 |
model. To do this we extended the line successive over relaxation |
model. To do this we extended the line successive over relaxation |
216 |
(LSOR) method of \citet{zhang97} for use in a parallel |
(LSOR) method of \citet{zhang97} for use in a parallel |
217 |
configuration. |
configuration. An EVP model and a free-drift implemtation can be |
218 |
|
selected with runtime flags. |
219 |
|
|
220 |
Note, that by default the seaice-package includes the orginial |
\paragraph{Compatibility with ice-thermodynamics package \code{thsice}\label{sec:pkg:seaice:thsice}}~\\ |
221 |
|
% |
222 |
|
Note, that by default the \code{seaice}-package includes the orginial |
223 |
so-called zero-layer thermodynamics following \citet{hib80} with a |
so-called zero-layer thermodynamics following \citet{hib80} with a |
224 |
snow cover as in \citet{zha98a}. The zero-layer thermodynamic model |
snow cover as in \citet{zha98a}. The zero-layer thermodynamic model |
225 |
assumes that ice does not store heat and, therefore, tends to |
assumes that ice does not store heat and, therefore, tends to |
226 |
exaggerate the seasonal variability in ice thickness. This |
exaggerate the seasonal variability in ice thickness. This |
227 |
exaggeration can be significantly reduced by using |
exaggeration can be significantly reduced by using |
228 |
\citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic model |
\citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic |
229 |
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 |
230 |
thermodynamic model has been reformulated by \citet{win00}. The |
\citet{win00}. The reformulation improves model physics by |
231 |
reformulation improves model physics by representing the brine content |
representing the brine content of the upper ice with a variable heat |
232 |
of the upper ice with a variable heat capacity. It also improves |
capacity. It also improves model numerics and consumes less computer |
233 |
model numerics and consumes less computer time and memory. The Winton |
time and memory. |
234 |
sea-ice thermodynamics have been ported to the MIT GCM; they currently |
|
235 |
reside under pkg/thsice. The package pkg/thsice is fully compatible |
The Winton sea-ice thermodynamics have been ported to the MIT GCM; |
236 |
with pkg/seaice and with pkg/exf. When turned on together with |
they currently reside under \code{pkg/thsice}. The package |
237 |
pkg/seaice, the zero-layer thermodynamics are replaced by the Winton |
\code{thsice} is described in section~\ref{sec:pkg:thsice}; it is |
238 |
thermodynamics. |
fully compatible with the packages \code{seaice} and \code{exf}. When |
239 |
|
turned on together with \code{seaice}, the zero-layer thermodynamics |
240 |
|
are replaced by the Winton thermodynamics. In order to use the |
241 |
|
\code{seaice}-package with the thermodynamics of \code{thsice}, |
242 |
|
compile both packages and turn both package on in \code{data.pkg}; see |
243 |
|
an example in \code{global\_ocean.cs32x15/input.icedyn}. Note, that |
244 |
|
once \code{thsice} is turned on, the variables and diagnostics |
245 |
|
associated to the default thermodynamics are meaningless, and the |
246 |
|
diagnostics of \code{thsice} have to be used instead. |
247 |
|
|
248 |
|
\paragraph{Surface forcing\label{sec:pkg:seaice:surfaceforcing}}~\\ |
249 |
|
% |
250 |
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 |
251 |
air temperature and specific humidity, downward longwave and shortwave |
air temperature and specific humidity, downward longwave and shortwave |
252 |
radiations, precipitation, evaporation, and river and glacier runoff. |
radiations, precipitation, evaporation, and river and glacier runoff. |
257 |
global: in ice-free regions bulk formulae are used to estimate oceanic |
global: in ice-free regions bulk formulae are used to estimate oceanic |
258 |
forcing from the atmospheric fields. |
forcing from the atmospheric fields. |
259 |
|
|
260 |
\paragraph{Dynamics\label{sec:pkg:seaice:dynamics}} |
\paragraph{Dynamics\label{sec:pkg:seaice:dynamics}}~\\ |
261 |
|
% |
262 |
\newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}} |
\newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}} |
263 |
\newcommand{\vtau}{\vek{\mathbf{\tau}}} |
\newcommand{\vtau}{\vek{\mathbf{\tau}}} |
264 |
The momentum equation of the sea-ice model is |
The momentum equation of the sea-ice model is |
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\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 |
303 |
viscous-plastic (VP) constitutive law \citep{hib79, zhang97}: |
viscous-plastic (VP) constitutive law \citep{hib79, zhang97}: |
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 |
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 so-called truncated ellipse method the shear viscosity $\eta$ |
Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in |
353 |
is capped to suppress any tensile stress \citep{hibler97, geiger98}: |
\code{SEAICE\_OPTIONS.h} before compiling replaces the method for |
354 |
|
bounding $\zeta$ by a smooth (differentiable) expression: |
355 |
\begin{equation} |
\begin{equation} |
356 |
\label{eq:etatem} |
\label{eq:zetaregsmooth} |
357 |
\eta = \min\left(\frac{\zeta}{e^2}, |
\begin{split} |
358 |
\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} |
\zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min}) |
359 |
{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2 |
\,\zeta_{\max}}\right)\\ |
360 |
+4\dot{\epsilon}_{12}^2}}\right). |
&= \frac{P}{2\Delta^*} |
361 |
|
\tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right) |
362 |
|
\end{split} |
363 |
\end{equation} |
\end{equation} |
364 |
To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in |
where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions |
365 |
\code{SEAICE\_OPTIONS.h} and turn it on with |
by zero. |
|
\code{SEAICEuseTEM=.TRUE.} in \code{data.seaice}. |
|
366 |
|
|
367 |
In the current implementation, the VP-model is integrated with the |
\paragraph{LSR and JFNK solver \label{sec:pkg:seaice:LSRJFNK}}~\\ |
368 |
semi-implicit line successive over relaxation (LSOR)-solver of |
% |
369 |
\citet{zhang97}, which allows for long time steps that, in our case, |
% By default, the VP-model is integrated by a Pickwith the |
370 |
are limited by the explicit treatment of the Coriolis term. The |
% semi-implicit line successive over relaxation (LSOR)-solver of |
371 |
explicit treatment of the Coriolis term does not represent a severe |
% \citet{zhang97}, which allows for long time steps that, in our case, |
372 |
limitation because it restricts the time step to approximately the |
% are limited by the explicit treatment of the Coriolis term. The |
373 |
same length as in the ocean model where the Coriolis term is also |
% explicit treatment of the Coriolis term does not represent a severe |
374 |
treated explicitly. |
% 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}}~\\ |
541 |
|
% |
542 |
\citet{hun97}'s introduced an elastic contribution to the strain |
\citet{hun97}'s introduced an elastic contribution to the strain |
543 |
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 |
544 |
the resulting elastic-viscous-plastic (EVP) and VP models are |
the resulting elastic-viscous-plastic (EVP) and VP models are |
556 |
%used and compared the present sea-ice model study. |
%used and compared the present sea-ice model study. |
557 |
The EVP-model uses an explicit time stepping scheme with a short |
The EVP-model uses an explicit time stepping scheme with a short |
558 |
timestep. According to the recommendation of \citet{hun97}, the |
timestep. According to the recommendation of \citet{hun97}, the |
559 |
EVP-model is stepped forward in time 120 times within the physical |
EVP-model should be stepped forward in time 120 times |
560 |
ocean model time step (although this parameter is under debate), to |
($\code{SEAICE\_deltaTevp} = \code{SEAICIE\_deltaTdyn}/120$) within |
561 |
allow for elastic waves to disappear. Because the scheme does not |
the physical ocean model time step (although this parameter is under |
562 |
require a matrix inversion it is fast in spite of the small internal |
debate), to allow for elastic waves to disappear. Because the scheme |
563 |
timestep and simple to implement on parallel computers |
does not require a matrix inversion it is fast in spite of the small |
564 |
|
internal timestep and simple to implement on parallel computers |
565 |
\citep{hun97}. For completeness, we repeat the equations for the |
\citep{hun97}. For completeness, we repeat the equations for the |
566 |
components of the stress tensor $\sigma_{1} = |
components of the stress tensor $\sigma_{1} = |
567 |
\sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and |
\sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and |
582 |
\frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T} |
\frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T} |
583 |
&= \frac{P}{4T\Delta} D_S |
&= \frac{P}{4T\Delta} D_S |
584 |
\end{align} |
\end{align} |
585 |
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 |
586 |
$T$ for elastic waves \[E=\frac{\zeta}{T}.\] |
timescale $T$ for elastic waves \[E=\frac{\zeta}{T}.\] |
587 |
$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 |
588 |
the external (long) timestep $\Delta{t}$. \citet{hun97} recommend |
(long) timestep $\Delta{t}$. $E_{0} = \frac{1}{3}$ is the default |
589 |
$E_{0} = \frac{1}{3}$ (which is the default value in the code). |
value in the code and close to what \citet{hun97} and |
590 |
|
\citet{hun01} recommend. |
591 |
|
|
592 |
To use the EVP solver, make sure that both \code{SEAICE\_CGRID} and |
To use the EVP solver, make sure that both \code{SEAICE\_CGRID} and |
593 |
\code{SEAICE\_ALLOW\_EVP} are defined in \code{SEAICE\_OPTIONS.h} |
\code{SEAICE\_ALLOW\_EVP} are defined in \code{SEAICE\_OPTIONS.h} |
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 variant of Elastic-Viscous-Plastic Dynamics: EVP*\label{sec:pkg:seaice:EVPstar}}~\\ |
606 |
|
% |
607 |
|
The genuine EVP schemes appears to give noisy solutions \citep{hun01, |
608 |
|
lemieux12, bouillon13}. This has lead to a modified EVP or EVP* |
609 |
|
\citep{lemieux12, bouillon13, kimmritz15}; here, refer to these |
610 |
|
variants by EVP*. The main idea is to modify the ``natural'' |
611 |
|
time-discretization of the momentum equations: |
612 |
|
\begin{equation} |
613 |
|
\label{eq:evpstar} |
614 |
|
m\frac{D\vec{u}}{Dt} \approx m\frac{u^{p+1}-u^{n}}{\Delta{t}} |
615 |
|
+ \beta^{*}\frac{u^{p+1}-u^{p}}{\Delta{t}_{\mathrm{EVP}}} |
616 |
|
\end{equation} |
617 |
|
where $n$ is the previous time step index, and $p$ is the previous |
618 |
|
sub-cycling index. The extra ``intertial'' term |
619 |
|
$m\,(u^{p+1}-u^{n})/\Delta{t})$ allows the definition of a residual |
620 |
|
$|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to |
621 |
|
$0$. In this way EVP can be re-interpreted as a pure iterative solver |
622 |
|
where the sub-cycling has no association with time-relation (through |
623 |
|
$\Delta{t}_{\mathrm{EVP}}$) \citep{bouillon13, kimmritz15}. Using the |
624 |
|
terminology of \citet{kimmritz15}, the evolution equations of stress |
625 |
|
$\sigma_{ij}$ and momentum $\vec{u}$ can be written as: |
626 |
|
\begin{align} |
627 |
|
\label{eq:evpstarsigma} |
628 |
|
\sigma_{ij}^{p+1}&=\sigma_{ij}^p+\frac{1}{\alpha} |
629 |
|
\Big(\sigma_{ij}(\vec{u}^p)-\sigma_{ij}^p\Big), |
630 |
|
\phantom{\int}\\ |
631 |
|
\label{eq:evpstarmom} |
632 |
|
\vec{u}^{p+1}&=\vec{u}^p+\frac{1}{\beta} |
633 |
|
\Big(\frac{\Delta t}{m}\nabla \cdot{\bf \sigma}^{p+1}+ |
634 |
|
\frac{\Delta t}{m}\vec{R}^{p}+\vec{u}_n-\vec{u}^p\Big). |
635 |
|
\end{align} |
636 |
|
$\vec{R}$ contains all terms in the momentum equations except for the |
637 |
|
rheology terms and the time derivative; $\alpha$ and $\beta$ are free |
638 |
|
parameters (\code{SEAICE\_evpAlpha}, \code{SEAICE\_evpBeta}) that |
639 |
|
replace the time stepping parameters \code{SEAICE\_deltaTevp} |
640 |
|
($\Delta{T}_{\mathrm{EVP}}$), \code{SEAICE\_elasticParm} ($E_{0}$), or |
641 |
|
\code{SEAICE\_evpTauRelax} ($T$). $\alpha$ and $\beta$ determine the |
642 |
|
speed of convergence and the stability. Usually, it makes sense to use |
643 |
|
$\alpha = \beta$, and \code{SEAICEnEVPstarSteps} $\gg |
644 |
|
(\alpha,\,\beta)$ \citep{kimmritz15}. Currently, there is no |
645 |
|
termination criterion and the number of EVP* iterations is fixed to |
646 |
|
\code{SEAICEnEVPstarSteps}. |
647 |
|
|
648 |
|
In order to use EVP* in the MITgcm, set \code{SEAICEuseEVPstar = |
649 |
|
.TRUE.,} in \code{data.seaice}. If \code{SEAICEuseEVPrev =.TRUE.,} |
650 |
|
the actual form of equations (\ref{eq:evpstarsigma}) and |
651 |
|
(\ref{eq:evpstarmom}) is used with fewer implicit terms and the factor |
652 |
|
of $e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2}) |
653 |
|
and (\ref{eq:evpstresstensor12}). Although this modifies the original |
654 |
|
EVP-equations, it turns out to improve convergence \citep{bouillon13}. |
655 |
|
|
656 |
|
Note, that for historical reasons, \code{SEAICE\_deltaTevp} needs to |
657 |
|
be set to some (any!) value in order to use also EVP*; this behavoir |
658 |
|
many change in the future. Also note, that |
659 |
|
probably because of the C-grid staggering of velocities and stresses, |
660 |
|
EVP* does not converge as successfully as in \citet{kimmritz15}. |
661 |
|
|
662 |
|
\paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\ |
663 |
|
% |
664 |
|
In the so-called truncated ellipse method the shear viscosity $\eta$ |
665 |
|
is capped to suppress any tensile stress \citep{hibler97, geiger98}: |
666 |
|
\begin{equation} |
667 |
|
\label{eq:etatem} |
668 |
|
\eta = \min\left(\frac{\zeta}{e^2}, |
669 |
|
\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} |
670 |
|
{\sqrt{\max(\Delta_{\min}^{2},(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2 |
671 |
|
+4\dot{\epsilon}_{12}^2})}\right). |
672 |
|
\end{equation} |
673 |
|
To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in |
674 |
|
\code{SEAICE\_OPTIONS.h} and turn it on with |
675 |
|
\code{SEAICEuseTEM} in \code{data.seaice}. |
676 |
|
|
677 |
|
\paragraph{Ice-Ocean stress \label{sec:pkg:seaice:iceoceanstress}}~\\ |
678 |
|
% |
679 |
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 |
680 |
the stress $\vtau_{ocean}$ in Eq.~\ref{eq:momseaice}. This stess is |
the stress $\vtau_{ocean}$ in Eq.~\ref{eq:momseaice}. This stess is |
681 |
applied directly to the surface layer of the ocean model. An |
applied directly to the surface layer of the ocean model. An |
705 |
% $P$ at vorticity points. |
% $P$ at vorticity points. |
706 |
|
|
707 |
\paragraph{Finite-volume discretization of the stress tensor |
\paragraph{Finite-volume discretization of the stress tensor |
708 |
divergence\label{sec:pkg:seaice:discretization}} |
divergence\label{sec:pkg:seaice:discretization}}~\\ |
709 |
|
% |
710 |
On an Arakawa C~grid, ice thickness and concentration and thus ice |
On an Arakawa C~grid, ice thickness and concentration and thus ice |
711 |
strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are |
strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are |
712 |
naturally defined a C-points in the center of the grid |
naturally defined a C-points in the center of the grid |
768 |
[(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2 |
[(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2 |
769 |
]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence |
]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence |
770 |
$(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is |
$(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is |
771 |
discretized in finite volumes. This conveniently avoids dealing with |
discretized in finite volumes \citep[see |
772 |
|
also][]{losch10:_mitsim}. This conveniently avoids dealing with |
773 |
further metric terms, as these are ``hidden'' in the differential cell |
further metric terms, as these are ``hidden'' in the differential cell |
774 |
widths. For the $u$-equation ($\alpha=1$) we have: |
widths. For the $u$-equation ($\alpha=1$) we have: |
775 |
\begin{align} |
\begin{align} |
868 |
analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set |
analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set |
869 |
$\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. |
870 |
|
|
871 |
\paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}} |
\paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\ |
872 |
|
% |
873 |
|
\noindent\textbf{NOTE: THIS SECTION IS TERRIBLY OUT OF DATE}\\ |
874 |
In its original formulation the sea ice model \citep{menemenlis05} |
In its original formulation the sea ice model \citep{menemenlis05} |
875 |
uses simple thermodynamics following the appendix of |
uses simple thermodynamics following the appendix of |
876 |
\citet{sem76}. This formulation does not allow storage of heat, |
\citet{sem76}. This formulation does not allow storage of heat, |
886 |
The conductive heat flux depends strongly on the ice thickness $h$. |
The conductive heat flux depends strongly on the ice thickness $h$. |
887 |
However, the ice thickness in the model represents a mean over a |
However, the ice thickness in the model represents a mean over a |
888 |
potentially very heterogeneous thickness distribution. In order to |
potentially very heterogeneous thickness distribution. In order to |
889 |
parameterize a sub-grid scale distribution for heat flux |
parameterize a sub-grid scale distribution for heat flux computations, |
890 |
computations, the mean ice thickness $h$ is split into seven thickness |
the mean ice thickness $h$ is split into $N$ thickness categories |
891 |
categories $H_{n}$ that are equally distributed between $2h$ and a |
$H_{n}$ that are equally distributed between $2h$ and a minimum |
892 |
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$ |
893 |
\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 |
894 |
thickness category is area-averaged to give the total heat flux |
is area-averaged to give the total heat flux \citep{hibler84}. To use |
895 |
\citep{hibler84}. To use this thickness category parameterization set |
this thickness category parameterization set \code{SEAICE\_multDim} to |
896 |
\code{\#define SEAICE\_MULTICATEGORY}; note that this requires |
the number of desired categories (7 is a good guess, for anything |
897 |
different restart files and switching this flag on in the middle of an |
larger than 7 modify \code{SEAICE\_SIZE.h}) in |
898 |
integration is not possible. |
\code{data.seaice}; note that this requires different restart files |
899 |
|
and switching this flag on in the middle of an integration is not |
900 |
|
advised. In order to include the same distribution for snow, set |
901 |
|
\code{SEAICE\_useMultDimSnow = .TRUE.}; only then, the |
902 |
|
parameterization of always having a fraction of thin ice is efficient |
903 |
|
and generally thicker ice is produced \citep{castro-morales14}. |
904 |
|
|
905 |
|
|
906 |
The atmospheric heat flux is balanced by an oceanic heat flux from |
The atmospheric heat flux is balanced by an oceanic heat flux from |
907 |
below. The oceanic flux is proportional to |
below. The oceanic flux is proportional to |
930 |
\code{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter |
\code{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter |
931 |
\code{SEAICEuseFlooding=.true.}. |
\code{SEAICEuseFlooding=.true.}. |
932 |
|
|
933 |
|
\paragraph{Advection of thermodynamic variables\label{sec:pkg:seaice:advection}}~\\ |
934 |
|
% |
935 |
Effective ice thickness (ice volume per unit area, |
Effective ice thickness (ice volume per unit area, |
936 |
$c\cdot{h}$), concentration $c$ and effective snow thickness |
$c\cdot{h}$), concentration $c$ and effective snow thickness |
937 |
($c\cdot{h}_{s}$) are advected by ice velocities: |
($c\cdot{h}_{s}$) are advected by ice velocities: |
951 |
(negative thickness or concentration). These schemes conserve volume |
(negative thickness or concentration). These schemes conserve volume |
952 |
and horizontal area and are unconditionally stable, so that we can set |
and horizontal area and are unconditionally stable, so that we can set |
953 |
$D_{X}=0$. Run-timeflags: \code{SEAICEadvScheme} (default=2, is the |
$D_{X}=0$. Run-timeflags: \code{SEAICEadvScheme} (default=2, is the |
954 |
historic 2nd-order, centered difference scheme), \code{DIFF1} |
historic 2nd-order, centered difference scheme), \code{DIFF1} = |
955 |
|
$D_{X}/\Delta{x}$ |
956 |
(default=0.004). |
(default=0.004). |
957 |
|
|
958 |
There is considerable doubt about the reliability of a ``zero-layer'' |
The MITgcm sea ice model provides the option to use |
959 |
thermodynamic model --- \citet{semtner84} found significant errors in |
the thermodynamics model of \citet{win00}, which in turn is based on |
960 |
phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in |
the 3-layer model of \citet{sem76} and which treats brine content by |
961 |
such models --- so that today many sea ice models employ more complex |
means of enthalpy conservation; the corresponding package |
962 |
thermodynamics. The MITgcm sea ice model provides the option to use |
\code{thsice} is described in section~\ref{sec:pkg:thsice}. This |
963 |
the thermodynamics model of \citet{win00}, which in turn is based |
scheme requires additional state variables, namely the enthalpy of the |
964 |
on the 3-layer model of \citet{sem76} and which treats brine |
two ice layers (instead of effective ice salinity), to be advected by |
965 |
content by means of enthalpy conservation. This scheme requires |
ice velocities. |
|
additional state variables, namely the enthalpy of the two ice layers |
|
|
(instead of effective ice salinity), to be advected by ice velocities. |
|
966 |
% |
% |
967 |
The internal sea ice temperature is inferred from ice enthalpy. To |
The internal sea ice temperature is inferred from ice enthalpy. To |
968 |
avoid unphysical (negative) values for ice thickness and |
avoid unphysical (negative) values for ice thickness and |
969 |
concentration, a positive 2nd-order advection scheme with a SuperBee |
concentration, a positive 2nd-order advection scheme with a SuperBee |
970 |
flux limiter \citep{roe:85} is used in this study to advect all |
flux limiter \citep{roe:85} should be used to advect all |
971 |
sea-ice-related quantities of the \citet{win00} thermodynamic |
sea-ice-related quantities of the \citet{win00} thermodynamic model |
972 |
model. Because of the non-linearity of the advection scheme, care |
(runtime flag \code{thSIceAdvScheme=77} and |
973 |
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 |
974 |
velocity to advect enthalpy, the total energy (i.e., the volume |
non-linearity of the advection scheme, care must be taken in advecting |
975 |
integral of enthalpy) is not conserved. Alternatively, one can advect |
these quantities: when simply using ice velocity to advect enthalpy, |
976 |
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 |
977 |
false enthalpy extrema can occur, which then leads to unrealistic ice |
conserved. Alternatively, one can advect the energy content (i.e., |
978 |
temperature. In the currently implemented solution, the sea-ice mass |
product of ice-volume and enthalpy) but then false enthalpy extrema |
979 |
flux is used to advect the enthalpy in order to ensure conservation of |
can occur, which then leads to unrealistic ice temperature. In the |
980 |
enthalpy and to prevent false enthalpy extrema. |
currently implemented solution, the sea-ice mass flux is used to |
981 |
|
advect the enthalpy in order to ensure conservation of enthalpy and to |
982 |
|
prevent false enthalpy extrema. % |
983 |
|
|
984 |
%---------------------------------------------------------------------- |
%---------------------------------------------------------------------- |
985 |
|
|
1045 |
Available output fields are summarized in |
Available output fields are summarized in |
1046 |
Table \ref{tab:pkg:seaice:diagnostics}. |
Table \ref{tab:pkg:seaice:diagnostics}. |
1047 |
|
|
1048 |
\begin{table}[!ht] |
\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} |
|
|
|
|
1049 |
|
|
1050 |
%\subsubsection{Package Reference} |
%\subsubsection{Package Reference} |
1051 |
|
|
1053 |
\label{sec:pkg:seaice:experiments} |
\label{sec:pkg:seaice:experiments} |
1054 |
|
|
1055 |
\begin{itemize} |
\begin{itemize} |
1056 |
\item{Labrador Sea experiment in lab\_sea verification directory. } |
\item{Labrador Sea experiment in \code{lab\_sea} verification directory. } |
1057 |
|
\item \code{seaice\_obcs}, based on \code{lab\_sea} |
1058 |
|
\item \code{offline\_exf\_seaice/input.seaicetd}, based on \code{lab\_sea} |
1059 |
|
\item \code{global\_ocean.cs32x15/input.icedyn} and |
1060 |
|
\code{global\_ocean.cs32x15/input.seaice}, global |
1061 |
|
cubed-sphere-experiment with combinations of \code{seaice} and |
1062 |
|
\code{thsice} |
1063 |
\end{itemize} |
\end{itemize} |
1064 |
|
|
1065 |
|
|