| 1 | 
% $Header: $ | 
| 2 | 
% $Name: $ | 
| 3 | 
 | 
| 4 | 
\section{Example: Barotropic Ocean Gyre In Cartesian Coordinates} | 
| 5 | 
 | 
| 6 | 
\bodytext{bgcolor="#FFFFFFFF"} | 
| 7 | 
 | 
| 8 | 
%\begin{center}  | 
| 9 | 
%{\Large \bf Using MITgcm to Simulate a Barotropic Ocean Gyre In Cartesian | 
| 10 | 
%Coordinates} | 
| 11 | 
% | 
| 12 | 
%\vspace*{4mm} | 
| 13 | 
% | 
| 14 | 
%\vspace*{3mm} | 
| 15 | 
%{\large May 2001} | 
| 16 | 
%\end{center} | 
| 17 | 
 | 
| 18 | 
\subsection{Introduction} | 
| 19 | 
 | 
| 20 | 
This document is the first in a series of documents describing | 
| 21 | 
example MITgcm numerical experiments. The example experiments  | 
| 22 | 
include both straightforward examples of idealised geophysical  | 
| 23 | 
fluid simulations and more involved cases encompassing | 
| 24 | 
large scale modeling and | 
| 25 | 
automatic differentiation. Both hydrostatic and non-hydrostatic  | 
| 26 | 
experiements are presented, as well as experiments employing | 
| 27 | 
cartesian, spherical-polar and cube-sphere coordinate systems. | 
| 28 | 
These ``case study'' documents include information describing | 
| 29 | 
the experimental configuration and detailed information on how to | 
| 30 | 
configure the MITgcm code and input files for each experiment. | 
| 31 | 
 | 
| 32 | 
\subsection{Experiment Overview} | 
| 33 | 
 | 
| 34 | 
This example experiment demonstrates using the MITgcm to simulate | 
| 35 | 
a barotropic, wind-forced, ocean gyre circulation. The experiment  | 
| 36 | 
is a numerical rendition of the gyre circulation problem simliar | 
| 37 | 
to the problems described analytically by Stommel in 1966  | 
| 38 | 
\cite{Stommel66} and numerically in Holland et. al \cite{Holland75}. | 
| 39 | 
 | 
| 40 | 
In this experiment the model  | 
| 41 | 
is configured to represent a rectangular enclosed box of fluid, | 
| 42 | 
$1200 \times 1200 $~km in lateral extent. The fluid is $5$~km deep and is forced | 
| 43 | 
by a constant in time zonal wind stress, $\tau_x$, that varies sinusoidally | 
| 44 | 
in the ``north-south'' direction. Topologically the grid is cartesian and  | 
| 45 | 
the coriolis parameter $f$ is defined according to a mid-latitude beta-plane  | 
| 46 | 
equation | 
| 47 | 
 | 
| 48 | 
\begin{equation} | 
| 49 | 
\label{EQ:fcori} | 
| 50 | 
f(y) = f_{0}+\beta y | 
| 51 | 
\end{equation} | 
| 52 | 
  | 
| 53 | 
\noindent where $y$ is the distance along the ``north-south'' axis of the  | 
| 54 | 
simulated domain. For this experiment $f_{0}$ is set to $10^{-4}s^{-1}$ in  | 
| 55 | 
(\ref{EQ:fcori}) and $\beta = 10^{-11}s^{-1}m^{-1}$.  | 
| 56 | 
\\ | 
| 57 | 
\\ | 
| 58 | 
 The sinusoidal wind-stress variations are defined according to  | 
| 59 | 
 | 
| 60 | 
\begin{equation} | 
| 61 | 
\label{EQ:taux} | 
| 62 | 
\tau_x(y) = \tau_{0}\sin(\pi \frac{y}{L_y}) | 
| 63 | 
\end{equation} | 
| 64 | 
  | 
| 65 | 
\noindent where $L_{y}$ is the lateral domain extent ($1200$~km) and  | 
| 66 | 
$\tau_0$ is set to $0.1N m^{-2}$.  | 
| 67 | 
\\ | 
| 68 | 
\\ | 
| 69 | 
Figure \ref{FIG:simulation_config} | 
| 70 | 
summarises the configuration simulated. | 
| 71 | 
 | 
| 72 | 
\begin{figure} | 
| 73 | 
\centerline{ | 
| 74 | 
 \resizebox{7.5in}{5.5in}{ | 
| 75 | 
   \includegraphics*[0.2in,0.7in][10.5in,10.5in] | 
| 76 | 
    {part3/case_studies/barotropic_gyre/simulation_config.eps} } | 
| 77 | 
} | 
| 78 | 
\caption{Schematic of simulation domain and wind-stress forcing function  | 
| 79 | 
for barotropic gyre numerical experiment. The domain is enclosed bu solid | 
| 80 | 
walls at $x=$~0,1200km and at $y=$~0,1200km.} | 
| 81 | 
\label{FIG:simulation_config} | 
| 82 | 
\end{figure} | 
| 83 | 
 | 
| 84 | 
\subsection{Discrete Numerical Configuration} | 
| 85 | 
 | 
| 86 | 
 The model is configured in hydrostatic form.  The domain is discretised with  | 
| 87 | 
a uniform grid spacing in the horizontal set to | 
| 88 | 
 $\Delta x=\Delta y=20$~km, so  | 
| 89 | 
that there are sixty grid cells in the $x$ and $y$ directions. Vertically the  | 
| 90 | 
model is configured with a single layer with depth, $\Delta z$, of $5000$~m.  | 
| 91 | 
The implicit free surface form of the  | 
| 92 | 
pressure equation described in Marshall et. al \cite{Marshall97a} is  | 
| 93 | 
employed.  | 
| 94 | 
A horizontal laplacian operator $\nabla_{h}^2$ provides viscous | 
| 95 | 
dissipation. The wind-stress momentum input is added to the momentum equation | 
| 96 | 
for the ``zonal flow'', $u$. Other terms in the model | 
| 97 | 
are explicitly switched off for this experiement configuration (see section | 
| 98 | 
\ref{SEC:code_config} ), yielding an active set of equations solved in this  | 
| 99 | 
configuration as follows | 
| 100 | 
 | 
| 101 | 
\begin{eqnarray} | 
| 102 | 
\label{EQ:model_equations} | 
| 103 | 
\frac{Du}{Dt} - fv +  | 
| 104 | 
              g\frac{\partial \eta}{\partial x} -  | 
| 105 | 
              A_{h}\nabla_{h}^2u  | 
| 106 | 
& = & | 
| 107 | 
\frac{\tau_{x}}{\rho_{0}\Delta z} | 
| 108 | 
\\ | 
| 109 | 
\frac{Dv}{Dt} + fu + g\frac{\partial \eta}{\partial y} - | 
| 110 | 
              A_{h}\nabla_{h}^2v  | 
| 111 | 
& = & | 
| 112 | 
0 | 
| 113 | 
\\ | 
| 114 | 
\frac{\partial \eta}{\partial t} + \nabla_{h}\cdot \vec{u} | 
| 115 | 
&=& | 
| 116 | 
0 | 
| 117 | 
\end{eqnarray} | 
| 118 | 
 | 
| 119 | 
\noindent where $u$ and $v$ and the $x$ and $y$ components of the | 
| 120 | 
flow vector $\vec{u}$.  | 
| 121 | 
\\ | 
| 122 | 
 | 
| 123 | 
\subsubsection{Numerical Stability Criteria} | 
| 124 | 
 | 
| 125 | 
The laplacian dissipation coefficient, $A_{h}$, is set to $400 m s^{-1}$. | 
| 126 | 
This value is chosen to yield a Munk layer width \cite{Adcroft_thesis}, | 
| 127 | 
 | 
| 128 | 
\begin{eqnarray} | 
| 129 | 
\label{EQ:munk_layer} | 
| 130 | 
M_{w} = \pi ( \frac { A_{h} }{ \beta } )^{\frac{1}{3}} | 
| 131 | 
\end{eqnarray} | 
| 132 | 
 | 
| 133 | 
\noindent  of $\approx 100$km. This is greater than the model | 
| 134 | 
resolution $\Delta x$, ensuring that the frictional boundary | 
| 135 | 
layer is well resolved. | 
| 136 | 
\\ | 
| 137 | 
 | 
| 138 | 
\noindent The model is stepped forward with a  | 
| 139 | 
time step $\delta t=1200$secs. With this time step the stability  | 
| 140 | 
parameter to the horizontal laplacian friction \cite{Adcroft_thesis} | 
| 141 | 
 | 
| 142 | 
 | 
| 143 | 
 | 
| 144 | 
\begin{eqnarray} | 
| 145 | 
\label{EQ:laplacian_stability} | 
| 146 | 
S_{l} = 4 \frac{A_{h} \delta t}{{\Delta x}^2} | 
| 147 | 
\end{eqnarray} | 
| 148 | 
 | 
| 149 | 
\noindent evaluates to 0.012, which is well below the 0.3 upper limit | 
| 150 | 
for stability.  | 
| 151 | 
\\ | 
| 152 | 
 | 
| 153 | 
\noindent The numerical stability for inertial oscillations   | 
| 154 | 
\cite{Adcroft_thesis}  | 
| 155 | 
 | 
| 156 | 
\begin{eqnarray} | 
| 157 | 
\label{EQ:inertial_stability} | 
| 158 | 
S_{i} = f^{2} {\delta t}^2 | 
| 159 | 
\end{eqnarray} | 
| 160 | 
 | 
| 161 | 
\noindent evaluates to $0.0144$, which is well below the $0.5$ upper  | 
| 162 | 
limit for stability. | 
| 163 | 
\\ | 
| 164 | 
 | 
| 165 | 
\noindent The advective CFL \cite{Adcroft_thesis} for an extreme maximum  | 
| 166 | 
horizontal flow speed of $ | \vec{u} | = 2 ms^{-1}$ | 
| 167 | 
 | 
| 168 | 
\begin{eqnarray} | 
| 169 | 
\label{EQ:cfl_stability} | 
| 170 | 
S_{a} = \frac{| \vec{u} | \delta t}{ \Delta x} | 
| 171 | 
\end{eqnarray} | 
| 172 | 
 | 
| 173 | 
\noindent evaluates to 0.12. This is approaching the stability limit | 
| 174 | 
of 0.5 and limits $\delta t$ to $1200s$. | 
| 175 | 
 | 
| 176 | 
\subsection{Code Configuration} | 
| 177 | 
\label{SEC:code_config} | 
| 178 | 
 | 
| 179 | 
The model configuration for this experiment resides under the  | 
| 180 | 
directory {\it verification/exp0/}.  The experiment files  | 
| 181 | 
\begin{itemize} | 
| 182 | 
\item {\it input/data} | 
| 183 | 
\item {\it input/data.pkg} | 
| 184 | 
\item {\it input/eedata}, | 
| 185 | 
\item {\it input/windx.sin\_y}, | 
| 186 | 
\item {\it input/topog.box}, | 
| 187 | 
\item {\it code/CPP\_EEOPTIONS.h} | 
| 188 | 
\item {\it code/CPP\_OPTIONS.h}, | 
| 189 | 
\item {\it code/SIZE.h}.  | 
| 190 | 
\end{itemize} | 
| 191 | 
contain the code customisations and parameter settings for this  | 
| 192 | 
experiements. Below we describe the customisations | 
| 193 | 
to these files associated with this experiment. | 
| 194 | 
 | 
| 195 | 
\subsubsection{File {\it input/data}} | 
| 196 | 
 | 
| 197 | 
This file, reproduced completely below, specifies the main parameters  | 
| 198 | 
for the experiment. The parameters that are significant for this configuration | 
| 199 | 
are | 
| 200 | 
 | 
| 201 | 
\begin{itemize} | 
| 202 | 
 | 
| 203 | 
\item Line 7, \begin{verbatim} viscAh=4.E2, \end{verbatim} this line sets | 
| 204 | 
the laplacian friction coefficient to $400 m^2s^{-1}$ | 
| 205 | 
\item Line 10, \begin{verbatim} beta=1.E-11, \end{verbatim} this line sets | 
| 206 | 
$\beta$ (the gradient of the coriolis parameter, $f$) to $10^{-11} s^{-1}m^{-1}$ | 
| 207 | 
 | 
| 208 | 
\item Lines 15 and 16 | 
| 209 | 
\begin{verbatim} | 
| 210 | 
rigidLid=.FALSE., | 
| 211 | 
implicitFreeSurface=.TRUE., | 
| 212 | 
\end{verbatim} | 
| 213 | 
these lines suppress the rigid lid formulation of the surface | 
| 214 | 
pressure inverter and activate the implicit free surface form | 
| 215 | 
of the pressure inverter. | 
| 216 | 
 | 
| 217 | 
\item Line 27, | 
| 218 | 
\begin{verbatim} | 
| 219 | 
startTime=0, | 
| 220 | 
\end{verbatim} | 
| 221 | 
this line indicates that the experiment should start from $t=0$ | 
| 222 | 
and implicitly supresses searching for checkpoint files associated | 
| 223 | 
with restarting an numerical integration from a previously saved state. | 
| 224 | 
 | 
| 225 | 
\item Line 29, | 
| 226 | 
\begin{verbatim} | 
| 227 | 
endTime=12000, | 
| 228 | 
\end{verbatim} | 
| 229 | 
this line indicates that the experiment should start finish at $t=12000s$. | 
| 230 | 
A restart file will be written at this time that will enable the | 
| 231 | 
simulation to be continued from this point. | 
| 232 | 
 | 
| 233 | 
\item Line 30, | 
| 234 | 
\begin{verbatim} | 
| 235 | 
deltaTmom=1200, | 
| 236 | 
\end{verbatim} | 
| 237 | 
This line sets the momentum equation timestep to $1200s$. | 
| 238 | 
 | 
| 239 | 
\item Line 39, | 
| 240 | 
\begin{verbatim} | 
| 241 | 
usingCartesianGrid=.TRUE., | 
| 242 | 
\end{verbatim} | 
| 243 | 
This line requests that the simulation be performed in a  | 
| 244 | 
cartesian coordinate system. | 
| 245 | 
 | 
| 246 | 
\item Line 41, | 
| 247 | 
\begin{verbatim} | 
| 248 | 
delX=60*20E3, | 
| 249 | 
\end{verbatim} | 
| 250 | 
This line sets the horizontal grid spacing between each x-coordinate line | 
| 251 | 
in the discrete grid. The syntax indicates that the discrete grid | 
| 252 | 
should be comprise of $60$ grid lines each separated by $20 \times 10^{3}m$ | 
| 253 | 
($20$~km). | 
| 254 | 
 | 
| 255 | 
\item Line 42, | 
| 256 | 
\begin{verbatim} | 
| 257 | 
delY=60*20E3, | 
| 258 | 
\end{verbatim} | 
| 259 | 
This line sets the horizontal grid spacing between each y-coordinate line | 
| 260 | 
in the discrete grid to $20 \times 10^{3}m$ ($20$~km). | 
| 261 | 
 | 
| 262 | 
\item Line 43, | 
| 263 | 
\begin{verbatim} | 
| 264 | 
delZ=5000, | 
| 265 | 
\end{verbatim} | 
| 266 | 
This line sets the vertical grid spacing between each z-coordinate line | 
| 267 | 
in the discrete grid to $5000m$ ($5$~km). | 
| 268 | 
 | 
| 269 | 
\item Line 46, | 
| 270 | 
\begin{verbatim} | 
| 271 | 
bathyFile='topog.box' | 
| 272 | 
\end{verbatim} | 
| 273 | 
This line specifies the name of the file from which the domain | 
| 274 | 
bathymetry is read. This file is a two-dimensional ($x,y$) map of | 
| 275 | 
depths. This file is assumed to contain 64-bit binary numbers  | 
| 276 | 
giving the depth of the model at each grid cell, ordered with the x  | 
| 277 | 
coordinate varying fastest. The points are ordered from low coordinate | 
| 278 | 
to high coordinate for both axes. The units and orientation of the | 
| 279 | 
depths in this file are the same as used in the MITgcm code. In this | 
| 280 | 
experiment, a depth of $0m$ indicates a solid wall and a depth | 
| 281 | 
of $-5000m$ indicates open ocean. The matlab program | 
| 282 | 
{\it input/gendata.m} shows an example of how to generate a | 
| 283 | 
bathymetry file. | 
| 284 | 
 | 
| 285 | 
 | 
| 286 | 
\item Line 49, | 
| 287 | 
\begin{verbatim} | 
| 288 | 
zonalWindFile='windx.sin_y' | 
| 289 | 
\end{verbatim} | 
| 290 | 
This line specifies the name of the file from which the x-direction | 
| 291 | 
surface wind stress is read. This file is also a two-dimensional | 
| 292 | 
($x,y$) map and is enumerated and formatted in the same manner as the  | 
| 293 | 
bathymetry file. The matlab program {\it input/gendata.m} includes example  | 
| 294 | 
code to generate a valid {\bf zonalWindFile} file.   | 
| 295 | 
 | 
| 296 | 
\end{itemize} | 
| 297 | 
 | 
| 298 | 
\noindent other lines in the file {\it input/data} are standard values | 
| 299 | 
that are described in the MITgcm Getting Started and MITgcm Parameters | 
| 300 | 
notes. | 
| 301 | 
 | 
| 302 | 
\begin{small} | 
| 303 | 
\input{part3/case_studies/barotropic_gyre/input/data} | 
| 304 | 
\end{small} | 
| 305 | 
 | 
| 306 | 
\subsubsection{File {\it input/data.pkg}} | 
| 307 | 
 | 
| 308 | 
This file uses standard default values and does not contain | 
| 309 | 
customisations for this experiment. | 
| 310 | 
 | 
| 311 | 
\subsubsection{File {\it input/eedata}} | 
| 312 | 
 | 
| 313 | 
This file uses standard default values and does not contain | 
| 314 | 
customisations for this experiment. | 
| 315 | 
 | 
| 316 | 
\subsubsection{File {\it input/windx.sin\_y}} | 
| 317 | 
 | 
| 318 | 
The {\it input/windx.sin\_y} file specifies a two-dimensional ($x,y$)  | 
| 319 | 
map of wind stress ,$\tau_{x}$, values. The units used are $Nm^{-2}$. | 
| 320 | 
Although $\tau_{x}$ is only a function of $y$n in this experiment | 
| 321 | 
this file must still define a complete two-dimensional map in order | 
| 322 | 
to be compatible with the standard code for loading forcing fields  | 
| 323 | 
in MITgcm. The included matlab program {\it input/gendata.m} gives a complete | 
| 324 | 
code for creating the {\it input/windx.sin\_y} file. | 
| 325 | 
 | 
| 326 | 
\subsubsection{File {\it input/topog.box}} | 
| 327 | 
 | 
| 328 | 
 | 
| 329 | 
The {\it input/topog.box} file specifies a two-dimensional ($x,y$)  | 
| 330 | 
map of depth values. For this experiment values are either | 
| 331 | 
$0m$ or {\bf -delZ}m, corresponding respectively to a wall or to deep | 
| 332 | 
ocean. The file contains a raw binary stream of data that is enumerated | 
| 333 | 
in the same way as standard MITgcm two-dimensional, horizontal arrays. | 
| 334 | 
The included matlab program {\it input/gendata.m} gives a complete | 
| 335 | 
code for creating the {\it input/topog.box} file. | 
| 336 | 
 | 
| 337 | 
\subsubsection{File {\it code/SIZE.h}} | 
| 338 | 
 | 
| 339 | 
Two lines are customized in this file for the current experiment | 
| 340 | 
 | 
| 341 | 
\begin{itemize} | 
| 342 | 
 | 
| 343 | 
\item Line 39,  | 
| 344 | 
\begin{verbatim} sNx=60, \end{verbatim} this line sets | 
| 345 | 
the lateral domain extent in grid points for the | 
| 346 | 
axis aligned with the x-coordinate. | 
| 347 | 
 | 
| 348 | 
\item Line 40,  | 
| 349 | 
\begin{verbatim} sNy=60, \end{verbatim} this line sets | 
| 350 | 
the lateral domain extent in grid points for the | 
| 351 | 
axis aligned with the y-coordinate. | 
| 352 | 
 | 
| 353 | 
\end{itemize} | 
| 354 | 
 | 
| 355 | 
\begin{small} | 
| 356 | 
\input{part3/case_studies/barotropic_gyre/code/SIZE.h} | 
| 357 | 
\end{small} | 
| 358 | 
 | 
| 359 | 
\subsubsection{File {\it code/CPP\_OPTIONS.h}} | 
| 360 | 
 | 
| 361 | 
This file uses standard default values and does not contain | 
| 362 | 
customisations for this experiment. | 
| 363 | 
 | 
| 364 | 
 | 
| 365 | 
\subsubsection{File {\it code/CPP\_EEOPTIONS.h}} | 
| 366 | 
 | 
| 367 | 
This file uses standard default values and does not contain | 
| 368 | 
customisations for this experiment. | 
| 369 | 
 |