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