/[MITgcm]/manual/s_overview/text/manual.tex
ViewVC logotype

Contents of /manual/s_overview/text/manual.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.25 - (show annotations) (download) (as text)
Sat Apr 8 01:50:49 2006 UTC (19 years, 3 months ago) by edhill
Branch: MAIN
Changes since 1.24: +44 -42 lines
File MIME type: application/x-tex
LaTeX syntax cleanups:
 - the degree symbols ("\circ") are now properly raised after latex2html
 - don't use $...$ when it should be \textit{...}

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

  ViewVC Help
Powered by ViewVC 1.1.22