| 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 | 
| 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 | 
 \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}[h!] | 
 \begin{table}[!ht] | 
| 70 | 
 \centering | 
 \centering | 
| 71 | 
   \label{tab:pkg:seaice:cpp} | 
   \label{tab:pkg:seaice:cpp} | 
| 72 | 
   {\footnotesize | 
   {\footnotesize | 
| 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 variants of Elastic-Viscous-Plastic Dynamics: | 
| 606 | 
  | 
   EVP* , mEVP, and aEVP \label{sec:pkg:seaice:EVPstar}}~\\  | 
| 607 | 
  | 
 % | 
| 608 | 
  | 
 The genuine EVP schemes appears to give noisy solutions \citep{hun01, | 
| 609 | 
  | 
   lemieux12, bouillon13}. This has lead to a modified EVP or EVP* | 
| 610 | 
  | 
 \citep{lemieux12, bouillon13, kimmritz15}; here, we refer to these | 
| 611 | 
  | 
 variants by modified EVP (mEVP) and adaptive EVP (aEVP) | 
| 612 | 
  | 
 \citep{kimmritz16}. The main idea is to modify the ``natural''  | 
| 613 | 
  | 
 time-discretization of the momentum equations: | 
| 614 | 
  | 
 \begin{equation} | 
| 615 | 
  | 
   \label{eq:evpstar} | 
| 616 | 
  | 
   m\frac{D\vec{u}}{Dt} \approx m\frac{u^{p+1}-u^{n}}{\Delta{t}} | 
| 617 | 
  | 
   + \beta^{*}\frac{u^{p+1}-u^{p}}{\Delta{t}_{\mathrm{EVP}}} | 
| 618 | 
  | 
 \end{equation} | 
| 619 | 
  | 
 where $n$ is the previous time step index, and $p$ is the previous | 
| 620 | 
  | 
 sub-cycling index. The extra ``intertial'' term | 
| 621 | 
  | 
 $m\,(u^{p+1}-u^{n})/\Delta{t})$ allows the definition of a residual | 
| 622 | 
  | 
 $|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to | 
| 623 | 
  | 
 $0$. In this way EVP can be re-interpreted as a pure iterative solver | 
| 624 | 
  | 
 where the sub-cycling has no association with time-relation (through | 
| 625 | 
  | 
 $\Delta{t}_{\mathrm{EVP}}$) \citep{bouillon13, kimmritz15}. Using the | 
| 626 | 
  | 
 terminology of \citet{kimmritz15}, the evolution equations of stress | 
| 627 | 
  | 
 $\sigma_{ij}$ and momentum $\vec{u}$ can be written as: | 
| 628 | 
  | 
 \begin{align} | 
| 629 | 
  | 
   \label{eq:evpstarsigma} | 
| 630 | 
  | 
   \sigma_{ij}^{p+1}&=\sigma_{ij}^p+\frac{1}{\alpha} | 
| 631 | 
  | 
   \Big(\sigma_{ij}(\vec{u}^p)-\sigma_{ij}^p\Big), | 
| 632 | 
  | 
   \phantom{\int}\\ | 
| 633 | 
  | 
   \label{eq:evpstarmom} | 
| 634 | 
  | 
   \vec{u}^{p+1}&=\vec{u}^p+\frac{1}{\beta} | 
| 635 | 
  | 
   \Big(\frac{\Delta t}{m}\nabla \cdot{\bf \sigma}^{p+1}+ | 
| 636 | 
  | 
   \frac{\Delta t}{m}\vec{R}^{p}+\vec{u}_n-\vec{u}^p\Big). | 
| 637 | 
  | 
 \end{align} | 
| 638 | 
  | 
 $\vec{R}$ contains all terms in the momentum equations except for the | 
| 639 | 
  | 
 rheology terms and the time derivative; $\alpha$ and $\beta$ are free | 
| 640 | 
  | 
 parameters (\code{SEAICE\_evpAlpha}, \code{SEAICE\_evpBeta}) that | 
| 641 | 
  | 
 replace the time stepping parameters \code{SEAICE\_deltaTevp} | 
| 642 | 
  | 
 ($\Delta{T}_{\mathrm{EVP}}$), \code{SEAICE\_elasticParm} ($E_{0}$), or | 
| 643 | 
  | 
 \code{SEAICE\_evpTauRelax} ($T$). $\alpha$ and $\beta$ determine the | 
| 644 | 
  | 
 speed of convergence and the stability. Usually, it makes sense to use | 
| 645 | 
  | 
 $\alpha = \beta$, and \code{SEAICEnEVPstarSteps} $\gg | 
| 646 | 
  | 
 (\alpha,\,\beta)$ \citep{kimmritz15}. Currently, there is no | 
| 647 | 
  | 
 termination criterion and the number of mEVP iterations is fixed to | 
| 648 | 
  | 
 \code{SEAICEnEVPstarSteps}. | 
| 649 | 
  | 
  | 
| 650 | 
  | 
 In order to use mEVP in the MITgcm, set \code{SEAICEuseEVPstar = | 
| 651 | 
  | 
   .TRUE.,} in \code{data.seaice}. If \code{SEAICEuseEVPrev =.TRUE.,} | 
| 652 | 
  | 
 the actual form of equations (\ref{eq:evpstarsigma}) and | 
| 653 | 
  | 
 (\ref{eq:evpstarmom}) is used with fewer implicit terms and the factor | 
| 654 | 
  | 
 of $e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2}) | 
| 655 | 
  | 
 and (\ref{eq:evpstresstensor12}). Although this modifies the original | 
| 656 | 
  | 
 EVP-equations, it turns out to improve convergence \citep{bouillon13}. | 
| 657 | 
  | 
  | 
| 658 | 
  | 
 Another variant is the aEVP scheme \citep{kimmritz16}, where the value | 
| 659 | 
  | 
 of $\alpha$ is set dynamically based on the stability criterion | 
| 660 | 
  | 
 \begin{equation} | 
| 661 | 
  | 
   \label{eq:aevpalpha} | 
| 662 | 
  | 
   \alpha = \beta = \max\left( \tilde{c}\pi\sqrt{c \frac{\zeta}{A_{c}} | 
| 663 | 
  | 
     \frac{\Delta{t}}{\max(m,10^{-4}\text{\,kg})}},\alpha_{\min} \right) | 
| 664 | 
  | 
 \end{equation} | 
| 665 | 
  | 
 with the grid cell area $A_c$ and the ice and snow mass $m$.  This | 
| 666 | 
  | 
 choice sacrifices speed of convergence for stability with the result | 
| 667 | 
  | 
 that aEVP converges quickly to VP where $\alpha$ can be small and more | 
| 668 | 
  | 
 slowly in areas where the equations are stiff. In practice, aEVP leads | 
| 669 | 
  | 
 to an overall better convergence than mEVP \citep{kimmritz16}. | 
| 670 | 
  | 
 %  | 
| 671 | 
  | 
 To use aEVP in the MITgcm set \code{SEAICEaEVPcoeff} $= \tilde{c}$; | 
| 672 | 
  | 
 this also sets the default values of \code{SEAICEaEVPcStar} ($c=4$) | 
| 673 | 
  | 
 and \code{SEAICEaEVPalphaMin} ($\alpha_{\min}=5$). Good convergence | 
| 674 | 
  | 
 has been obtained with setting these values \citep{kimmritz16}: | 
| 675 | 
  | 
 \code{SEAICEaEVPcoeff = 0.5, SEAICEnEVPstarSteps = 500, | 
| 676 | 
  | 
   SEAICEuseEVPstar = .TRUE., SEAICEuseEVPrev = .TRUE.} | 
| 677 | 
  | 
  | 
| 678 | 
  | 
 Note, that probably because of the C-grid staggering of velocities and | 
| 679 | 
  | 
 stresses, mEVP may not converge as successfully as in | 
| 680 | 
  | 
 \citet{kimmritz15}, and that convergence at very high resolution | 
| 681 | 
  | 
 (order 5\,km) has not been studied yet. | 
| 682 | 
  | 
  | 
| 683 | 
  | 
 \paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\ | 
| 684 | 
  | 
 % | 
| 685 | 
  | 
 In the so-called truncated ellipse method the shear viscosity $\eta$ | 
| 686 | 
  | 
 is capped to suppress any tensile stress \citep{hibler97, geiger98}: | 
| 687 | 
  | 
 \begin{equation} | 
| 688 | 
  | 
   \label{eq:etatem} | 
| 689 | 
  | 
   \eta = \min\left(\frac{\zeta}{e^2}, | 
| 690 | 
  | 
   \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})} | 
| 691 | 
  | 
   {\sqrt{\max(\Delta_{\min}^{2},(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2 | 
| 692 | 
  | 
       +4\dot{\epsilon}_{12}^2})}\right). | 
| 693 | 
  | 
 \end{equation} | 
| 694 | 
  | 
 To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in | 
| 695 | 
  | 
 \code{SEAICE\_OPTIONS.h} and turn it on with | 
| 696 | 
  | 
 \code{SEAICEuseTEM} in \code{data.seaice}.  | 
| 697 | 
  | 
  | 
| 698 | 
  | 
 \paragraph{Ice-Ocean stress \label{sec:pkg:seaice:iceoceanstress}}~\\ | 
| 699 | 
  | 
 % | 
| 700 | 
 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 | 
| 701 | 
 the stress $\vtau_{ocean}$ in Eq.~\ref{eq:momseaice}. This stess is | 
 the stress $\vtau_{ocean}$ in Eq.~\ref{eq:momseaice}. This stess is | 
| 702 | 
 applied directly to the surface layer of the ocean model. An | 
 applied directly to the surface layer of the ocean model. An | 
| 726 | 
 % $P$ at vorticity points. | 
 % $P$ at vorticity points. | 
| 727 | 
  | 
  | 
| 728 | 
 \paragraph{Finite-volume discretization of the stress tensor | 
 \paragraph{Finite-volume discretization of the stress tensor | 
| 729 | 
   divergence\label{sec:pkg:seaice:discretization}} | 
   divergence\label{sec:pkg:seaice:discretization}}~\\ | 
| 730 | 
  | 
 % | 
| 731 | 
 On an Arakawa C~grid, ice thickness and concentration and thus ice | 
 On an Arakawa C~grid, ice thickness and concentration and thus ice | 
| 732 | 
 strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are | 
 strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are | 
| 733 | 
 naturally defined a C-points in the center of the grid | 
 naturally defined a C-points in the center of the grid | 
| 789 | 
 [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2 | 
 [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2 | 
| 790 | 
 ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence | 
 ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence | 
| 791 | 
 $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is | 
 $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is | 
| 792 | 
 discretized in finite volumes. This conveniently avoids dealing with | 
 discretized in finite volumes \citep[see | 
| 793 | 
  | 
 also][]{losch10:_mitsim}. This conveniently avoids dealing with | 
| 794 | 
 further metric terms, as these are ``hidden'' in the differential cell | 
 further metric terms, as these are ``hidden'' in the differential cell | 
| 795 | 
 widths. For the $u$-equation ($\alpha=1$) we have: | 
 widths. For the $u$-equation ($\alpha=1$) we have: | 
| 796 | 
 \begin{align} | 
 \begin{align} | 
| 813 | 
   \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{} | 
   \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{} | 
| 814 | 
   + (\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 | 
| 815 | 
   \biggr\} | 
   \biggr\} | 
| 816 | 
   \intertext{with} | 
 \end{align} | 
| 817 | 
  | 
 with | 
| 818 | 
  | 
 \begin{align} | 
| 819 | 
   (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+} | 
   (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+} | 
| 820 | 
   \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j} | 
   \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j} | 
| 821 | 
   \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 | 
| 857 | 
   \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{} | 
   \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{} | 
| 858 | 
   + (\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 | 
| 859 | 
   \biggr\}  | 
   \biggr\}  | 
| 860 | 
   \intertext{with} | 
 \end{align} | 
| 861 | 
  | 
 with | 
| 862 | 
  | 
 \begin{align} | 
| 863 | 
   (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+} | 
   (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+} | 
| 864 | 
   \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j} | 
   \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j} | 
| 865 | 
   \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}}  | 
   \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}}  | 
| 889 | 
 analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set | 
 analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set | 
| 890 | 
 $\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. | 
| 891 | 
  | 
  | 
| 892 | 
 \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}} | 
 \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\ | 
| 893 | 
  | 
 % | 
| 894 | 
  | 
 \noindent\textbf{NOTE: THIS SECTION IS TERRIBLY OUT OF DATE}\\ | 
| 895 | 
 In its original formulation the sea ice model \citep{menemenlis05} | 
 In its original formulation the sea ice model \citep{menemenlis05} | 
| 896 | 
 uses simple thermodynamics following the appendix of | 
 uses simple thermodynamics following the appendix of | 
| 897 | 
 \citet{sem76}. This formulation does not allow storage of heat, | 
 \citet{sem76}. This formulation does not allow storage of heat, | 
| 907 | 
 The conductive heat flux depends strongly on the ice thickness $h$. | 
 The conductive heat flux depends strongly on the ice thickness $h$. | 
| 908 | 
 However, the ice thickness in the model represents a mean over a | 
 However, the ice thickness in the model represents a mean over a | 
| 909 | 
 potentially very heterogeneous thickness distribution.  In order to | 
 potentially very heterogeneous thickness distribution.  In order to | 
| 910 | 
 parameterize a sub-grid scale distribution for heat flux | 
 parameterize a sub-grid scale distribution for heat flux computations, | 
| 911 | 
 computations, the mean ice thickness $h$ is split into seven thickness | 
 the mean ice thickness $h$ is split into $N$ thickness categories | 
| 912 | 
 categories $H_{n}$ that are equally distributed between $2h$ and a | 
 $H_{n}$ that are equally distributed between $2h$ and a minimum | 
| 913 | 
 minimum imposed ice thickness of $5\text{\,cm}$ by $H_n= | 
 imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$ | 
| 914 | 
 \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each | 
 for $n\in[1,N]$. The heat fluxes computed for each thickness category | 
| 915 | 
 thickness category is area-averaged to give the total heat flux | 
 is area-averaged to give the total heat flux \citep{hibler84}. To use | 
| 916 | 
 \citep{hibler84}. To use this thickness category parameterization set | 
 this thickness category parameterization set \code{SEAICE\_multDim} to | 
| 917 | 
 \code{\#define SEAICE\_MULTICATEGORY}; note that this requires | 
 the number of desired categories (7 is a good guess, for anything | 
| 918 | 
 different restart files and switching this flag on in the middle of an | 
 larger than 7 modify \code{SEAICE\_SIZE.h}) in | 
| 919 | 
 integration is not possible. | 
 \code{data.seaice}; note that this requires different restart files | 
| 920 | 
  | 
 and switching this flag on in the middle of an integration is not | 
| 921 | 
  | 
 advised. In order to include the same distribution for snow, set | 
| 922 | 
  | 
 \code{SEAICE\_useMultDimSnow = .TRUE.}; only then, the | 
| 923 | 
  | 
 parameterization of always having a fraction of thin ice is efficient | 
| 924 | 
  | 
 and generally thicker ice is produced \citep{castro-morales14}. | 
| 925 | 
  | 
  | 
| 926 | 
  | 
  | 
| 927 | 
 The atmospheric heat flux is balanced by an oceanic heat flux from | 
 The atmospheric heat flux is balanced by an oceanic heat flux from | 
| 928 | 
 below.  The oceanic flux is proportional to | 
 below.  The oceanic flux is proportional to | 
| 951 | 
 \code{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter | 
 \code{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter | 
| 952 | 
 \code{SEAICEuseFlooding=.true.}. | 
 \code{SEAICEuseFlooding=.true.}. | 
| 953 | 
  | 
  | 
| 954 | 
  | 
 \paragraph{Advection of thermodynamic variables\label{sec:pkg:seaice:advection}}~\\ | 
| 955 | 
  | 
 % | 
| 956 | 
 Effective ice thickness (ice volume per unit area, | 
 Effective ice thickness (ice volume per unit area, | 
| 957 | 
 $c\cdot{h}$), concentration $c$ and effective snow thickness | 
 $c\cdot{h}$), concentration $c$ and effective snow thickness | 
| 958 | 
 ($c\cdot{h}_{s}$) are advected by ice velocities: | 
 ($c\cdot{h}_{s}$) are advected by ice velocities: | 
| 965 | 
 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})$. | 
| 966 | 
 % | 
 % | 
| 967 | 
 From the various advection scheme that are available in the MITgcm, we | 
 From the various advection scheme that are available in the MITgcm, we | 
| 968 | 
 choose flux-limited schemes \citep[multidimensional 2nd and 3rd-order | 
 recommend flux-limited schemes \citep[multidimensional 2nd and | 
| 969 | 
 advection scheme with flux limiter][]{roe:85, hundsdorfer94} to | 
 3rd-order advection scheme with flux limiter][]{roe:85, hundsdorfer94} | 
| 970 | 
 preserve sharp gradients and edges that are typical of sea ice | 
 to preserve sharp gradients and edges that are typical of sea ice | 
| 971 | 
 distributions and to rule out unphysical over- and undershoots | 
 distributions and to rule out unphysical over- and undershoots | 
| 972 | 
 (negative thickness or concentration). These scheme conserve volume | 
 (negative thickness or concentration). These schemes conserve volume | 
| 973 | 
 and horizontal area and are unconditionally stable, so that we can set | 
 and horizontal area and are unconditionally stable, so that we can set | 
| 974 | 
 $D_{X}=0$. Run-timeflags: \code{SEAICEadvScheme} (default=2), | 
 $D_{X}=0$. Run-timeflags: \code{SEAICEadvScheme} (default=2, is the | 
| 975 | 
 \code{DIFF1} (default=0.004). | 
 historic 2nd-order, centered difference scheme), \code{DIFF1} = | 
| 976 | 
  | 
 $D_{X}/\Delta{x}$ | 
| 977 | 
 There is considerable doubt about the reliability of a ``zero-layer'' | 
 (default=0.004). | 
| 978 | 
 thermodynamic model --- \citet{semtner84} found significant errors in | 
  | 
| 979 | 
 phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in | 
 The MITgcm sea ice model provides the option to use | 
| 980 | 
 such models --- so that today many sea ice models employ more complex | 
 the thermodynamics model of \citet{win00}, which in turn is based on | 
| 981 | 
 thermodynamics. The MITgcm sea ice model provides the option to use | 
 the 3-layer model of \citet{sem76} and which treats brine content by | 
| 982 | 
 the thermodynamics model of \citet{win00}, which in turn is based | 
 means of enthalpy conservation; the corresponding package | 
| 983 | 
 on the 3-layer model of \citet{sem76} and which treats brine | 
 \code{thsice} is described in section~\ref{sec:pkg:thsice}. This | 
| 984 | 
 content by means of enthalpy conservation. This scheme requires | 
 scheme requires additional state variables, namely the enthalpy of the | 
| 985 | 
 additional state variables, namely the enthalpy of the two ice layers | 
 two ice layers (instead of effective ice salinity), to be advected by | 
| 986 | 
 (instead of effective ice salinity), to be advected by ice velocities. | 
 ice velocities. | 
| 987 | 
 % | 
 % | 
| 988 | 
 The internal sea ice temperature is inferred from ice enthalpy.  To | 
 The internal sea ice temperature is inferred from ice enthalpy.  To | 
| 989 | 
 avoid unphysical (negative) values for ice thickness and | 
 avoid unphysical (negative) values for ice thickness and | 
| 990 | 
 concentration, a positive 2nd-order advection scheme with a SuperBee | 
 concentration, a positive 2nd-order advection scheme with a SuperBee | 
| 991 | 
 flux limiter \citep{roe:85} is used in this study to advect all | 
 flux limiter \citep{roe:85} should be used to advect all | 
| 992 | 
 sea-ice-related quantities of the \citet{win00} thermodynamic | 
 sea-ice-related quantities of the \citet{win00} thermodynamic model | 
| 993 | 
 model.  Because of the non-linearity of the advection scheme, care | 
 (runtime flag \code{thSIceAdvScheme=77} and | 
| 994 | 
 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 | 
| 995 | 
 velocity to advect enthalpy, the total energy (i.e., the volume | 
 non-linearity of the advection scheme, care must be taken in advecting | 
| 996 | 
 integral of enthalpy) is not conserved. Alternatively, one can advect | 
 these quantities: when simply using ice velocity to advect enthalpy, | 
| 997 | 
 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 | 
| 998 | 
 false enthalpy extrema can occur, which then leads to unrealistic ice | 
 conserved. Alternatively, one can advect the energy content (i.e., | 
| 999 | 
 temperature.  In the currently implemented solution, the sea-ice mass | 
 product of ice-volume and enthalpy) but then false enthalpy extrema | 
| 1000 | 
 flux is used to advect the enthalpy in order to ensure conservation of | 
 can occur, which then leads to unrealistic ice temperature.  In the | 
| 1001 | 
 enthalpy and to prevent false enthalpy extrema. | 
 currently implemented solution, the sea-ice mass flux is used to | 
| 1002 | 
  | 
 advect the enthalpy in order to ensure conservation of enthalpy and to | 
| 1003 | 
  | 
 prevent false enthalpy extrema. % | 
| 1004 | 
  | 
  | 
| 1005 | 
 %---------------------------------------------------------------------- | 
 %---------------------------------------------------------------------- | 
| 1006 | 
  | 
  | 
| 1066 | 
 Available output fields are summarized in  | 
 Available output fields are summarized in  | 
| 1067 | 
 Table \ref{tab:pkg:seaice:diagnostics}. | 
 Table \ref{tab:pkg:seaice:diagnostics}. | 
| 1068 | 
  | 
  | 
| 1069 | 
 \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} | 
  | 
 | 
  | 
  | 
| 1070 | 
  | 
  | 
| 1071 | 
 %\subsubsection{Package Reference} | 
 %\subsubsection{Package Reference} | 
| 1072 | 
  | 
  | 
| 1074 | 
 \label{sec:pkg:seaice:experiments} | 
 \label{sec:pkg:seaice:experiments} | 
| 1075 | 
  | 
  | 
| 1076 | 
 \begin{itemize} | 
 \begin{itemize} | 
| 1077 | 
 \item{Labrador Sea experiment in lab\_sea verification directory. } | 
 \item{Labrador Sea experiment in \code{lab\_sea} verification directory. } | 
| 1078 | 
  | 
 \item \code{seaice\_obcs}, based on \code{lab\_sea} | 
| 1079 | 
  | 
 \item \code{offline\_exf\_seaice/input.seaicetd}, based on \code{lab\_sea} | 
| 1080 | 
  | 
 \item \code{global\_ocean.cs32x15/input.icedyn} and | 
| 1081 | 
  | 
   \code{global\_ocean.cs32x15/input.seaice}, global | 
| 1082 | 
  | 
   cubed-sphere-experiment with combinations of \code{seaice} and | 
| 1083 | 
  | 
   \code{thsice} | 
| 1084 | 
 \end{itemize} | 
 \end{itemize} | 
| 1085 | 
  | 
  | 
| 1086 | 
  | 
  | 
| 1087 | 
 %%% Local Variables:  | 
 %%% Local Variables:  | 
| 1088 | 
 %%% mode: latex | 
 %%% mode: latex | 
| 1089 | 
 %%% TeX-master: "../manual" | 
 %%% TeX-master: "../../manual" | 
| 1090 | 
 %%% End:  | 
 %%% End:  |