/[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.27 - (show annotations) (download) (as text)
Thu Jan 17 21:28:22 2008 UTC (17 years, 5 months ago) by jmc
Branch: MAIN
Changes since 1.26: +10 -8 lines
File MIME type: application/x-tex
add few references

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

  ViewVC Help
Powered by ViewVC 1.1.22