| 2 |
% $Name$ |
% $Name$ |
| 3 |
|
|
| 4 |
\section{Example: Four layer Baroclinic Ocean Gyre In Spherical Coordinates} |
\section{Example: Four layer Baroclinic Ocean Gyre In Spherical Coordinates} |
| 5 |
|
\label{sec:eg-fourlayer} |
| 6 |
|
|
| 7 |
\bodytext{bgcolor="#FFFFFFFF"} |
\bodytext{bgcolor="#FFFFFFFF"} |
| 8 |
|
|
| 16 |
%{\large May 2001} |
%{\large May 2001} |
| 17 |
%\end{center} |
%\end{center} |
| 18 |
|
|
| 19 |
\subsection{Introduction} |
This document describes an example experiment using MITgcm |
| 20 |
|
to simulate a baroclinic ocean gyre in spherical |
| 21 |
This document describes the second example MITgcm experiment. The first |
polar coordinates. The barotropic |
| 22 |
example experiment ilustrated how to configure the code for a single layer |
example experiment in section \ref{sec:eg-baro} |
| 23 |
|
ilustrated how to configure the code for a single layer |
| 24 |
simulation in a cartesian grid. In this example a similar physical problem |
simulation in a cartesian grid. In this example a similar physical problem |
| 25 |
is simulated, but the code is now configured |
is simulated, but the code is now configured |
| 26 |
for four layers and in a spherical polar coordinate system. |
for four layers and in a spherical polar coordinate system. |
| 37 |
In this experiment the model is configured to represent a mid-latitude |
In this experiment the model is configured to represent a mid-latitude |
| 38 |
enclosed sector of fluid on a sphere, $60^{\circ} \times 60^{\circ}$ in |
enclosed sector of fluid on a sphere, $60^{\circ} \times 60^{\circ}$ in |
| 39 |
lateral extent. The fluid is $2$~km deep and is forced |
lateral extent. The fluid is $2$~km deep and is forced |
| 40 |
by a constant in time zonal wind stress, $\tau_x$, that varies sinusoidally |
by a constant in time zonal wind stress, $\tau_{\lambda}$, that varies |
| 41 |
in the north-south direction. Topologically the simulated |
sinusoidally in the north-south direction. Topologically the simulated |
| 42 |
domain is a sector on a sphere and the coriolis parameter, $f$, is defined |
domain is a sector on a sphere and the coriolis parameter, $f$, is defined |
| 43 |
according to latitude, $\phi$ |
according to latitude, $\varphi$ |
| 44 |
|
|
| 45 |
\begin{equation} |
\begin{equation} |
| 46 |
\label{EQ:fcori} |
\label{EQ:fcori} |
| 47 |
f(\phi) = 2 \Omega \sin( \phi ) |
f(\varphi) = 2 \Omega \sin( \varphi ) |
| 48 |
\end{equation} |
\end{equation} |
| 49 |
|
|
| 50 |
\noindent with the rotation rate, $\Omega$ set to $\frac{2 \pi}{86400s}$. |
\noindent with the rotation rate, $\Omega$ set to $\frac{2 \pi}{86400s}$. |
| 54 |
|
|
| 55 |
\begin{equation} |
\begin{equation} |
| 56 |
\label{EQ:taux} |
\label{EQ:taux} |
| 57 |
\tau_x(\phi) = \tau_{0}\sin(\pi \frac{\phi}{L_{\phi}}) |
\tau_{\lambda}(\varphi) = \tau_{0}\sin(\pi \frac{\varphi}{L_{\varphi}}) |
| 58 |
\end{equation} |
\end{equation} |
| 59 |
|
|
| 60 |
\noindent where $L_{\phi}$ is the lateral domain extent ($60^{\circ}$) and |
\noindent where $L_{\varphi}$ is the lateral domain extent ($60^{\circ}$) and |
| 61 |
$\tau_0$ is set to $0.1N m^{-2}$. |
$\tau_0$ is set to $0.1N m^{-2}$. |
| 62 |
\\ |
\\ |
| 63 |
|
|
| 64 |
Figure \ref{FIG:simulation_config} |
Figure \ref{FIG:simulation_config} |
| 65 |
summarises the configuration simulated. |
summarises the configuration simulated. |
| 66 |
In contrast to example (1) \cite{baro_gyre_case_study}, the current |
In contrast to the example in section \ref{sec:eg-baro}, the |
| 67 |
experiment simulates a spherical polar domain. However, as indicated |
current experiment simulates a spherical polar domain. As indicated |
| 68 |
by the axes in the lower left of the figure the model code works internally |
by the axes in the lower left of the figure the model code works internally |
| 69 |
in a locally orthoganal coordinate $(x,y,z)$. In the remainder of this |
in a locally orthoganal coordinate $(x,y,z)$. For this experiment description |
| 70 |
document the local coordinate $(x,y,z)$ will be adopted. |
of this document the local orthogonal model coordinate $(x,y,z)$ is synonomous |
| 71 |
|
with the spherical polar coordinate shown in figure |
| 72 |
|
\ref{fig:spherical-polar-coord} |
| 73 |
\\ |
\\ |
| 74 |
|
|
| 75 |
The experiment has four levels in the vertical, each of equal thickness, |
The experiment has four levels in the vertical, each of equal thickness, |
| 95 |
|
|
| 96 |
\noindent with $\rho_{0}=999.8\,{\rm kg\,m}^{-3}$ and |
\noindent with $\rho_{0}=999.8\,{\rm kg\,m}^{-3}$ and |
| 97 |
$\alpha_{\theta}=2\times10^{-4}\,{\rm degrees}^{-1} $. Integrated forward in |
$\alpha_{\theta}=2\times10^{-4}\,{\rm degrees}^{-1} $. Integrated forward in |
| 98 |
this configuration the model state variable {\bf theta} is synonomous with |
this configuration the model state variable {\bf theta} is equivalent to |
| 99 |
either in-situ temperature, $T$, or potential temperature, $\theta$. For |
either in-situ temperature, $T$, or potential temperature, $\theta$. For |
| 100 |
consistency with later examples, in which the equation of state is |
consistency with later examples, in which the equation of state is |
| 101 |
non-linear, we use $\theta$ to represent temperature here. This is |
non-linear, we use $\theta$ to represent temperature here. This is |
| 110 |
\caption{Schematic of simulation domain and wind-stress forcing function |
\caption{Schematic of simulation domain and wind-stress forcing function |
| 111 |
for the four-layer gyre numerical experiment. The domain is enclosed by solid |
for the four-layer gyre numerical experiment. The domain is enclosed by solid |
| 112 |
walls at $0^{\circ}$~E, $60^{\circ}$~E, $0^{\circ}$~N and $60^{\circ}$~N. |
walls at $0^{\circ}$~E, $60^{\circ}$~E, $0^{\circ}$~N and $60^{\circ}$~N. |
| 113 |
In the four-layer case an initial temperature stratification is |
An initial stratification is |
| 114 |
imposed by setting the potential temperature, $\theta$, in each layer. |
imposed by setting the potential temperature, $\theta$, in each layer. |
| 115 |
The vertical spacing, $\Delta z$, is constant and equal to $500$m. |
The vertical spacing, $\Delta z$, is constant and equal to $500$m. |
| 116 |
} |
} |
| 117 |
\label{FIG:simulation_config} |
\label{FIG:simulation_config} |
| 118 |
\end{figure} |
\end{figure} |
| 119 |
|
|
| 120 |
\subsection{Discrete Numerical Configuration} |
\subsection{Equations solved} |
| 121 |
|
|
| 122 |
The model is configured in hydrostatic form. The domain is discretised with |
The implicit free surface {\bf HPE} form of the |
| 123 |
a uniform grid spacing in latitude and longitude |
equations described in Marshall et. al \cite{Marshall97a} is |
| 124 |
$\Delta x=\Delta y=1^{\circ}$, so |
employed. The flow is three-dimensional with just temperature, $\theta$, as |
| 125 |
that there are sixty grid cells in the $x$ and $y$ directions. Vertically the |
an active tracer. The equation of state is linear. |
|
model is configured with a four layers with constant depth, |
|
|
$\Delta z$, of $500$~m. |
|
|
The implicit free surface form of the |
|
|
pressure equation described in Marshall et. al \cite{Marshall97a} is |
|
|
employed. |
|
| 126 |
A horizontal laplacian operator $\nabla_{h}^2$ provides viscous |
A horizontal laplacian operator $\nabla_{h}^2$ provides viscous |
| 127 |
dissipation. The wind-stress momentum input is added to the momentum equation |
dissipation and provides a diffusive sub-grid scale closure for the |
| 128 |
for the ``zonal flow'', $u$. Other terms in the model |
temperature equation. A wind-stress momentum forcing is added to the momentum |
| 129 |
|
equation for the zonal flow, $u$. Other terms in the model |
| 130 |
are explicitly switched off for this experiement configuration (see section |
are explicitly switched off for this experiement configuration (see section |
| 131 |
\ref{SEC:code_config} ), yielding an active set of equations solved in this |
\ref{SEC:eg_fourl_code_config} ). This yields an active set of equations |
| 132 |
configuration as follows |
solved in this configuration, written in spherical polar coordinates as |
| 133 |
|
follows |
| 134 |
|
|
| 135 |
\begin{eqnarray} |
\begin{eqnarray} |
| 136 |
\label{EQ:model_equations} |
\label{EQ:model_equations} |
| 137 |
\frac{Du}{Dt} - fv + |
\frac{Du}{Dt} - fv + |
| 138 |
\frac{1}{\rho}\frac{\partial p^{'}}{\partial x} - |
\frac{1}{\rho}\frac{\partial p^{\prime}}{\partial \lambda} - |
| 139 |
A_{h}\nabla_{h}^2u - A_{z}\frac{\partial^{2}u}{\partial z^{2}} |
A_{h}\nabla_{h}^2u - A_{z}\frac{\partial^{2}u}{\partial z^{2}} |
| 140 |
& = & |
& = & |
| 141 |
\cal{F} |
\cal{F}_{\lambda} |
| 142 |
\\ |
\\ |
| 143 |
\frac{Dv}{Dt} + fu + |
\frac{Dv}{Dt} + fu + |
| 144 |
\frac{1}{\rho}\frac{\partial p^{'}}{\partial y} - |
\frac{1}{\rho}\frac{\partial p^{\prime}}{\partial \varphi} - |
| 145 |
A_{h}\nabla_{h}^2v - A_{z}\frac{\partial^{2}v}{\partial z^{2}} |
A_{h}\nabla_{h}^2v - A_{z}\frac{\partial^{2}v}{\partial z^{2}} |
| 146 |
& = & |
& = & |
| 147 |
0 |
0 |
| 148 |
\\ |
\\ |
| 149 |
\frac{\partial \eta}{\partial t} + \nabla_{h}\cdot \vec{u} |
\frac{\partial \eta}{\partial t} + \frac{\partial H \widehat{u}}{\partial \lambda} + |
| 150 |
|
\frac{\partial H \widehat{v}}{\partial \varphi} |
| 151 |
&=& |
&=& |
| 152 |
0 |
0 |
| 153 |
|
\label{eq:fourl_example_continuity} |
| 154 |
\\ |
\\ |
| 155 |
\frac{D\theta}{Dt} - |
\frac{D\theta}{Dt} - |
| 156 |
K_{h}\nabla_{h}^2\theta - K_{z}\frac{\partial^{2}\theta}{\partial z^{2}} |
K_{h}\nabla_{h}^2\theta - K_{z}\frac{\partial^{2}\theta}{\partial z^{2}} |
| 157 |
& = & |
& = & |
| 158 |
0 |
0 |
| 159 |
|
\label{eq:eg_fourl_theta} |
| 160 |
|
\\ |
| 161 |
|
p^{\prime} & = & |
| 162 |
|
g\rho_{0} \eta + \int^{0}_{-z}\rho^{\prime} dz |
| 163 |
\\ |
\\ |
| 164 |
g\rho_{0} \eta + \int^{0}_{-z}\rho^{'} dz & = & p^{'} |
\rho^{\prime} & = &- \alpha_{\theta}\rho_{0}\theta^{\prime} |
| 165 |
\\ |
\\ |
| 166 |
{\cal F} |_{s} & = & \frac{\tau_{x}}{\rho_{0}\Delta z_{s}} |
{\cal F}_{\lambda} |_{s} & = & \frac{\tau_{\lambda}}{\rho_{0}\Delta z_{s}} |
| 167 |
\\ |
\\ |
| 168 |
{\cal F} |_{i} & = & 0 |
{\cal F}_{\lambda} |_{i} & = & 0 |
| 169 |
\end{eqnarray} |
\end{eqnarray} |
| 170 |
|
|
| 171 |
\noindent where $u$ and $v$ are the $x$ and $y$ components of the |
\noindent where $u$ and $v$ are the components of the horizontal |
| 172 |
flow vector $\vec{u}$. The suffices ${s},{i}$ indicate surface and |
flow vector $\vec{u}$ on the sphere ($u=\dot{\lambda},v=\dot{\varphi}$). |
| 173 |
interior model levels respectively. As described in |
The terms $H\widehat{u}$ and $H\widehat{v}$ are the components of the vertical |
| 174 |
MITgcm Numerical Solution Procedure \cite{MITgcm_Numerical_Scheme}, the time |
integral term given in equation \ref{eq:free-surface} and |
| 175 |
evolution of potential temperature, $\theta$, equation is solved prognostically. |
explained in more detail in section \ref{sect:pressure-method-linear-backward}. |
| 176 |
The total pressure, $p$, is diagnosed by summing pressure due to surface |
However, for the problem presented here, the continuity relation (equation |
| 177 |
elevation $\eta$ and the hydrostatic pressure. |
\ref{eq:fourl_example_continuity}) differs from the general form given |
| 178 |
\\ |
in section \ref{sect:pressure-method-linear-backward}, |
| 179 |
|
equation \ref{eq:linear-free-surface=P-E+R}, because the source terms |
| 180 |
|
${\cal P}-{\cal E}+{\cal R}$ |
| 181 |
|
are all $0$. |
| 182 |
|
|
| 183 |
|
The pressure field, $p^{\prime}$, is separated into a barotropic part |
| 184 |
|
due to variations in sea-surface height, $\eta$, and a hydrostatic |
| 185 |
|
part due to variations in density, $\rho^{\prime}$, integrated |
| 186 |
|
through the water column. |
| 187 |
|
|
| 188 |
|
The suffices ${s},{i}$ indicate surface and interior of the domain. |
| 189 |
|
The windstress forcing, ${\cal F}_{\lambda}$, is applied in the surface layer |
| 190 |
|
by a source term in the zonal momentum equation. In the ocean interior |
| 191 |
|
this term is zero. |
| 192 |
|
|
| 193 |
|
In the momentum equations |
| 194 |
|
lateral and vertical boundary conditions for the $\nabla_{h}^{2}$ |
| 195 |
|
and $\frac{\partial^{2}}{\partial z^{2}}$ operators are specified |
| 196 |
|
when the numerical simulation is run - see section |
| 197 |
|
\ref{SEC:eg_fourl_code_config}. For temperature |
| 198 |
|
the boundary condition is ``zero-flux'' |
| 199 |
|
e.g. $\frac{\partial \theta}{\partial \varphi}= |
| 200 |
|
\frac{\partial \theta}{\partial \lambda}=\frac{\partial \theta}{\partial z}=0$. |
| 201 |
|
|
| 202 |
|
|
| 203 |
|
|
| 204 |
|
\subsection{Discrete Numerical Configuration} |
| 205 |
|
|
| 206 |
|
The domain is discretised with |
| 207 |
|
a uniform grid spacing in latitude and longitude |
| 208 |
|
$\Delta \lambda=\Delta \varphi=1^{\circ}$, so |
| 209 |
|
that there are sixty grid cells in the zonal and meridional directions. |
| 210 |
|
Vertically the |
| 211 |
|
model is configured with four layers with constant depth, |
| 212 |
|
$\Delta z$, of $500$~m. The internal, locally orthogonal, model coordinate |
| 213 |
|
variables $x$ and $y$ are initialised from the values of |
| 214 |
|
$\lambda$, $\varphi$, $\Delta \lambda$ and $\Delta \varphi$ in |
| 215 |
|
radians according to |
| 216 |
|
|
| 217 |
|
\begin{eqnarray} |
| 218 |
|
x=r\cos(\varphi)\lambda,~\Delta x & = &r\cos(\varphi)\Delta \lambda \\ |
| 219 |
|
y=r\varphi,~\Delta y &= &r\Delta \varphi |
| 220 |
|
\end{eqnarray} |
| 221 |
|
|
| 222 |
|
The procedure for generating a set of internal grid variables from a |
| 223 |
|
spherical polar grid specification is discussed in section |
| 224 |
|
\ref{sec:spatial_discrete_horizontal_grid}. |
| 225 |
|
|
| 226 |
|
\noindent\fbox{ \begin{minipage}{5.5in} |
| 227 |
|
{\em S/R INI\_SPHERICAL\_POLAR\_GRID} ({\em |
| 228 |
|
model/src/ini\_spherical\_polar\_grid.F}) |
| 229 |
|
|
| 230 |
|
$A_c$, $A_\zeta$, $A_w$, $A_s$: {\bf rAc}, {\bf rAz}, {\bf rAw}, {\bf rAs} |
| 231 |
|
({\em GRID.h}) |
| 232 |
|
|
| 233 |
|
$\Delta x_g$, $\Delta y_g$: {\bf DXg}, {\bf DYg} ({\em GRID.h}) |
| 234 |
|
|
| 235 |
|
$\Delta x_c$, $\Delta y_c$: {\bf DXc}, {\bf DYc} ({\em GRID.h}) |
| 236 |
|
|
| 237 |
|
$\Delta x_f$, $\Delta y_f$: {\bf DXf}, {\bf DYf} ({\em GRID.h}) |
| 238 |
|
|
| 239 |
|
$\Delta x_v$, $\Delta y_u$: {\bf DXv}, {\bf DYu} ({\em GRID.h}) |
| 240 |
|
|
| 241 |
|
\end{minipage} }\\ |
| 242 |
|
|
| 243 |
|
|
| 244 |
|
|
| 245 |
|
As described in \ref{sec:tracer_equations}, the time evolution of potential |
| 246 |
|
temperature, |
| 247 |
|
$\theta$, (equation \ref{eq:eg_fourl_theta}) |
| 248 |
|
is evaluated prognostically. The centered second-order scheme with |
| 249 |
|
Adams-Bashforth time stepping described in section |
| 250 |
|
\ref{sec:tracer_equations_abII} is used to step forward the temperature |
| 251 |
|
equation. The pressure forces that drive the fluid motions, ( |
| 252 |
|
$\frac{\partial p^{'}}{\partial \lambda}$ and $\frac{\partial p^{'}}{\partial \varphi}$), are found by summing pressure due to surface |
| 253 |
|
elevation $\eta$ and the hydrostatic pressure. The hydrostatic part of the |
| 254 |
|
pressure is evaluated explicitly by integrating density. The sea-surface |
| 255 |
|
height, $\eta$, is solved for implicitly as described in section |
| 256 |
|
\ref{sect:pressure-method-linear-backward}. |
| 257 |
|
|
| 258 |
\subsubsection{Numerical Stability Criteria} |
\subsubsection{Numerical Stability Criteria} |
| 259 |
|
|
| 334 |
stability limit of 0.25. |
stability limit of 0.25. |
| 335 |
|
|
| 336 |
\subsection{Code Configuration} |
\subsection{Code Configuration} |
| 337 |
\label{SEC:code_config} |
\label{SEC:eg_fourl_code_config} |
| 338 |
|
|
| 339 |
The model configuration for this experiment resides under the |
The model configuration for this experiment resides under the |
| 340 |
directory {\it verification/exp1/}. The experiment files |
directory {\it verification/exp1/}. The experiment files |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch