1 |
% $Header: /u/gcmpack/mitgcmdoc/part1/manual.tex,v 1.3 2001/10/10 16:48:45 cnh Exp $ |
2 |
% $Name: $ |
3 |
|
4 |
\documentclass[12pt]{book} |
5 |
\usepackage{amsmath} |
6 |
\usepackage{html} |
7 |
\usepackage{epsfig} |
8 |
\usepackage{graphics,subfigure} |
9 |
\usepackage{array} |
10 |
\usepackage{multirow} |
11 |
\usepackage{fancyhdr} |
12 |
\usepackage{psfrag} |
13 |
|
14 |
%TCIDATA{OutputFilter=Latex.dll} |
15 |
%TCIDATA{LastRevised=Thursday, October 04, 2001 14:41:22} |
16 |
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">} |
17 |
%TCIDATA{Language=American English} |
18 |
|
19 |
\fancyhead{} |
20 |
\fancyhead[LO]{\slshape \rightmark} |
21 |
\fancyhead[RE]{\slshape \leftmark} |
22 |
\fancyhead[RO,LE]{\thepage} |
23 |
\fancyfoot[CO,CE]{\today} |
24 |
\fancyfoot[RO,LE]{ } |
25 |
\renewcommand{\headrulewidth}{0.4pt} |
26 |
\renewcommand{\footrulewidth}{0.4pt} |
27 |
\setcounter{secnumdepth}{3} |
28 |
\input{tcilatex} |
29 |
|
30 |
\begin{document} |
31 |
|
32 |
\tableofcontents |
33 |
|
34 |
|
35 |
\part{MIT GCM basics} |
36 |
|
37 |
% Section: Overview |
38 |
|
39 |
% $Header: /u/gcmpack/mitgcmdoc/part1/manual.tex,v 1.3 2001/10/10 16:48:45 cnh Exp $ |
40 |
% $Name: $ |
41 |
|
42 |
\section{Introduction} |
43 |
|
44 |
This documentation provides the reader with the information necessary to |
45 |
carry out numerical experiments using MITgcm. It gives a comprehensive |
46 |
description of the continuous equations on which the model is based, the |
47 |
numerical algorithms the model employs and a description of the associated |
48 |
program code. Along with the hydrodynamical kernel, physical and |
49 |
biogeochemical parameterizations of key atmospheric and oceanic processes |
50 |
are available. A number of examples illustrating the use of the model in |
51 |
both process and general circulation studies of the atmosphere and ocean are |
52 |
also presented. |
53 |
|
54 |
MITgcm has a number of novel aspects: |
55 |
|
56 |
\begin{itemize} |
57 |
\item it can be used to study both atmospheric and oceanic phenomena; one |
58 |
hydrodynamical kernel is used to drive forward both atmospheric and oceanic |
59 |
models - see fig |
60 |
\marginpar{ |
61 |
Fig.1 One model}\ref{fig:onemodel} |
62 |
|
63 |
%% CNHbegin |
64 |
%notci%\input{part1/one_model_figure} |
65 |
%% CNHend |
66 |
|
67 |
\item it has a non-hydrostatic capability and so can be used to study both |
68 |
small-scale and large scale processes - see fig |
69 |
\marginpar{ |
70 |
Fig.2 All scales}\ref{fig:all-scales} |
71 |
|
72 |
%% CNHbegin |
73 |
%notci%\input{part1/all_scales_figure} |
74 |
%% CNHend |
75 |
|
76 |
\item finite volume techniques are employed yielding an intuitive |
77 |
discretization and support for the treatment of irregular geometries using |
78 |
orthogonal curvilinear grids and shaved cells - see fig |
79 |
\marginpar{ |
80 |
Fig.3 Finite volumes}\ref{fig:finite-volumes} |
81 |
|
82 |
%% CNHbegin |
83 |
%notci%\input{part1/fvol_figure} |
84 |
%% CNHend |
85 |
|
86 |
\item tangent linear and adjoint counterparts are automatically maintained |
87 |
along with the forward model, permitting sensitivity and optimization |
88 |
studies. |
89 |
|
90 |
\item the model is developed to perform efficiently on a wide variety of |
91 |
computational platforms. |
92 |
\end{itemize} |
93 |
|
94 |
Key publications reporting on and charting the development of the model are |
95 |
listed in an Appendix. |
96 |
|
97 |
We begin by briefly showing some of the results of the model in action to |
98 |
give a feel for the wide range of problems that can be addressed using it. |
99 |
\pagebreak |
100 |
|
101 |
% $Header: /u/gcmpack/mitgcmdoc/part1/manual.tex,v 1.3 2001/10/10 16:48:45 cnh Exp $ |
102 |
% $Name: $ |
103 |
|
104 |
\section{Illustrations of the model in action} |
105 |
|
106 |
The MITgcm has been designed and used to model a wide range of phenomena, |
107 |
from convection on the scale of meters in the ocean to the global pattern of |
108 |
atmospheric winds - see fig.2\ref{fig:all-scales}. To give a flavor of the |
109 |
kinds of problems the model has been used to study, we briefly describe some |
110 |
of them here. A more detailed description of the underlying formulation, |
111 |
numerical algorithm and implementation that lie behind these calculations is |
112 |
given later. Indeed many of the illustrative examples shown below can be |
113 |
easily reproduced: simply download the model (the minimum you need is a PC |
114 |
running linux, together with a FORTRAN\ 77 compiler) and follow the examples |
115 |
described in detail in the documentation. |
116 |
|
117 |
\subsection{Global atmosphere: `Held-Suarez' benchmark} |
118 |
|
119 |
A novel feature of MITgcm is its ability to simulate both atmospheric and |
120 |
oceanographic flows at both small and large scales. |
121 |
|
122 |
Fig.E1a.\ref{fig:eddy_cs} shows an instantaneous plot of the 500$mb$ |
123 |
temperature field obtained using the atmospheric isomorph of MITgcm run at |
124 |
2.8$^{\circ }$ resolution on the cubed sphere. We see cold air over the pole |
125 |
(blue) and warm air along an equatorial band (red). Fully developed |
126 |
baroclinic eddies spawned in the northern hemisphere storm track are |
127 |
evident. There are no mountains or land-sea contrast in this calculation, |
128 |
but you can easily put them in. The model is driven by relaxation to a |
129 |
radiative-convective equilibrium profile, following the description set out |
130 |
in Held and Suarez; 1994 designed to test atmospheric hydrodynamical cores - |
131 |
there are no mountains or land-sea contrast. |
132 |
|
133 |
%% CNHbegin |
134 |
%notci%\input{part1/cubic_eddies_figure} |
135 |
%% CNHend |
136 |
|
137 |
As described in Adcroft (2001), a `cubed sphere' is used to discretize the |
138 |
globe permitting a uniform gridding and obviated the need to fourier filter. |
139 |
The `vector-invariant' form of MITgcm supports any orthogonal curvilinear |
140 |
grid, of which the cubed sphere is just one of many choices. |
141 |
|
142 |
Fig.E1b shows the 5-year mean, zonally averaged potential temperature, zonal |
143 |
wind and meridional overturning streamfunction from a 20-level version of |
144 |
the model. It compares favorable with more conventional spatial |
145 |
discretization approaches. |
146 |
|
147 |
A regular spherical lat-lon grid can also be used. |
148 |
|
149 |
%% CNHbegin |
150 |
%notci%\input{part1/hs_zave_u_figure} |
151 |
%% CNHend |
152 |
|
153 |
\subsection{Ocean gyres} |
154 |
|
155 |
Baroclinic instability is a ubiquitous process in the ocean, as well as the |
156 |
atmosphere. Ocean eddies play an important role in modifying the |
157 |
hydrographic structure and current systems of the oceans. Coarse resolution |
158 |
models of the oceans cannot resolve the eddy field and yield rather broad, |
159 |
diffusive patterns of ocean currents. But if the resolution of our models is |
160 |
increased until the baroclinic instability process is resolved, numerical |
161 |
solutions of a different and much more realistic kind, can be obtained. |
162 |
|
163 |
Fig. ?.? shows the surface temperature and velocity field obtained from |
164 |
MITgcm run at $\frac{1}{6}^{\circ }$ horizontal resolution on a $lat-lon$ |
165 |
grid in which the pole has been rotated by 90$^{\circ }$ on to the equator |
166 |
(to avoid the converging of meridian in northern latitudes). 21 vertical |
167 |
levels are used in the vertical with a `lopped cell' representation of |
168 |
topography. The development and propagation of anomalously warm and cold |
169 |
eddies can be clearly been seen in the Gulf Stream region. The transport of |
170 |
warm water northward by the mean flow of the Gulf Stream is also clearly |
171 |
visible. |
172 |
|
173 |
%% CNHbegin |
174 |
%notci%\input{part1/ocean_gyres_figure} |
175 |
%% CNHend |
176 |
|
177 |
|
178 |
\subsection{Global ocean circulation} |
179 |
|
180 |
Fig.E2a shows the pattern of ocean currents at the surface of a 4$^{\circ }$ |
181 |
global ocean model run with 15 vertical levels. Lopped cells are used to |
182 |
represent topography on a regular $lat-lon$ grid extending from 70$^{\circ |
183 |
}N $ to 70$^{\circ }S$. The model is driven using monthly-mean winds with |
184 |
mixed boundary conditions on temperature and salinity at the surface. The |
185 |
transfer properties of ocean eddies, convection and mixing is parameterized |
186 |
in this model. |
187 |
|
188 |
Fig.E2b shows the meridional overturning circulation of the global ocean in |
189 |
Sverdrups. |
190 |
|
191 |
%%CNHbegin |
192 |
%notci%\input{part1/global_circ_figure} |
193 |
%%CNHend |
194 |
|
195 |
\subsection{Convection and mixing over topography} |
196 |
|
197 |
Dense plumes generated by localized cooling on the continental shelf of the |
198 |
ocean may be influenced by rotation when the deformation radius is smaller |
199 |
than the width of the cooling region. Rather than gravity plumes, the |
200 |
mechanism for moving dense fluid down the shelf is then through geostrophic |
201 |
eddies. The simulation shown in the figure (blue is cold dense fluid, red is |
202 |
warmer, lighter fluid) employs the non-hydrostatic capability of MITgcm to |
203 |
trigger convection by surface cooling. The cold, dense water falls down the |
204 |
slope but is deflected along the slope by rotation. It is found that |
205 |
entrainment in the vertical plane is reduced when rotational control is |
206 |
strong, and replaced by lateral entrainment due to the baroclinic |
207 |
instability of the along-slope current. |
208 |
|
209 |
%%CNHbegin |
210 |
%notci%\input{part1/convect_and_topo} |
211 |
%%CNHend |
212 |
|
213 |
\subsection{Boundary forced internal waves} |
214 |
|
215 |
The unique ability of MITgcm to treat non-hydrostatic dynamics in the |
216 |
presence of complex geometry makes it an ideal tool to study internal wave |
217 |
dynamics and mixing in oceanic canyons and ridges driven by large amplitude |
218 |
barotropic tidal currents imposed through open boundary conditions. |
219 |
|
220 |
Fig. ?.? shows the influence of cross-slope topographic variations on |
221 |
internal wave breaking - the cross-slope velocity is in color, the density |
222 |
contoured. The internal waves are excited by application of open boundary |
223 |
conditions on the left.\ They propagate to the sloping boundary (represented |
224 |
using MITgcm's finite volume spatial discretization) where they break under |
225 |
nonhydrostatic dynamics. |
226 |
|
227 |
%%CNHbegin |
228 |
%notci%\input{part1/boundary_forced_waves} |
229 |
%%CNHend |
230 |
|
231 |
\subsection{Parameter sensitivity using the adjoint of MITgcm} |
232 |
|
233 |
Forward and tangent linear counterparts of MITgcm are supported using an |
234 |
`automatic adjoint compiler'. These can be used in parameter sensitivity and |
235 |
data assimilation studies. |
236 |
|
237 |
As one example of application of the MITgcm adjoint, Fig.E4 maps the |
238 |
gradient $\frac{\partial J}{\partial \mathcal{H}}$where $J$ is the magnitude |
239 |
of the overturning streamfunction shown in fig?.? at 40$^{\circ }$N and $ |
240 |
\mathcal{H}$ is the air-sea heat flux 100 years before. We see that $J$ is |
241 |
sensitive to heat fluxes over the Labrador Sea, one of the important sources |
242 |
of deep water for the thermohaline circulations. This calculation also |
243 |
yields sensitivities to all other model parameters. |
244 |
|
245 |
%%CNHbegin |
246 |
%notci%\input{part1/adj_hf_ocean_figure} |
247 |
%%CNHend |
248 |
|
249 |
\subsection{Global state estimation of the ocean} |
250 |
|
251 |
An important application of MITgcm is in state estimation of the global |
252 |
ocean circulation. An appropriately defined `cost function', which measures |
253 |
the departure of the model from observations (both remotely sensed and |
254 |
insitu) over an interval of time, is minimized by adjusting `control |
255 |
parameters' such as air-sea fluxes, the wind field, the initial conditions |
256 |
etc. Figure ?.? shows an estimate of the time-mean surface elevation of the |
257 |
ocean obtained by bringing the model in to consistency with altimetric and |
258 |
in-situ observations over the period 1992-1997. |
259 |
|
260 |
%% CNHbegin |
261 |
%notci%\input{part1/globes_figure} |
262 |
%% CNHend |
263 |
|
264 |
\subsection{Ocean biogeochemical cycles} |
265 |
|
266 |
MITgcm is being used to study global biogeochemical cycles in the ocean. For |
267 |
example one can study the effects of interannual changes in meteorological |
268 |
forcing and upper ocean circulation on the fluxes of carbon dioxide and |
269 |
oxygen between the ocean and atmosphere. The figure shows the annual air-sea |
270 |
flux of oxygen and its relation to density outcrops in the southern oceans |
271 |
from a single year of a global, interannually varying simulation. |
272 |
|
273 |
%%CNHbegin |
274 |
%notci%\input{part1/biogeo_figure} |
275 |
%%CNHend |
276 |
|
277 |
\subsection{Simulations of laboratory experiments} |
278 |
|
279 |
Figure ?.? shows MITgcm being used to simulate a laboratory experiment |
280 |
enquiring in to the dynamics of the Antarctic Circumpolar Current (ACC). An |
281 |
initially homogeneous tank of water ($1m$ in diameter) is driven from its |
282 |
free surface by a rotating heated disk. The combined action of mechanical |
283 |
and thermal forcing creates a lens of fluid which becomes baroclinically |
284 |
unstable. The stratification and depth of penetration of the lens is |
285 |
arrested by its instability in a process analogous to that whic sets the |
286 |
stratification of the ACC. |
287 |
|
288 |
%%CNHbegin |
289 |
%notci%\input{part1/lab_figure} |
290 |
%%CNHend |
291 |
|
292 |
% $Header: /u/gcmpack/mitgcmdoc/part1/manual.tex,v 1.3 2001/10/10 16:48:45 cnh Exp $ |
293 |
% $Name: $ |
294 |
|
295 |
\section{Continuous equations in `r' coordinates} |
296 |
|
297 |
To render atmosphere and ocean models from one dynamical core we exploit |
298 |
`isomorphisms' between equation sets that govern the evolution of the |
299 |
respective fluids - see fig.4 |
300 |
\marginpar{ |
301 |
Fig.4. Isomorphisms}. One system of hydrodynamical equations is written down |
302 |
and encoded. The model variables have different interpretations depending on |
303 |
whether the atmosphere or ocean is being studied. Thus, for example, the |
304 |
vertical coordinate `$r$' is interpreted as pressure, $p$, if we are |
305 |
modeling the atmosphere and height, $z$, if we are modeling the ocean. |
306 |
|
307 |
%%CNHbegin |
308 |
%notci%\input{part1/zandpcoord_figure.tex} |
309 |
%%CNHend |
310 |
|
311 |
The state of the fluid at any time is characterized by the distribution of |
312 |
velocity $\vec{\mathbf{v}}$, active tracers $\theta $ and $S$, a |
313 |
`geopotential' $\phi $ and density $\rho =\rho (\theta ,S,p)$ which may |
314 |
depend on $\theta $, $S$, and $p$. The equations that govern the evolution |
315 |
of these fields, obtained by applying the laws of classical mechanics and |
316 |
thermodynamics to a Boussinesq, Navier-Stokes fluid are, written in terms of |
317 |
a generic vertical coordinate, $r$, see fig.5 |
318 |
\marginpar{ |
319 |
Fig.5 The vertical coordinate of model}: |
320 |
|
321 |
%%CNHbegin |
322 |
%notci%\input{part1/vertcoord_figure.tex} |
323 |
%%CNHend |
324 |
|
325 |
\begin{equation*} |
326 |
\frac{D\vec{\mathbf{v}_{h}}}{Dt}+\left( 2\vec{\Omega}\times \vec{\mathbf{v}} |
327 |
\right) _{h}+\mathbf{\nabla }_{h}\phi =\mathcal{F}_{\vec{\mathbf{v}_{h}}} |
328 |
\text{ horizontal mtm} |
329 |
\end{equation*} |
330 |
|
331 |
\begin{equation*} |
332 |
\frac{D\dot{r}}{Dt}+\widehat{k}\cdot \left( 2\vec{\Omega}\times \vec{\mathbf{ |
333 |
v}}\right) +\frac{\partial \phi }{\partial r}+b=\mathcal{F}_{\dot{r}}\text{ |
334 |
vertical mtm} |
335 |
\end{equation*} |
336 |
|
337 |
\begin{equation} |
338 |
\mathbf{\nabla }_{h}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial \dot{r}}{ |
339 |
\partial r}=0\text{ continuity} \label{eq:continuous} |
340 |
\end{equation} |
341 |
|
342 |
\begin{equation*} |
343 |
b=b(\theta ,S,r)\text{ equation of state} |
344 |
\end{equation*} |
345 |
|
346 |
\begin{equation*} |
347 |
\frac{D\theta }{Dt}=\mathcal{Q}_{\theta }\text{ potential temperature} |
348 |
\end{equation*} |
349 |
|
350 |
\begin{equation*} |
351 |
\frac{DS}{Dt}=\mathcal{Q}_{S}\text{ humidity/salinity} |
352 |
\end{equation*} |
353 |
|
354 |
Here: |
355 |
|
356 |
\begin{equation*} |
357 |
r\text{ is the vertical coordinate} |
358 |
\end{equation*} |
359 |
|
360 |
\begin{equation*} |
361 |
\frac{D}{Dt}=\frac{\partial }{\partial t}+\vec{\mathbf{v}}\cdot \nabla \text{ |
362 |
is the total derivative} |
363 |
\end{equation*} |
364 |
|
365 |
\begin{equation*} |
366 |
\mathbf{\nabla }=\mathbf{\nabla }_{h}+\widehat{k}\frac{\partial }{\partial r} |
367 |
\text{ is the `grad' operator} |
368 |
\end{equation*} |
369 |
with $\mathbf{\nabla }_{h}$ operating in the horizontal and $\widehat{k} |
370 |
\frac{\partial }{\partial r}$ operating in the vertical, where $\widehat{k}$ |
371 |
is a unit vector in the vertical |
372 |
|
373 |
\begin{equation*} |
374 |
t\text{ is time} |
375 |
\end{equation*} |
376 |
|
377 |
\begin{equation*} |
378 |
\vec{\mathbf{v}}=(u,v,\dot{r})=(\vec{\mathbf{v}}_{h},\dot{r})\text{ is the |
379 |
velocity} |
380 |
\end{equation*} |
381 |
|
382 |
\begin{equation*} |
383 |
\phi \text{ is the `pressure'/`geopotential'} |
384 |
\end{equation*} |
385 |
|
386 |
\begin{equation*} |
387 |
\vec{\Omega}\text{ is the Earth's rotation} |
388 |
\end{equation*} |
389 |
|
390 |
\begin{equation*} |
391 |
b\text{ is the `buoyancy'} |
392 |
\end{equation*} |
393 |
|
394 |
\begin{equation*} |
395 |
\theta \text{ is potential temperature} |
396 |
\end{equation*} |
397 |
|
398 |
\begin{equation*} |
399 |
S\text{ is specific humidity in the atmosphere; salinity in the ocean} |
400 |
\end{equation*} |
401 |
|
402 |
\begin{equation*} |
403 |
\mathcal{F}_{\vec{\mathbf{v}}}\text{ are forcing and dissipation of }\vec{ |
404 |
\mathbf{v}} |
405 |
\end{equation*} |
406 |
|
407 |
\begin{equation*} |
408 |
\mathcal{Q}_{\theta }\mathcal{\ }\text{are forcing and dissipation of }\theta |
409 |
\end{equation*} |
410 |
|
411 |
\begin{equation*} |
412 |
\mathcal{Q}_{S}\mathcal{\ }\text{are forcing and dissipation of }S |
413 |
\end{equation*} |
414 |
|
415 |
The $\mathcal{F}^{\prime }s$ and $\mathcal{Q}^{\prime }s$ are provided by |
416 |
extensive `physics' packages for atmosphere and ocean described in Chapter 6. |
417 |
|
418 |
\subsection{Kinematic Boundary conditions} |
419 |
|
420 |
\subsubsection{vertical} |
421 |
|
422 |
at fixed and moving $r$ surfaces we set (see fig.5): |
423 |
|
424 |
\begin{equation} |
425 |
\dot{r}=0atr=R_{fixed}(x,y)\text{ (ocean bottom, top of the atmosphere)} |
426 |
\label{eq:fixedbc} |
427 |
\end{equation} |
428 |
|
429 |
\begin{equation} |
430 |
\dot{r}=\frac{Dr}{Dt}atr=R_{moving}\text{ \ |
431 |
(oceansurface,bottomoftheatmosphere)} \label{eq:movingbc} |
432 |
\end{equation} |
433 |
|
434 |
Here |
435 |
|
436 |
\begin{equation*} |
437 |
R_{moving}=R_{o}+\eta |
438 |
\end{equation*} |
439 |
where $R_{o}(x,y)$ is the `$r-$value' (height or pressure, depending on |
440 |
whether we are in the atmosphere or ocean) of the `moving surface' in the |
441 |
resting fluid and $\eta $ is the departure from $R_{o}(x,y)$ in the presence |
442 |
of motion. |
443 |
|
444 |
\subsubsection{horizontal} |
445 |
|
446 |
\begin{equation} |
447 |
\vec{\mathbf{v}}\cdot \vec{\mathbf{n}}=0 \label{eq:noflow} |
448 |
\end{equation} |
449 |
where $\vec{\mathbf{n}}$ is the normal to a solid boundary. |
450 |
|
451 |
\subsection{Atmosphere} |
452 |
|
453 |
In the atmosphere, see fig.5, we interpret: |
454 |
|
455 |
\begin{equation} |
456 |
r=p\text{ is the pressure} \label{eq:atmos-r} |
457 |
\end{equation} |
458 |
|
459 |
\begin{equation} |
460 |
\dot{r}=\frac{Dp}{Dt}=\omega \text{ is the vertical velocity in }p\text{ |
461 |
coordinates} \label{eq:atmos-omega} |
462 |
\end{equation} |
463 |
|
464 |
\begin{equation} |
465 |
\phi =g\,z\text{ is the geopotential height} \label{eq:atmos-phi} |
466 |
\end{equation} |
467 |
|
468 |
\begin{equation} |
469 |
b=\frac{\partial \Pi }{\partial p}\theta \text{ is the buoyancy} |
470 |
\label{eq:atmos-b} |
471 |
\end{equation} |
472 |
|
473 |
\begin{equation} |
474 |
\theta =T(\frac{p_{c}}{p})^{\kappa }\text{ is potential temperature} |
475 |
\label{eq:atmos-theta} |
476 |
\end{equation} |
477 |
|
478 |
\begin{equation} |
479 |
S=q,\text{ is the specific humidity} \label{eq:atmos-s} |
480 |
\end{equation} |
481 |
where |
482 |
|
483 |
\begin{equation*} |
484 |
T\text{ is absolute temperature} |
485 |
\end{equation*} |
486 |
\begin{equation*} |
487 |
p\text{ is the pressure} |
488 |
\end{equation*} |
489 |
\begin{eqnarray*} |
490 |
&&z\text{ is the height of the pressure surface} \\ |
491 |
&&g\text{ is the acceleration due to gravity} |
492 |
\end{eqnarray*} |
493 |
|
494 |
In the above the ideal gas law, $p=\rho RT$, has been expressed in terms of |
495 |
the Exner function $\Pi (p)$ given by (see Appendix Atmosphere) |
496 |
\begin{equation} |
497 |
\Pi (p)=c_{p}(\frac{p}{p_{c}})^{\kappa } \label{eq:exner} |
498 |
\end{equation} |
499 |
where $p_{c}$ is a reference pressure and $\kappa =R/c_{p}$ with $R$ the gas |
500 |
constant and $c_{p}$ the specific heat of air at constant pressure. |
501 |
|
502 |
At the top of the atmosphere (which is `fixed' in our $r$ coordinate): |
503 |
|
504 |
\begin{equation*} |
505 |
R_{fixed}=p_{top}=0 |
506 |
\end{equation*} |
507 |
In a resting atmosphere the elevation of the mountains at the bottom is |
508 |
given by |
509 |
\begin{equation*} |
510 |
R_{moving}=R_{o}(x,y)=p_{o}(x,y) |
511 |
\end{equation*} |
512 |
i.e. the (hydrostatic) pressure at the top of the mountains in a resting |
513 |
atmosphere. |
514 |
|
515 |
The boundary conditions at top and bottom are given by: |
516 |
|
517 |
\begin{eqnarray} |
518 |
&&\omega =0~\text{at }r=R_{fixed} \text{ (top of the atmosphere)} |
519 |
\label{eq:fixed-bc-atmos} \\ |
520 |
\omega &=&\frac{Dp_{s}}{Dt}\text{; at }r=R_{moving}\text{ (bottom of the |
521 |
atmosphere)} \label{eq:moving-bc-atmos} |
522 |
\end{eqnarray} |
523 |
|
524 |
Then the (hydrostatic form of) eq(\ref{eq:continuous}) yields a consistent |
525 |
set of atmospheric equations which, for convenience, are written out in $p$ |
526 |
coordinates in Appendix Atmosphere - see eqs(\ref{eq:atmos-prime}). |
527 |
|
528 |
\subsection{Ocean} |
529 |
|
530 |
In the ocean we interpret: |
531 |
\begin{eqnarray} |
532 |
r &=&z\text{ is the height} \label{eq:ocean-z} \\ |
533 |
\dot{r} &=&\frac{Dz}{Dt}=w\text{ is the vertical velocity} |
534 |
\label{eq:ocean-w} \\ |
535 |
\phi &=&\frac{p}{\rho _{c}}\text{ is the pressure} \label{eq:ocean-p} \\ |
536 |
b(\theta ,S,r) &=&\frac{g}{\rho _{c}}\left( \rho (\theta ,S,r)-\rho |
537 |
_{c}\right) \text{ is the buoyancy} \label{eq:ocean-b} |
538 |
\end{eqnarray} |
539 |
where $\rho _{c}$ is a fixed reference density of water and $g$ is the |
540 |
acceleration due to gravity.\noindent |
541 |
|
542 |
In the above |
543 |
|
544 |
At the bottom of the ocean: $R_{fixed}(x,y)=-H(x,y)$. |
545 |
|
546 |
The surface of the ocean is given by: $R_{moving}=\eta $ |
547 |
|
548 |
The position of the resting free surface of the ocean is given by $ |
549 |
R_{o}=Z_{o}=0$. |
550 |
|
551 |
Boundary conditions are: |
552 |
|
553 |
\begin{eqnarray} |
554 |
w &=&0~\text{at }r=R_{fixed}\text{ (ocean bottom)} \label{eq:fixed-bc-ocean} |
555 |
\\ |
556 |
w &=&\frac{D\eta }{Dt}\text{ at }r=R_{moving}=\eta \text{ (ocean surface) |
557 |
\label{eq:moving-bc-ocean}} |
558 |
\end{eqnarray} |
559 |
where $\eta $ is the elevation of the free surface. |
560 |
|
561 |
Then eq(\ref{eq:continuous}) yields a consistent set of oceanic equations |
562 |
which, for convenience, are written out in $z$ coordinates in Appendix Ocean |
563 |
- see eqs(\ref{eq:ocean-mom}) to (\ref{eq:ocean-salt}). |
564 |
|
565 |
\subsection{Hydrostatic, Quasi-hydrostatic, Quasi-nonhydrostatic and |
566 |
Non-hydrostatic forms} |
567 |
|
568 |
Let us separate $\phi $ in to surface, hydrostatic and non-hydrostatic terms: |
569 |
|
570 |
\begin{equation} |
571 |
\phi (x,y,r)=\phi _{s}(x,y)+\phi _{hyd}(x,y,r)+\phi _{nh}(x,y,r) |
572 |
\label{eq:phi-split} |
573 |
\end{equation} |
574 |
and write eq(\ref{incompressible}a,b) in the form: |
575 |
|
576 |
\begin{equation} |
577 |
\frac{\partial \vec{\mathbf{v}_{h}}}{\partial t}+\mathbf{\nabla }_{h}\phi |
578 |
_{s}+\mathbf{\nabla }_{h}\phi _{hyd}+\epsilon _{nh}\mathbf{\nabla }_{h}\phi |
579 |
_{nh}=\vec{\mathbf{G}}_{\vec{v}_{h}} \label{eq:mom-h} |
580 |
\end{equation} |
581 |
|
582 |
\begin{equation} |
583 |
\frac{\partial \phi _{hyd}}{\partial r}=-b \label{eq:hydrostatic} |
584 |
\end{equation} |
585 |
|
586 |
\begin{equation} |
587 |
\epsilon _{nh}\frac{\partial \dot{r}}{\partial t}+\frac{\partial \phi _{nh}}{ |
588 |
\partial r}=G_{\dot{r}} \label{eq:mom-w} |
589 |
\end{equation} |
590 |
Here $\epsilon _{nh}$ is a non-hydrostatic parameter. |
591 |
|
592 |
The $\left( \vec{\mathbf{G}}_{\vec{v}},G_{\dot{r}}\right) $ in eq(\ref |
593 |
{eq:mom-h}) and (\ref{eq:mom-w}) represent advective, metric and Coriolis |
594 |
terms in the momentum equations. In spherical coordinates they take the form |
595 |
\footnote{ |
596 |
In the hydrostatic primitive equations (\textbf{HPE}) all underlined terms |
597 |
in (\ref{eq:gu-speherical}), (\ref{eq:gv-spherical}) and (\ref |
598 |
{eq:gw-spherical}) are omitted; the singly-underlined terms are included in |
599 |
the quasi-hydrostatic model (\textbf{QH}). The fully non-hydrostatic model ( |
600 |
\textbf{NH}) includes all terms.} - see Marshall et al 1997a for a full |
601 |
discussion: |
602 |
|
603 |
\begin{equation} |
604 |
\left. |
605 |
\begin{tabular}{l} |
606 |
$G_{u}=-\vec{\mathbf{v}}.\nabla u$ \\ |
607 |
$-\left\{ \underline{\frac{u\dot{r}}{{r}}}-\frac{uv\tan lat}{{r}}\right\} $ |
608 |
\\ |
609 |
$-\left\{ -2\Omega v\sin lat+\underline{2\Omega \dot{r}\cos lat}\right\} $ |
610 |
\\ |
611 |
$+\mathcal{F}_{u}$ |
612 |
\end{tabular} |
613 |
\ \right\} \left\{ |
614 |
\begin{tabular}{l} |
615 |
\textit{advection} \\ |
616 |
\textit{metric} \\ |
617 |
\textit{Coriolis} \\ |
618 |
\textit{\ Forcing/Dissipation} |
619 |
\end{tabular} |
620 |
\ \right. \qquad \label{eq:gu-speherical} |
621 |
\end{equation} |
622 |
|
623 |
\begin{equation} |
624 |
\left. |
625 |
\begin{tabular}{l} |
626 |
$G_{v}=-\vec{\mathbf{v}}.\nabla v$ \\ |
627 |
$-\left\{ \underline{\frac{v\dot{r}}{{r}}}-\frac{u^{2}\tan lat}{{r}}\right\} |
628 |
$ \\ |
629 |
$-\left\{ -2\Omega u\sin lat\right\} $ \\ |
630 |
$+\mathcal{F}_{v}$ |
631 |
\end{tabular} |
632 |
\ \right\} \left\{ |
633 |
\begin{tabular}{l} |
634 |
\textit{advection} \\ |
635 |
\textit{metric} \\ |
636 |
\textit{Coriolis} \\ |
637 |
\textit{\ Forcing/Dissipation} |
638 |
\end{tabular} |
639 |
\ \right. \qquad \label{eq:gv-spherical} |
640 |
\end{equation} |
641 |
\qquad \qquad \qquad \qquad \qquad |
642 |
|
643 |
\begin{equation} |
644 |
\left. |
645 |
\begin{tabular}{l} |
646 |
$G_{\dot{r}}=-\underline{\underline{\vec{\mathbf{v}}.\nabla \dot{r}}}$ \\ |
647 |
$+\left\{ \underline{\frac{u^{_{^{2}}}+v^{2}}{{r}}}\right\} $ \\ |
648 |
${+}\underline{{2\Omega u\cos lat}}$ \\ |
649 |
$\underline{\underline{\mathcal{F}_{\dot{r}}}}$ |
650 |
\end{tabular} |
651 |
\ \right\} \left\{ |
652 |
\begin{tabular}{l} |
653 |
\textit{advection} \\ |
654 |
\textit{metric} \\ |
655 |
\textit{Coriolis} \\ |
656 |
\textit{\ Forcing/Dissipation} |
657 |
\end{tabular} |
658 |
\ \right. \label{eq:gw-spherical} |
659 |
\end{equation} |
660 |
\qquad \qquad \qquad \qquad \qquad |
661 |
|
662 |
In the above `${r}$' is the distance from the center of the earth and `$lat$ |
663 |
' is latitude. |
664 |
|
665 |
Grad and div operators in spherical coordinates are defined in appendix |
666 |
OPERATORS. |
667 |
\marginpar{ |
668 |
Fig.6 Spherical polar coordinate system.} |
669 |
|
670 |
%%CNHbegin |
671 |
%notci%\input{part1/sphere_coord_figure.tex} |
672 |
%%CNHend |
673 |
|
674 |
\subsubsection{Shallow atmosphere approximation} |
675 |
|
676 |
Most models are based on the `hydrostatic primitive equations' (HPE's) in |
677 |
which the vertical momentum equation is reduced to a statement of |
678 |
hydrostatic balance and the `traditional approximation' is made in which the |
679 |
Coriolis force is treated approximately and the shallow atmosphere |
680 |
approximation is made.\ The MITgcm need not make the `traditional |
681 |
approximation'. To be able to support consistent non-hydrostatic forms the |
682 |
shallow atmosphere approximation can be relaxed - when dividing through by $ |
683 |
r $ in, for example, (\ref{eq:gu-speherical}), we do not replace $r$ by $a$, |
684 |
the radius of the earth. |
685 |
|
686 |
\subsubsection{Hydrostatic and quasi-hydrostatic forms} |
687 |
|
688 |
These are discussed at length in Marshall et al (1997a). |
689 |
|
690 |
In the `hydrostatic primitive equations' (\textbf{HPE)} all the underlined |
691 |
terms in Eqs. (\ref{eq:gu-speherical} $\rightarrow $\ \ref{eq:gw-spherical}) |
692 |
are neglected and `${r}$' is replaced by `$a$', the mean radius of the |
693 |
earth. Once the pressure is found at one level - e.g. by inverting a 2-d |
694 |
Elliptic equation for $\phi _{s}$ at $r=R_{moving}$ - the pressure can be |
695 |
computed at all other levels by integration of the hydrostatic relation, eq( |
696 |
\ref{eq:hydrostatic}). |
697 |
|
698 |
In the `quasi-hydrostatic' equations (\textbf{QH)} strict balance between |
699 |
gravity and vertical pressure gradients is not imposed. The $2\Omega u\cos |
700 |
\phi $ Coriolis term are not neglected and are balanced by a non-hydrostatic |
701 |
contribution to the pressure field: only the terms underlined twice in Eqs. ( |
702 |
\ref{eq:gu-speherical}$\rightarrow $\ \ref{eq:gw-spherical}) are set to zero |
703 |
and, simultaneously, the shallow atmosphere approximation is relaxed. In |
704 |
\textbf{QH}\ \textit{all} the metric terms are retained and the full |
705 |
variation of the radial position of a particle monitored. The \textbf{QH}\ |
706 |
vertical momentum equation (\ref{eq:mom-w}) becomes: |
707 |
|
708 |
\begin{equation*} |
709 |
\frac{\partial \phi _{nh}}{\partial r}=2\Omega u\cos lat |
710 |
\end{equation*} |
711 |
making a small correction to the hydrostatic pressure. |
712 |
|
713 |
\textbf{QH} has good energetic credentials - they are the same as for |
714 |
\textbf{HPE}. Importantly, however, it has the same angular momentum |
715 |
principle as the full non-hydrostatic model (\textbf{NH)} - see Marshall |
716 |
et.al., 1997a. As in \textbf{HPE }only a 2-d elliptic problem need be solved. |
717 |
|
718 |
\subsubsection{Non-hydrostatic and quasi-nonhydrostatic forms} |
719 |
|
720 |
The MIT model presently supports a full non-hydrostatic ocean isomorph, but |
721 |
only a quasi-non-hydrostatic atmospheric isomorph. |
722 |
|
723 |
\paragraph{Non-hydrostatic Ocean} |
724 |
|
725 |
In the non-hydrostatic ocean model all terms in equations Eqs.(\ref |
726 |
{eq:gu-speherical} $\rightarrow $\ \ref{eq:gw-spherical}) are retained. A |
727 |
three dimensional elliptic equation must be solved subject to Neumann |
728 |
boundary conditions (see below). It is important to note that use of the |
729 |
full \textbf{NH} does not admit any new `fast' waves in to the system - the |
730 |
incompressible condition eq(\ref{eq:continuous})c has already filtered out |
731 |
acoustic modes. It does, however, ensure that the gravity waves are treated |
732 |
accurately with an exact dispersion relation. The \textbf{NH} set has a |
733 |
complete angular momentum principle and consistent energetics - see White |
734 |
and Bromley, 1995; Marshall et.al.\ 1997a. |
735 |
|
736 |
\paragraph{Quasi-nonhydrostatic Atmosphere} |
737 |
|
738 |
In the non-hydrostatic version of our atmospheric model we approximate $\dot{ |
739 |
r}$ in the vertical momentum eqs(\ref{eq:mom-w}) and (\ref{eq:gv-spherical}) |
740 |
(but only here) by: |
741 |
|
742 |
\begin{equation} |
743 |
\dot{r}=\frac{Dp}{Dt}=\frac{1}{g}\frac{D\phi }{Dt} \label{eq:quasi-nh-w} |
744 |
\end{equation} |
745 |
where $p_{hy}$ is the hydrostatic pressure. |
746 |
|
747 |
\subsubsection{Summary of equation sets supported by model} |
748 |
|
749 |
\paragraph{Atmosphere} |
750 |
|
751 |
Hydrostatic, and quasi-hydrostatic and quasi non-hydrostatic forms of the |
752 |
compressible non-Boussinesq equations in $p-$coordinates are supported. |
753 |
|
754 |
\subparagraph{Hydrostatic and quasi-hydrostatic} |
755 |
|
756 |
The hydrostatic set is written out in $p-$coordinates in appendix Atmosphere |
757 |
- see eq(\ref{eq:atmos-prime}). |
758 |
|
759 |
\subparagraph{Quasi-nonhydrostatic} |
760 |
|
761 |
A quasi-nonhydrostatic form is also supported. |
762 |
|
763 |
\paragraph{Ocean} |
764 |
|
765 |
\subparagraph{Hydrostatic and quasi-hydrostatic} |
766 |
|
767 |
Hydrostatic, and quasi-hydrostatic forms of the incompressible Boussinesq |
768 |
equations in $z-$coordinates are supported. |
769 |
|
770 |
\subparagraph{Non-hydrostatic} |
771 |
|
772 |
Non-hydrostatic forms of the incompressible Boussinesq equations in $z-$ |
773 |
coordinates are supported - see eqs(\ref{eq:ocean-mom}) to (\ref |
774 |
{eq:ocean-salt}). |
775 |
|
776 |
\subsection{Solution strategy} |
777 |
|
778 |
The method of solution employed in the \textbf{HPE}, \textbf{QH} and \textbf{ |
779 |
NH} models is summarized in Fig.7. |
780 |
\marginpar{ |
781 |
Fig.7 Solution strategy} Under all dynamics, a 2-d elliptic equation is |
782 |
first solved to find the surface pressure and the hydrostatic pressure at |
783 |
any level computed from the weight of fluid above. Under \textbf{HPE} and |
784 |
\textbf{QH} dynamics, the horizontal momentum equations are then stepped |
785 |
forward and $\dot{r}$ found from continuity. Under \textbf{NH} dynamics a |
786 |
3-d elliptic equation must be solved for the non-hydrostatic pressure before |
787 |
stepping forward the horizontal momentum equations; $\dot{r}$ is found by |
788 |
stepping forward the vertical momentum equation. |
789 |
|
790 |
%%CNHbegin |
791 |
%notci%\input{part1/solution_strategy_figure.tex} |
792 |
%%CNHend |
793 |
|
794 |
There is no penalty in implementing \textbf{QH} over \textbf{HPE} except, of |
795 |
course, some complication that goes with the inclusion of $\cos \phi \ $ |
796 |
Coriolis terms and the relaxation of the shallow atmosphere approximation. |
797 |
But this leads to negligible increase in computation. In \textbf{NH}, in |
798 |
contrast, one additional elliptic equation - a three-dimensional one - must |
799 |
be inverted for $p_{nh}$. However the `overhead' of the \textbf{NH} model is |
800 |
essentially negligible in the hydrostatic limit (see detailed discussion in |
801 |
Marshall et al, 1997) resulting in a non-hydrostatic algorithm that, in the |
802 |
hydrostatic limit, is as computationally economic as the \textbf{HPEs}. |
803 |
|
804 |
\subsection{Finding the pressure field} |
805 |
|
806 |
Unlike the prognostic variables $u$, $v$, $w$, $\theta $ and $S$, the |
807 |
pressure field must be obtained diagnostically. We proceed, as before, by |
808 |
dividing the total (pressure/geo) potential in to three parts, a surface |
809 |
part, $\phi _{s}(x,y)$, a hydrostatic part $\phi _{hyd}(x,y,r)$ and a |
810 |
non-hydrostatic part $\phi _{nh}(x,y,r)$, as in (\ref{eq:phi-split}), and |
811 |
writing the momentum equation as in (\ref{eq:mom-h}). |
812 |
|
813 |
\subsubsection{Hydrostatic pressure} |
814 |
|
815 |
Hydrostatic pressure is obtained by integrating (\ref{eq:hydrostatic}) |
816 |
vertically from $r=R_{o}$ where $\phi _{hyd}(r=R_{o})=0$, to yield: |
817 |
|
818 |
\begin{equation*} |
819 |
\int_{r}^{R_{o}}\frac{\partial \phi _{hyd}}{\partial r}dr=\left[ \phi _{hyd} |
820 |
\right] _{r}^{R_{o}}=\int_{r}^{R_{o}}-bdr |
821 |
\end{equation*} |
822 |
and so |
823 |
|
824 |
\begin{equation} |
825 |
\phi _{hyd}(x,y,r)=\int_{r}^{R_{o}}bdr \label{eq:hydro-phi} |
826 |
\end{equation} |
827 |
|
828 |
The model can be easily modified to accommodate a loading term (e.g |
829 |
atmospheric pressure pushing down on the ocean's surface) by setting: |
830 |
|
831 |
\begin{equation} |
832 |
\phi _{hyd}(r=R_{o})=loading \label{eq:loading} |
833 |
\end{equation} |
834 |
|
835 |
\subsubsection{Surface pressure} |
836 |
|
837 |
The surface pressure equation can be obtained by integrating continuity, ( |
838 |
\ref{eq:continuous})c, vertically from $r=R_{fixed}$ to $r=R_{moving}$ |
839 |
|
840 |
\begin{equation*} |
841 |
\int_{R_{fixed}}^{R_{moving}}\left( \mathbf{\nabla }_{h}\cdot \vec{\mathbf{v} |
842 |
}_{h}+\partial _{r}\dot{r}\right) dr=0 |
843 |
\end{equation*} |
844 |
|
845 |
Thus: |
846 |
|
847 |
\begin{equation*} |
848 |
\frac{\partial \eta }{\partial t}+\vec{\mathbf{v}}.\nabla \eta |
849 |
+\int_{R_{fixed}}^{R_{moving}}\mathbf{\nabla }_{h}\cdot \vec{\mathbf{v}} |
850 |
_{h}dr=0 |
851 |
\end{equation*} |
852 |
where $\eta =R_{moving}-R_{o}$ is the free-surface $r$-anomaly in units of $ |
853 |
r $. The above can be rearranged to yield, using Leibnitz's theorem: |
854 |
|
855 |
\begin{equation} |
856 |
\frac{\partial \eta }{\partial t}+\mathbf{\nabla }_{h}\cdot |
857 |
\int_{R_{fixed}}^{R_{moving}}\vec{\mathbf{v}}_{h}dr=\text{source} |
858 |
\label{eq:free-surface} |
859 |
\end{equation} |
860 |
where we have incorporated a source term. |
861 |
|
862 |
Whether $\phi $ is pressure (ocean model, $p/\rho _{c}$) or geopotential |
863 |
(atmospheric model), in (\ref{mtm-split}), the horizontal gradient term can |
864 |
be written |
865 |
\begin{equation} |
866 |
\mathbf{\nabla }_{h}\phi _{s}=\mathbf{\nabla }_{h}\left( b_{s}\eta \right) |
867 |
\label{eq:phi-surf} |
868 |
\end{equation} |
869 |
where $b_{s}$ is the buoyancy at the surface. |
870 |
|
871 |
In the hydrostatic limit ($\epsilon _{nh}=0$), Eqs(\ref{eq:mom-h}), (\ref |
872 |
{eq:free-surface}) and (\ref{eq:phi-surf}) can be solved by inverting a 2-d |
873 |
elliptic equation for $\phi _{s}$ as described in Chapter 2. Both `free |
874 |
surface' and `rigid lid' approaches are available. |
875 |
|
876 |
\subsubsection{Non-hydrostatic pressure} |
877 |
|
878 |
Taking the horizontal divergence of (\ref{hor-mtm}) and adding $\frac{ |
879 |
\partial }{\partial r}$ of (\ref{vertmtm}), invoking the continuity equation |
880 |
(\ref{incompressible}), we deduce that: |
881 |
|
882 |
\begin{equation} |
883 |
\nabla _{3}^{2}\phi _{nh}=\nabla .\vec{\mathbf{G}}_{\vec{v}}-\left( \mathbf{ |
884 |
\nabla }_{h}^{2}\phi _{s}+\mathbf{\nabla }^{2}\phi _{hyd}\right) =\nabla . |
885 |
\vec{\mathbf{F}} \label{eq:3d-invert} |
886 |
\end{equation} |
887 |
|
888 |
For a given rhs this 3-d elliptic equation must be inverted for $\phi _{nh}$ |
889 |
subject to appropriate choice of boundary conditions. This method is usually |
890 |
called \textit{The Pressure Method} [Harlow and Welch, 1965; Williams, 1969; |
891 |
Potter, 1976]. In the hydrostatic primitive equations case (\textbf{HPE}), |
892 |
the 3-d problem does not need to be solved. |
893 |
|
894 |
\paragraph{Boundary Conditions} |
895 |
|
896 |
We apply the condition of no normal flow through all solid boundaries - the |
897 |
coasts (in the ocean) and the bottom: |
898 |
|
899 |
\begin{equation} |
900 |
\vec{\mathbf{v}}.\widehat{n}=0 \label{nonormalflow} |
901 |
\end{equation} |
902 |
where $\widehat{n}$ is a vector of unit length normal to the boundary. The |
903 |
kinematic condition (\ref{nonormalflow}) is also applied to the vertical |
904 |
velocity at $r=R_{moving}$. No-slip $\left( v_{T}=0\right) \ $or slip $ |
905 |
\left( \partial v_{T}/\partial n=0\right) \ $conditions are employed on the |
906 |
tangential component of velocity, $v_{T}$, at all solid boundaries, |
907 |
depending on the form chosen for the dissipative terms in the momentum |
908 |
equations - see below. |
909 |
|
910 |
Eq.(\ref{nonormalflow}) implies, making use of (\ref{mtm-split}), that: |
911 |
|
912 |
\begin{equation} |
913 |
\widehat{n}.\nabla \phi _{nh}=\widehat{n}.\vec{\mathbf{F}} |
914 |
\label{eq:inhom-neumann-nh} |
915 |
\end{equation} |
916 |
where |
917 |
|
918 |
\begin{equation*} |
919 |
\vec{\mathbf{F}}=\vec{\mathbf{G}}_{\vec{v}}-\left( \mathbf{\nabla }_{h}\phi |
920 |
_{s}+\mathbf{\nabla }\phi _{hyd}\right) |
921 |
\end{equation*} |
922 |
presenting inhomogeneous Neumann boundary conditions to the Elliptic problem |
923 |
(\ref{eq:3d-invert}). As shown, for example, by Williams (1969), one can |
924 |
exploit classical 3D potential theory and, by introducing an appropriately |
925 |
chosen $\delta $-function sheet of `source-charge', replace the |
926 |
inhomogeneous boundary condition on pressure by a homogeneous one. The |
927 |
source term $rhs$ in (\ref{eq:3d-invert}) is the divergence of the vector $ |
928 |
\vec{\mathbf{F}}.$ By simultaneously setting $ |
929 |
\begin{array}{l} |
930 |
\widehat{n}.\vec{\mathbf{F}} |
931 |
\end{array} |
932 |
=0$\ and $\widehat{n}.\nabla \phi _{nh}=0\ $on the boundary the following |
933 |
self-consistent but simpler homogenized Elliptic problem is obtained: |
934 |
|
935 |
\begin{equation*} |
936 |
\nabla ^{2}\phi _{nh}=\nabla .\widetilde{\vec{\mathbf{F}}}\qquad |
937 |
\end{equation*} |
938 |
where $\widetilde{\vec{\mathbf{F}}}$ is a modified $\vec{\mathbf{F}}$ such |
939 |
that $\widetilde{\vec{\mathbf{F}}}.\widehat{n}=0$. As is implied by (\ref |
940 |
{eq:inhom-neumann-nh}) the modified boundary condition becomes: |
941 |
|
942 |
\begin{equation} |
943 |
\widehat{n}.\nabla \phi _{nh}=0 \label{eq:hom-neumann-nh} |
944 |
\end{equation} |
945 |
|
946 |
If the flow is `close' to hydrostatic balance then the 3-d inversion |
947 |
converges rapidly because $\phi _{nh}\ $is then only a small correction to |
948 |
the hydrostatic pressure field (see the discussion in Marshall et al, a,b). |
949 |
|
950 |
The solution $\phi _{nh}\ $to (\ref{eq:3d-invert}) and (\ref{homneuman}) |
951 |
does not vanish at $r=R_{moving}$, and so refines the pressure there. |
952 |
|
953 |
\subsection{Forcing/dissipation} |
954 |
|
955 |
\subsubsection{Forcing} |
956 |
|
957 |
The forcing terms $\mathcal{F}$ on the rhs of the equations are provided by |
958 |
`physics packages' described in detail in chapter ??. |
959 |
|
960 |
\subsubsection{Dissipation} |
961 |
|
962 |
\paragraph{Momentum} |
963 |
|
964 |
Many forms of momentum dissipation are available in the model. Laplacian and |
965 |
biharmonic frictions are commonly used: |
966 |
|
967 |
\begin{equation} |
968 |
D_{V}=A_{h}\nabla _{h}^{2}v+A_{v}\frac{\partial ^{2}v}{\partial z^{2}} |
969 |
+A_{4}\nabla _{h}^{4}v \label{eq:dissipation} |
970 |
\end{equation} |
971 |
where $A_{h}$ and $A_{v}\ $are (constant) horizontal and vertical viscosity |
972 |
coefficients and $A_{4}\ $is the horizontal coefficient for biharmonic |
973 |
friction. These coefficients are the same for all velocity components. |
974 |
|
975 |
\paragraph{Tracers} |
976 |
|
977 |
The mixing terms for the temperature and salinity equations have a similar |
978 |
form to that of momentum except that the diffusion tensor can be |
979 |
non-diagonal and have varying coefficients. $\qquad $ |
980 |
\begin{equation} |
981 |
D_{T,S}=\nabla .[\underline{\underline{K}}\nabla (T,S)]+K_{4}\nabla |
982 |
_{h}^{4}(T,S) \label{eq:diffusion} |
983 |
\end{equation} |
984 |
where $\underline{\underline{K}}\ $is the diffusion tensor and the $K_{4}\ $ |
985 |
horizontal coefficient for biharmonic diffusion. In the simplest case where |
986 |
the subgrid-scale fluxes of heat and salt are parameterized with constant |
987 |
horizontal and vertical diffusion coefficients, $\underline{\underline{K}}$, |
988 |
reduces to a diagonal matrix with constant coefficients: |
989 |
|
990 |
\begin{equation} |
991 |
\qquad \qquad \qquad \qquad K=\left( |
992 |
\begin{array}{ccc} |
993 |
K_{h} & 0 & 0 \\ |
994 |
0 & K_{h} & 0 \\ |
995 |
0 & 0 & K_{v} |
996 |
\end{array} |
997 |
\right) \qquad \qquad \qquad \label{eq:diagonal-diffusion-tensor} |
998 |
\end{equation} |
999 |
where $K_{h}\ $and $K_{v}\ $are the horizontal and vertical diffusion |
1000 |
coefficients. These coefficients are the same for all tracers (temperature, |
1001 |
salinity ... ). |
1002 |
|
1003 |
\subsection{Vector invariant form} |
1004 |
|
1005 |
For some purposes it is advantageous to write momentum advection in eq(\ref |
1006 |
{hor-mtm}) and (\ref{vertmtm}) in the (so-called) `vector invariant' form: |
1007 |
|
1008 |
\begin{equation} |
1009 |
\frac{D\vec{\mathbf{v}}}{Dt}=\frac{\partial \vec{\mathbf{v}}}{\partial t} |
1010 |
+\left( \nabla \times \vec{\mathbf{v}}\right) \times \vec{\mathbf{v}}+\nabla |
1011 |
\left[ \frac{1}{2}(\vec{\mathbf{v}}\cdot \vec{\mathbf{v}})\right] |
1012 |
\label{eq:vi-identity} |
1013 |
\end{equation} |
1014 |
This permits alternative numerical treatments of the non-linear terms based |
1015 |
on their representation as a vorticity flux. Because gradients of coordinate |
1016 |
vectors no longer appear on the rhs of (\ref{eq:vi-identity}), explicit |
1017 |
representation of the metric terms in (\ref{eq:gu-speherical}), (\ref |
1018 |
{eq:gv-spherical}) and (\ref{eq:gw-spherical}), can be avoided: information |
1019 |
about the geometry is contained in the areas and lengths of the volumes used |
1020 |
to discretize the model. |
1021 |
|
1022 |
\subsection{Adjoint} |
1023 |
|
1024 |
Tangent linear and adjoint counterparts of the forward model and described |
1025 |
in Chapter 5. |
1026 |
|
1027 |
% $Header: /u/gcmpack/mitgcmdoc/part1/manual.tex,v 1.3 2001/10/10 16:48:45 cnh Exp $ |
1028 |
% $Name: $ |
1029 |
|
1030 |
\section{Appendix ATMOSPHERE} |
1031 |
|
1032 |
\subsection{Hydrostatic Primitive Equations for the Atmosphere in pressure |
1033 |
coordinates} |
1034 |
|
1035 |
\label{sect-hpe-p} |
1036 |
|
1037 |
The hydrostatic primitive equations (HPEs) in p-coordinates are: |
1038 |
\begin{eqnarray} |
1039 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1040 |
_{h}+\mathbf{\nabla }_{p}\phi &=&\vec{\mathbf{\mathcal{F}}} |
1041 |
\label{eq:atmos-mom} \\ |
1042 |
\frac{\partial \phi }{\partial p}+\alpha &=&0 \label{eq-p-hydro-start} \\ |
1043 |
\mathbf{\nabla }_{p}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial \omega }{ |
1044 |
\partial p} &=&0 \label{eq:atmos-cont} \\ |
1045 |
p\alpha &=&RT \label{eq:atmos-eos} \\ |
1046 |
c_{v}\frac{DT}{Dt}+p\frac{D\alpha }{Dt} &=&\mathcal{Q} \label{eq:atmos-heat} |
1047 |
\end{eqnarray} |
1048 |
where $\vec{\mathbf{v}}_{h}=(u,v,0)$ is the `horizontal' (on pressure |
1049 |
surfaces) component of velocity,$\frac{D}{Dt}=\vec{\mathbf{v}}_{h}\cdot |
1050 |
\mathbf{\nabla }_{p}+\omega \frac{\partial }{\partial p}$ is the total |
1051 |
derivative, $f=2\Omega \sin lat$ is the Coriolis parameter, $\phi =gz$ is |
1052 |
the geopotential, $\alpha =1/\rho $ is the specific volume, $\omega =\frac{Dp |
1053 |
}{Dt}$ is the vertical velocity in the $p-$coordinate. Equation(\ref |
1054 |
{eq:atmos-heat}) is the first law of thermodynamics where internal energy $ |
1055 |
e=c_{v}T$, $T$ is temperature, $Q$ is the rate of heating per unit mass and $ |
1056 |
p\frac{D\alpha }{Dt}$ is the work done by the fluid in compressing. |
1057 |
|
1058 |
It is convenient to cast the heat equation in terms of potential temperature |
1059 |
$\theta $ so that it looks more like a generic conservation law. |
1060 |
Differentiating (\ref{eq:atmos-eos}) we get: |
1061 |
\begin{equation*} |
1062 |
p\frac{D\alpha }{Dt}+\alpha \frac{Dp}{Dt}=R\frac{DT}{Dt} |
1063 |
\end{equation*} |
1064 |
which, when added to the heat equation (\ref{eq:atmos-heat}) and using $ |
1065 |
c_{p}=c_{v}+R$, gives: |
1066 |
\begin{equation} |
1067 |
c_{p}\frac{DT}{Dt}-\alpha \frac{Dp}{Dt}=\mathcal{Q} |
1068 |
\label{eq-p-heat-interim} |
1069 |
\end{equation} |
1070 |
Potential temperature is defined: |
1071 |
\begin{equation} |
1072 |
\theta =T(\frac{p_{c}}{p})^{\kappa } \label{eq:potential-temp} |
1073 |
\end{equation} |
1074 |
where $p_{c}$ is a reference pressure and $\kappa =R/c_{p}$. For convenience |
1075 |
we will make use of the Exner function $\Pi (p)$ which defined by: |
1076 |
\begin{equation} |
1077 |
\Pi (p)=c_{p}(\frac{p}{p_{c}})^{\kappa } \label{Exner} |
1078 |
\end{equation} |
1079 |
The following relations will be useful and are easily expressed in terms of |
1080 |
the Exner function: |
1081 |
\begin{equation*} |
1082 |
c_{p}T=\Pi \theta \;\;;\;\;\frac{\partial \Pi }{\partial p}=\frac{\kappa \Pi |
1083 |
}{p}\;\;;\;\;\alpha =\frac{\kappa \Pi \theta }{p}=\frac{\partial \ \Pi }{ |
1084 |
\partial p}\theta \;\;;\;\;\frac{D\Pi }{Dt}=\frac{\partial \Pi }{\partial p} |
1085 |
\frac{Dp}{Dt} |
1086 |
\end{equation*} |
1087 |
where $b=\frac{\partial \ \Pi }{\partial p}\theta $ is the buoyancy. |
1088 |
|
1089 |
The heat equation is obtained by noting that |
1090 |
\begin{equation*} |
1091 |
c_{p}\frac{DT}{Dt}=\frac{D(\Pi \theta )}{Dt}=\Pi \frac{D\theta }{Dt}+\theta |
1092 |
\frac{D\Pi }{Dt}=\Pi \frac{D\theta }{Dt}+\alpha \frac{Dp}{Dt} |
1093 |
\end{equation*} |
1094 |
and on substituting into (\ref{eq-p-heat-interim}) gives: |
1095 |
\begin{equation} |
1096 |
\Pi \frac{D\theta }{Dt}=\mathcal{Q} |
1097 |
\label{eq:potential-temperature-equation} |
1098 |
\end{equation} |
1099 |
which is in conservative form. |
1100 |
|
1101 |
For convenience in the model we prefer to step forward (\ref |
1102 |
{eq:potential-temperature-equation}) rather than (\ref{eq:atmos-heat}). |
1103 |
|
1104 |
\subsubsection{Boundary conditions} |
1105 |
|
1106 |
The upper and lower boundary conditions are : |
1107 |
\begin{eqnarray} |
1108 |
\mbox{at the top:}\;\;p=0 &&\text{, }\omega =\frac{Dp}{Dt}=0 \\ |
1109 |
\mbox{at the surface:}\;\;p=p_{s} &&\text{, }\phi =\phi _{topo}=g~Z_{topo} |
1110 |
\label{eq:boundary-condition-atmosphere} |
1111 |
\end{eqnarray} |
1112 |
In $p$-coordinates, the upper boundary acts like a solid boundary ($\omega |
1113 |
=0 $); in $z$-coordinates and the lower boundary is analogous to a free |
1114 |
surface ($\phi $ is imposed and $\omega \neq 0$). |
1115 |
|
1116 |
\subsubsection{Splitting the geo-potential} |
1117 |
|
1118 |
For the purposes of initialization and reducing round-off errors, the model |
1119 |
deals with perturbations from reference (or ``standard'') profiles. For |
1120 |
example, the hydrostatic geopotential associated with the resting atmosphere |
1121 |
is not dynamically relevant and can therefore be subtracted from the |
1122 |
equations. The equations written in terms of perturbations are obtained by |
1123 |
substituting the following definitions into the previous model equations: |
1124 |
\begin{eqnarray} |
1125 |
\theta &=&\theta _{o}+\theta ^{\prime } \label{eq:atmos-ref-prof-theta} \\ |
1126 |
\alpha &=&\alpha _{o}+\alpha ^{\prime } \label{eq:atmos-ref-prof-alpha} \\ |
1127 |
\phi &=&\phi _{o}+\phi ^{\prime } \label{eq:atmos-ref-prof-phi} |
1128 |
\end{eqnarray} |
1129 |
The reference state (indicated by subscript ``0'') corresponds to |
1130 |
horizontally homogeneous atmosphere at rest ($\theta _{o},\alpha _{o},\phi |
1131 |
_{o}$) with surface pressure $p_{o}(x,y)$ that satisfies $\phi |
1132 |
_{o}(p_{o})=g~Z_{topo}$, defined: |
1133 |
\begin{eqnarray*} |
1134 |
\theta _{o}(p) &=&f^{n}(p) \\ |
1135 |
\alpha _{o}(p) &=&\Pi _{p}\theta _{o} \\ |
1136 |
\phi _{o}(p) &=&\phi _{topo}-\int_{p_{0}}^{p}\alpha _{o}dp |
1137 |
\end{eqnarray*} |
1138 |
%\begin{eqnarray*} |
1139 |
%\phi'_\alpha & = & \int^p_{p_o} (\alpha_o -\alpha) dp \\ |
1140 |
%\phi'_s(x,y,t) & = & \int_{p_o}^{p_s} \alpha dp |
1141 |
%\end{eqnarray*} |
1142 |
|
1143 |
The final form of the HPE's in p coordinates is then: |
1144 |
\begin{eqnarray} |
1145 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1146 |
_{h}+\mathbf{\nabla }_{p}\phi ^{\prime } &=&\vec{\mathbf{\mathcal{F}}} \\ |
1147 |
\frac{\partial \phi ^{\prime }}{\partial p}+\alpha ^{\prime } &=&0 \\ |
1148 |
\mathbf{\nabla }_{p}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial \omega }{ |
1149 |
\partial p} &=&0 \\ |
1150 |
\frac{\partial \Pi }{\partial p}\theta ^{\prime } &=&\alpha ^{\prime } \\ |
1151 |
\frac{D\theta }{Dt} &=&\frac{\mathcal{Q}}{\Pi } \label{eq:atmos-prime} |
1152 |
\end{eqnarray} |
1153 |
|
1154 |
% $Header: /u/gcmpack/mitgcmdoc/part1/manual.tex,v 1.3 2001/10/10 16:48:45 cnh Exp $ |
1155 |
% $Name: $ |
1156 |
|
1157 |
\section{Appendix OCEAN} |
1158 |
|
1159 |
\subsection{Equations of motion for the ocean} |
1160 |
|
1161 |
We review here the method by which the standard (Boussinesq, incompressible) |
1162 |
HPE's for the ocean written in z-coordinates are obtained. The |
1163 |
non-Boussinesq equations for oceanic motion are: |
1164 |
\begin{eqnarray} |
1165 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1166 |
_{h}+\frac{1}{\rho }\mathbf{\nabla }_{z}p &=&\vec{\mathbf{\mathcal{F}}} \\ |
1167 |
\epsilon _{nh}\frac{Dw}{Dt}+g+\frac{1}{\rho }\frac{\partial p}{\partial z} |
1168 |
&=&\epsilon _{nh}\mathcal{F}_{w} \\ |
1169 |
\frac{1}{\rho }\frac{D\rho }{Dt}+\mathbf{\nabla }_{z}\cdot \vec{\mathbf{v}} |
1170 |
_{h}+\frac{\partial w}{\partial z} &=&0 \\ |
1171 |
\rho &=&\rho (\theta ,S,p) \\ |
1172 |
\frac{D\theta }{Dt} &=&\mathcal{Q}_{\theta } \\ |
1173 |
\frac{DS}{Dt} &=&\mathcal{Q}_{s} \label{eq:non-boussinesq} |
1174 |
\end{eqnarray} |
1175 |
These equations permit acoustics modes, inertia-gravity waves, |
1176 |
non-hydrostatic motions, a geostrophic (Rossby) mode and a thermo-haline |
1177 |
mode. As written, they cannot be integrated forward consistently - if we |
1178 |
step $\rho $ forward in (\ref{eq-zns-cont}), the answer will not be |
1179 |
consistent with that obtained by stepping (\ref{eq-zns-heat}) and (\ref |
1180 |
{eq-zns-salt}) and then using (\ref{eq-zns-eos}) to yield $\rho $. It is |
1181 |
therefore necessary to manipulate the system as follows. Differentiating the |
1182 |
EOS (equation of state) gives: |
1183 |
|
1184 |
\begin{equation} |
1185 |
\frac{D\rho }{Dt}=\left. \frac{\partial \rho }{\partial \theta }\right| |
1186 |
_{S,p}\frac{D\theta }{Dt}+\left. \frac{\partial \rho }{\partial S}\right| |
1187 |
_{\theta ,p}\frac{DS}{Dt}+\left. \frac{\partial \rho }{\partial p}\right| |
1188 |
_{\theta ,S}\frac{Dp}{Dt} \label{EOSexpansion} |
1189 |
\end{equation} |
1190 |
|
1191 |
Note that $\frac{\partial \rho }{\partial p}=\frac{1}{c_{s}^{2}}$ is the |
1192 |
reciprocal of the sound speed ($c_{s}$) squared. Substituting into \ref |
1193 |
{eq-zns-cont} gives: |
1194 |
\begin{equation} |
1195 |
\frac{1}{\rho c_{s}^{2}}\frac{Dp}{Dt}+\mathbf{\nabla }_{z}\cdot \vec{\mathbf{ |
1196 |
v}}+\partial _{z}w\approx 0 \label{eq-zns-pressure} |
1197 |
\end{equation} |
1198 |
where we have used an approximation sign to indicate that we have assumed |
1199 |
adiabatic motion, dropping the $\frac{D\theta }{Dt}$ and $\frac{DS}{Dt}$. |
1200 |
Replacing \ref{eq-zns-cont} with \ref{eq-zns-pressure} yields a system that |
1201 |
can be explicitly integrated forward: |
1202 |
\begin{eqnarray} |
1203 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1204 |
_{h}+\frac{1}{\rho }\mathbf{\nabla }_{z}p &=&\vec{\mathbf{\mathcal{F}}} |
1205 |
\label{eq-cns-hmom} \\ |
1206 |
\epsilon _{nh}\frac{Dw}{Dt}+g+\frac{1}{\rho }\frac{\partial p}{\partial z} |
1207 |
&=&\epsilon _{nh}\mathcal{F}_{w} \label{eq-cns-hydro} \\ |
1208 |
\frac{1}{\rho c_{s}^{2}}\frac{Dp}{Dt}+\mathbf{\nabla }_{z}\cdot \vec{\mathbf{ |
1209 |
v}}_{h}+\frac{\partial w}{\partial z} &=&0 \label{eq-cns-cont} \\ |
1210 |
\rho &=&\rho (\theta ,S,p) \label{eq-cns-eos} \\ |
1211 |
\frac{D\theta }{Dt} &=&\mathcal{Q}_{\theta } \label{eq-cns-heat} \\ |
1212 |
\frac{DS}{Dt} &=&\mathcal{Q}_{s} \label{eq-cns-salt} |
1213 |
\end{eqnarray} |
1214 |
|
1215 |
\subsubsection{Compressible z-coordinate equations} |
1216 |
|
1217 |
Here we linearize the acoustic modes by replacing $\rho $ with $\rho _{o}(z)$ |
1218 |
wherever it appears in a product (ie. non-linear term) - this is the |
1219 |
`Boussinesq assumption'. The only term that then retains the full variation |
1220 |
in $\rho $ is the gravitational acceleration: |
1221 |
\begin{eqnarray} |
1222 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1223 |
_{h}+\frac{1}{\rho _{o}}\mathbf{\nabla }_{z}p &=&\vec{\mathbf{\mathcal{F}}} |
1224 |
\label{eq-zcb-hmom} \\ |
1225 |
\epsilon _{nh}\frac{Dw}{Dt}+\frac{g\rho }{\rho _{o}}+\frac{1}{\rho _{o}} |
1226 |
\frac{\partial p}{\partial z} &=&\epsilon _{nh}\mathcal{F}_{w} |
1227 |
\label{eq-zcb-hydro} \\ |
1228 |
\frac{1}{\rho _{o}c_{s}^{2}}\frac{Dp}{Dt}+\mathbf{\nabla }_{z}\cdot \vec{ |
1229 |
\mathbf{v}}_{h}+\frac{\partial w}{\partial z} &=&0 \label{eq-zcb-cont} \\ |
1230 |
\rho &=&\rho (\theta ,S,p) \label{eq-zcb-eos} \\ |
1231 |
\frac{D\theta }{Dt} &=&\mathcal{Q}_{\theta } \label{eq-zcb-heat} \\ |
1232 |
\frac{DS}{Dt} &=&\mathcal{Q}_{s} \label{eq-zcb-salt} |
1233 |
\end{eqnarray} |
1234 |
These equations still retain acoustic modes. But, because the |
1235 |
``compressible'' terms are linearized, the pressure equation \ref |
1236 |
{eq-zcb-cont} can be integrated implicitly with ease (the time-dependent |
1237 |
term appears as a Helmholtz term in the non-hydrostatic pressure equation). |
1238 |
These are the \emph{truly} compressible Boussinesq equations. Note that the |
1239 |
EOS must have the same pressure dependency as the linearized pressure term, |
1240 |
ie. $\left. \frac{\partial \rho }{\partial p}\right| _{\theta ,S}=\frac{1}{ |
1241 |
c_{s}^{2}}$, for consistency. |
1242 |
|
1243 |
\subsubsection{`Anelastic' z-coordinate equations} |
1244 |
|
1245 |
The anelastic approximation filters the acoustic mode by removing the |
1246 |
time-dependency in the continuity (now pressure-) equation (\ref{eq-zcb-cont} |
1247 |
). This could be done simply by noting that $\frac{Dp}{Dt}\approx -g\rho _{o} |
1248 |
\frac{Dz}{Dt}=-g\rho _{o}w$, but this leads to an inconsistency between |
1249 |
continuity and EOS. A better solution is to change the dependency on |
1250 |
pressure in the EOS by splitting the pressure into a reference function of |
1251 |
height and a perturbation: |
1252 |
\begin{equation*} |
1253 |
\rho =\rho (\theta ,S,p_{o}(z)+\epsilon _{s}p^{\prime }) |
1254 |
\end{equation*} |
1255 |
Remembering that the term $\frac{Dp}{Dt}$ in continuity comes from |
1256 |
differentiating the EOS, the continuity equation then becomes: |
1257 |
\begin{equation*} |
1258 |
\frac{1}{\rho _{o}c_{s}^{2}}\left( \frac{Dp_{o}}{Dt}+\epsilon _{s}\frac{ |
1259 |
Dp^{\prime }}{Dt}\right) +\mathbf{\nabla }_{z}\cdot \vec{\mathbf{v}}_{h}+ |
1260 |
\frac{\partial w}{\partial z}=0 |
1261 |
\end{equation*} |
1262 |
If the time- and space-scales of the motions of interest are longer than |
1263 |
those of acoustic modes, then $\frac{Dp^{\prime }}{Dt}<<(\frac{Dp_{o}}{Dt}, |
1264 |
\mathbf{\nabla }\cdot \vec{\mathbf{v}}_{h})$ in the continuity equations and |
1265 |
$\left. \frac{\partial \rho }{\partial p}\right| _{\theta ,S}\frac{ |
1266 |
Dp^{\prime }}{Dt}<<\left. \frac{\partial \rho }{\partial p}\right| _{\theta |
1267 |
,S}\frac{Dp_{o}}{Dt}$ in the EOS (\ref{EOSexpansion}). Thus we set $\epsilon |
1268 |
_{s}=0$, removing the dependency on $p^{\prime }$ in the continuity equation |
1269 |
and EOS. Expanding $\frac{Dp_{o}(z)}{Dt}=-g\rho _{o}w$ then leads to the |
1270 |
anelastic continuity equation: |
1271 |
\begin{equation} |
1272 |
\mathbf{\nabla }_{z}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial w}{\partial z}- |
1273 |
\frac{g}{c_{s}^{2}}w=0 \label{eq-za-cont1} |
1274 |
\end{equation} |
1275 |
A slightly different route leads to the quasi-Boussinesq continuity equation |
1276 |
where we use the scaling $\frac{\partial \rho ^{\prime }}{\partial t}+ |
1277 |
\mathbf{\nabla }_{3}\cdot \rho ^{\prime }\vec{\mathbf{v}}<<\mathbf{\nabla } |
1278 |
_{3}\cdot \rho _{o}\vec{\mathbf{v}}$ yielding: |
1279 |
\begin{equation} |
1280 |
\mathbf{\nabla }_{z}\cdot \vec{\mathbf{v}}_{h}+\frac{1}{\rho _{o}}\frac{ |
1281 |
\partial \left( \rho _{o}w\right) }{\partial z}=0 \label{eq-za-cont2} |
1282 |
\end{equation} |
1283 |
Equations \ref{eq-za-cont1} and \ref{eq-za-cont2} are in fact the same |
1284 |
equation if: |
1285 |
\begin{equation} |
1286 |
\frac{1}{\rho _{o}}\frac{\partial \rho _{o}}{\partial z}=\frac{-g}{c_{s}^{2}} |
1287 |
\end{equation} |
1288 |
Again, note that if $\rho _{o}$ is evaluated from prescribed $\theta _{o}$ |
1289 |
and $S_{o}$ profiles, then the EOS dependency on $p_{o}$ and the term $\frac{ |
1290 |
g}{c_{s}^{2}}$ in continuity should be referred to those same profiles. The |
1291 |
full set of `quasi-Boussinesq' or `anelastic' equations for the ocean are |
1292 |
then: |
1293 |
\begin{eqnarray} |
1294 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1295 |
_{h}+\frac{1}{\rho _{o}}\mathbf{\nabla }_{z}p &=&\vec{\mathbf{\mathcal{F}}} |
1296 |
\label{eq-zab-hmom} \\ |
1297 |
\epsilon _{nh}\frac{Dw}{Dt}+\frac{g\rho }{\rho _{o}}+\frac{1}{\rho _{o}} |
1298 |
\frac{\partial p}{\partial z} &=&\epsilon _{nh}\mathcal{F}_{w} |
1299 |
\label{eq-zab-hydro} \\ |
1300 |
\mathbf{\nabla }_{z}\cdot \vec{\mathbf{v}}_{h}+\frac{1}{\rho _{o}}\frac{ |
1301 |
\partial \left( \rho _{o}w\right) }{\partial z} &=&0 \label{eq-zab-cont} \\ |
1302 |
\rho &=&\rho (\theta ,S,p_{o}(z)) \label{eq-zab-eos} \\ |
1303 |
\frac{D\theta }{Dt} &=&\mathcal{Q}_{\theta } \label{eq-zab-heat} \\ |
1304 |
\frac{DS}{Dt} &=&\mathcal{Q}_{s} \label{eq-zab-salt} |
1305 |
\end{eqnarray} |
1306 |
|
1307 |
\subsubsection{Incompressible z-coordinate equations} |
1308 |
|
1309 |
Here, the objective is to drop the depth dependence of $\rho _{o}$ and so, |
1310 |
technically, to also remove the dependence of $\rho $ on $p_{o}$. This would |
1311 |
yield the ``truly'' incompressible Boussinesq equations: |
1312 |
\begin{eqnarray} |
1313 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1314 |
_{h}+\frac{1}{\rho _{c}}\mathbf{\nabla }_{z}p &=&\vec{\mathbf{\mathcal{F}}} |
1315 |
\label{eq-ztb-hmom} \\ |
1316 |
\epsilon _{nh}\frac{Dw}{Dt}+\frac{g\rho }{\rho _{c}}+\frac{1}{\rho _{c}} |
1317 |
\frac{\partial p}{\partial z} &=&\epsilon _{nh}\mathcal{F}_{w} |
1318 |
\label{eq-ztb-hydro} \\ |
1319 |
\mathbf{\nabla }_{z}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial w}{\partial z} |
1320 |
&=&0 \label{eq-ztb-cont} \\ |
1321 |
\rho &=&\rho (\theta ,S) \label{eq-ztb-eos} \\ |
1322 |
\frac{D\theta }{Dt} &=&\mathcal{Q}_{\theta } \label{eq-ztb-heat} \\ |
1323 |
\frac{DS}{Dt} &=&\mathcal{Q}_{s} \label{eq-ztb-salt} |
1324 |
\end{eqnarray} |
1325 |
where $\rho _{c}$ is a constant reference density of water. |
1326 |
|
1327 |
\subsubsection{Compressible non-divergent equations} |
1328 |
|
1329 |
The above ``incompressible'' equations are incompressible in both the flow |
1330 |
and the density. In many oceanic applications, however, it is important to |
1331 |
retain compressibility effects in the density. To do this we must split the |
1332 |
density thus: |
1333 |
\begin{equation*} |
1334 |
\rho =\rho _{o}+\rho ^{\prime } |
1335 |
\end{equation*} |
1336 |
We then assert that variations with depth of $\rho _{o}$ are unimportant |
1337 |
while the compressible effects in $\rho ^{\prime }$ are: |
1338 |
\begin{equation*} |
1339 |
\rho _{o}=\rho _{c} |
1340 |
\end{equation*} |
1341 |
\begin{equation*} |
1342 |
\rho ^{\prime }=\rho (\theta ,S,p_{o}(z))-\rho _{o} |
1343 |
\end{equation*} |
1344 |
This then yields what we can call the semi-compressible Boussinesq |
1345 |
equations: |
1346 |
\begin{eqnarray} |
1347 |
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}} |
1348 |
_{h}+\frac{1}{\rho _{c}}\mathbf{\nabla }_{z}p^{\prime } &=&\vec{\mathbf{ |
1349 |
\mathcal{F}}} \label{eq:ocean-mom} \\ |
1350 |
\epsilon _{nh}\frac{Dw}{Dt}+\frac{g\rho ^{\prime }}{\rho _{c}}+\frac{1}{\rho |
1351 |
_{c}}\frac{\partial p^{\prime }}{\partial z} &=&\epsilon _{nh}\mathcal{F}_{w} |
1352 |
\label{eq:ocean-wmom} \\ |
1353 |
\mathbf{\nabla }_{z}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial w}{\partial z} |
1354 |
&=&0 \label{eq:ocean-cont} \\ |
1355 |
\rho ^{\prime } &=&\rho (\theta ,S,p_{o}(z))-\rho _{c} \label{eq:ocean-eos} |
1356 |
\\ |
1357 |
\frac{D\theta }{Dt} &=&\mathcal{Q}_{\theta } \label{eq:ocean-theta} \\ |
1358 |
\frac{DS}{Dt} &=&\mathcal{Q}_{s} \label{eq:ocean-salt} |
1359 |
\end{eqnarray} |
1360 |
Note that the hydrostatic pressure of the resting fluid, including that |
1361 |
associated with $\rho _{c}$, is subtracted out since it has no effect on the |
1362 |
dynamics. |
1363 |
|
1364 |
Though necessary, the assumptions that go into these equations are messy |
1365 |
since we essentially assume a different EOS for the reference density and |
1366 |
the perturbation density. Nevertheless, it is the hydrostatic ($\epsilon |
1367 |
_{nh}=0$ form of these equations that are used throughout the ocean modeling |
1368 |
community and referred to as the primitive equations (HPE). |
1369 |
|
1370 |
% $Header: /u/gcmpack/mitgcmdoc/part1/manual.tex,v 1.3 2001/10/10 16:48:45 cnh Exp $ |
1371 |
% $Name: $ |
1372 |
|
1373 |
\section{Appendix:OPERATORS} |
1374 |
|
1375 |
\subsection{Coordinate systems} |
1376 |
|
1377 |
\subsubsection{Spherical coordinates} |
1378 |
|
1379 |
In spherical coordinates, the velocity components in the zonal, meridional |
1380 |
and vertical direction respectively, are given by (see Fig.2) : |
1381 |
|
1382 |
\begin{equation*} |
1383 |
u=r\cos \phi \frac{D\lambda }{Dt} |
1384 |
\end{equation*} |
1385 |
|
1386 |
\begin{equation*} |
1387 |
v=r\frac{D\phi }{Dt}\qquad |
1388 |
\end{equation*} |
1389 |
$\qquad \qquad \qquad \qquad $ |
1390 |
|
1391 |
\begin{equation*} |
1392 |
\dot{r}=\frac{Dr}{Dt} |
1393 |
\end{equation*} |
1394 |
|
1395 |
Here $\phi $ is the latitude, $\lambda $ the longitude, $r$ the radial |
1396 |
distance of the particle from the center of the earth, $\Omega $ is the |
1397 |
angular speed of rotation of the Earth and $D/Dt$ is the total derivative. |
1398 |
|
1399 |
The `grad' ($\nabla $) and `div' ($\nabla $.) operators are defined by, in |
1400 |
spherical coordinates: |
1401 |
|
1402 |
\begin{equation*} |
1403 |
\nabla \equiv \left( \frac{1}{r\cos \phi }\frac{\partial }{\partial \lambda } |
1404 |
,\frac{1}{r}\frac{\partial }{\partial \phi },\frac{\partial }{\partial r} |
1405 |
\right) |
1406 |
\end{equation*} |
1407 |
|
1408 |
\begin{equation*} |
1409 |
\nabla .v\equiv \frac{1}{r\cos \phi }\left\{ \frac{\partial u}{\partial |
1410 |
\lambda }+\frac{\partial }{\partial \phi }\left( v\cos \phi \right) \right\} |
1411 |
+\frac{1}{r^{2}}\frac{\partial \left( r^{2}\dot{r}\right) }{\partial r} |
1412 |
\end{equation*} |
1413 |
|
1414 |
\end{document} |