/[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.15 - (show annotations) (download) (as text)
Wed Nov 21 16:33:17 2001 UTC (23 years, 7 months ago) by cnh
Branch: MAIN
Changes since 1.14: +11 -10 lines
File MIME type: application/x-tex
Caoption and figure reference updates

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

  ViewVC Help
Powered by ViewVC 1.1.22