/[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.17 - (show annotations) (download) (as text)
Thu Aug 7 18:27:51 2003 UTC (21 years, 11 months ago) by edhill
Branch: MAIN
Changes since 1.16: +12 -12 lines
File MIME type: application/x-tex
Two improvements:
  - the manual build now works on a "dirty" (post-verification-runs)
    source tree
  - various small fixes to the manual (mostly Y. Hu's comments)

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

  ViewVC Help
Powered by ViewVC 1.1.22