8 |
centered finite difference) in the fluid interior but allows |
centered finite difference) in the fluid interior but allows |
9 |
boundaries to intersect a regular grid allowing a more accurate |
boundaries to intersect a regular grid allowing a more accurate |
10 |
representation of the position of the boundary. We treat the |
representation of the position of the boundary. We treat the |
11 |
horizontal and veritical directions as seperable and differently. |
horizontal and vertical directions as separable and differently. |
12 |
|
|
13 |
\input{part2/notation} |
\input{part2/notation} |
14 |
|
|
17 |
|
|
18 |
The finite volume method is used to discretize the equations in |
The finite volume method is used to discretize the equations in |
19 |
space. The expression ``finite volume'' actually has two meanings; one |
space. The expression ``finite volume'' actually has two meanings; one |
20 |
is the method of cut or instecting boundaries (shaved or lopped cells |
is the method of embedded or intersecting boundaries (shaved or lopped |
21 |
in our terminology) and the other is non-linear interpolation methods |
cells in our terminology) and the other is non-linear interpolation |
22 |
that can deal with non-smooth solutions such as shocks (i.e. flux |
methods that can deal with non-smooth solutions such as shocks |
23 |
limiters for advection). Both make use of the integral form of the |
(i.e. flux limiters for advection). Both make use of the integral form |
24 |
conservation laws to which the {\it weak solution} is a solution on |
of the conservation laws to which the {\it weak solution} is a |
25 |
each finite volume of (sub-domain). The weak solution can be |
solution on each finite volume of (sub-domain). The weak solution can |
26 |
constructed outof piece-wise constant elements or be |
be constructed out of piece-wise constant elements or be |
27 |
differentiable. The differentiable equations can not be satisfied by |
differentiable. The differentiable equations can not be satisfied by |
28 |
piece-wise constant functions. |
piece-wise constant functions. |
29 |
|
|
37 |
\begin{displaymath} |
\begin{displaymath} |
38 |
\Delta x \partial_t \theta + \delta_i ( F ) = 0 |
\Delta x \partial_t \theta + \delta_i ( F ) = 0 |
39 |
\end{displaymath} |
\end{displaymath} |
40 |
is exact if $\theta(x)$ is peice-wise constant over the interval |
is exact if $\theta(x)$ is piece-wise constant over the interval |
41 |
$\Delta x_i$ or more generally if $\theta_i$ is defined as the average |
$\Delta x_i$ or more generally if $\theta_i$ is defined as the average |
42 |
over the interval $\Delta x_i$. |
over the interval $\Delta x_i$. |
43 |
|
|
57 |
interior of a fluid. Differences arise at boundaries where a boundary |
interior of a fluid. Differences arise at boundaries where a boundary |
58 |
is not positioned on a regular or smoothly varying grid. This method |
is not positioned on a regular or smoothly varying grid. This method |
59 |
is used to represent the topography using lopped cell, see |
is used to represent the topography using lopped cell, see |
60 |
\cite{Adcroft98}. Subtle difference also appear in more than one |
\cite{adcroft:97}. Subtle difference also appear in more than one |
61 |
dimension away from boundaries. This happens because the each |
dimension away from boundaries. This happens because the each |
62 |
direction is discretized independantly in the finite difference method |
direction is discretized independently in the finite difference method |
63 |
while the integrating over finite volume implicitly treats all |
while the integrating over finite volume implicitly treats all |
64 |
directions simultaneously. Illustration of this is given in |
directions simultaneously. Illustration of this is given in |
65 |
\cite{Adcroft02}. |
\cite{ac:02}. |
66 |
|
|
67 |
\subsection{C grid staggering of variables} |
\subsection{C grid staggering of variables} |
68 |
|
|
69 |
\begin{figure} |
\begin{figure} |
70 |
\centerline{ \resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} } |
\begin{center} |
71 |
|
\resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} |
72 |
|
\end{center} |
73 |
\caption{Three dimensional staggering of velocity components. This |
\caption{Three dimensional staggering of velocity components. This |
74 |
facilitates the natural discretization of the continuity and tracer |
facilitates the natural discretization of the continuity and tracer |
75 |
equations. } |
equations. } |
79 |
The basic algorithm employed for stepping forward the momentum |
The basic algorithm employed for stepping forward the momentum |
80 |
equations is based on retaining non-divergence of the flow at all |
equations is based on retaining non-divergence of the flow at all |
81 |
times. This is most naturally done if the components of flow are |
times. This is most naturally done if the components of flow are |
82 |
staggered in space in the form of an Arakawa C grid \cite{Arakawa70}. |
staggered in space in the form of an Arakawa C grid \cite{arakawa:77}. |
83 |
|
|
84 |
Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$) |
Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$) |
85 |
staggered in space such that the zonal component falls on the |
staggered in space such that the zonal component falls on the |
86 |
interface between continiuty cells in the zonal direction. Similarly |
interface between continuity cells in the zonal direction. Similarly |
87 |
for the meridional and vertical directions. The continiuty cell is |
for the meridional and vertical directions. The continuity cell is |
88 |
synonymous with tracer cells (they are one and the same). |
synonymous with tracer cells (they are one and the same). |
89 |
|
|
90 |
|
|
115 |
|
|
116 |
|
|
117 |
\subsection{Horizontal grid} |
\subsection{Horizontal grid} |
118 |
|
\label{sec:spatial_discrete_horizontal_grid} |
119 |
|
|
120 |
\begin{figure} |
\begin{figure} |
121 |
\centerline{ \begin{tabular}{cc} |
\begin{center} |
122 |
|
\begin{tabular}{cc} |
123 |
\raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}} |
\raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}} |
124 |
& \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}} |
& \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}} |
125 |
\\ |
\\ |
126 |
\raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}} |
\raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}} |
127 |
& \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}} |
& \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}} |
128 |
\end{tabular} } |
\end{tabular} |
129 |
|
\end{center} |
130 |
\caption{ |
\caption{ |
131 |
Staggering of horizontal grid descriptors (lengths and areas). The |
Staggering of horizontal grid descriptors (lengths and areas). The |
132 |
grid lines indicate the tracer cell boundaries and are the reference |
grid lines indicate the tracer cell boundaries and are the reference |
141 |
|
|
142 |
The model domain is decomposed into tiles and within each tile a |
The model domain is decomposed into tiles and within each tile a |
143 |
quasi-regular grid is used. A tile is the basic unit of domain |
quasi-regular grid is used. A tile is the basic unit of domain |
144 |
decomposition for parallelization but may be used whether parallized |
decomposition for parallelization but may be used whether parallelized |
145 |
or not; see section \ref{sect:tiles} for more details. Although the |
or not; see section \ref{sect:tiles} for more details. Although the |
146 |
tiles may be patched together in an unstructured manner |
tiles may be patched together in an unstructured manner |
147 |
(i.e. irregular or non-tessilating pattern), the interior of tiles is |
(i.e. irregular or non-tessilating pattern), the interior of tiles is |
148 |
a structered grid of quadrilateral cells. The horizontal coordinate |
a structured grid of quadrilateral cells. The horizontal coordinate |
149 |
system is orthogonal curvilinear meaning we can not necessarily treat |
system is orthogonal curvilinear meaning we can not necessarily treat |
150 |
the two horizontal directions as seperable. Instead, each cell in the |
the two horizontal directions as separable. Instead, each cell in the |
151 |
horizontal grid is described by the length of it's sides and it's |
horizontal grid is described by the length of it's sides and it's |
152 |
area. |
area. |
153 |
|
|
154 |
The grid information is quite general and describes any of the |
The grid information is quite general and describes any of the |
155 |
available coordinates systems, cartesian, spherical-polar or |
available coordinates systems, cartesian, spherical-polar or |
156 |
curvilinear. All that is necessary to distinguish between the |
curvilinear. All that is necessary to distinguish between the |
157 |
coordinate systems is to initialize the grid data (discriptors) |
coordinate systems is to initialize the grid data (descriptors) |
158 |
appropriately. |
appropriately. |
159 |
|
|
160 |
In the following, we refer to the orientation of quantities on the |
In the following, we refer to the orientation of quantities on the |
293 |
spacing can be set to uniform via scalars {\bf dXspacing} and {\bf |
spacing can be set to uniform via scalars {\bf dXspacing} and {\bf |
294 |
dYspacing} in namelist {\em PARM04} or to variable resolution by the |
dYspacing} in namelist {\em PARM04} or to variable resolution by the |
295 |
vectors {\bf DELX} and {\bf DELY}. Units are normally |
vectors {\bf DELX} and {\bf DELY}. Units are normally |
296 |
meters. Non-dimensional coordinates can be used by interpretting the |
meters. Non-dimensional coordinates can be used by interpreting the |
297 |
gravitational constant as the Rayleigh number. |
gravitational constant as the Rayleigh number. |
298 |
|
|
299 |
\subsubsection{Spherical-polar coordinates} |
\subsubsection{Spherical-polar coordinates} |
318 |
\subsection{Vertical grid} |
\subsection{Vertical grid} |
319 |
|
|
320 |
\begin{figure} |
\begin{figure} |
321 |
\centerline{ \begin{tabular}{cc} |
\begin{center} |
322 |
|
\begin{tabular}{cc} |
323 |
\raisebox{4in}{a)} \resizebox{!}{4in}{ |
\raisebox{4in}{a)} \resizebox{!}{4in}{ |
324 |
\includegraphics{part2/vgrid-cellcentered.eps}} & \raisebox{4in}{b)} |
\includegraphics{part2/vgrid-cellcentered.eps}} & \raisebox{4in}{b)} |
325 |
\resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}} |
\resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}} |
326 |
\end{tabular} } |
\end{tabular} |
327 |
|
\end{center} |
328 |
\caption{Two versions of the vertical grid. a) The cell centered |
\caption{Two versions of the vertical grid. a) The cell centered |
329 |
approach where the interface depths are specified and the tracer |
approach where the interface depths are specified and the tracer |
330 |
points centered in between the interfaces. b) The interface centered |
points centered in between the interfaces. b) The interface centered |
347 |
INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in |
INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in |
348 |
namelist {\em PARM04}. The units of ``r'' are either meters or Pascals |
namelist {\em PARM04}. The units of ``r'' are either meters or Pascals |
349 |
depending on the isomorphism being used which in turn is dependent |
depending on the isomorphism being used which in turn is dependent |
350 |
only on the choise of equation of state. |
only on the choice of equation of state. |
351 |
|
|
352 |
There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which |
There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which |
353 |
dictate whether z- or |
dictate whether z- or |
389 |
\subsection{Topography: partially filled cells} |
\subsection{Topography: partially filled cells} |
390 |
|
|
391 |
\begin{figure} |
\begin{figure} |
392 |
\centerline{ |
\begin{center} |
393 |
\resizebox{4.5in}{!}{\includegraphics{part2/vgrid-xz.eps}} |
\resizebox{4.5in}{!}{\includegraphics{part2/vgrid-xz.eps}} |
394 |
} |
\end{center} |
395 |
\caption{ |
\caption{ |
396 |
A schematic of the x-r plane showing the location of the |
A schematic of the x-r plane showing the location of the |
397 |
non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a |
non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a |
400 |
\label{fig:hfacs} |
\label{fig:hfacs} |
401 |
\end{figure} |
\end{figure} |
402 |
|
|
403 |
\cite{Adcroft97} presented two alternatives to the step-wise finite |
\cite{adcroft:97} presented two alternatives to the step-wise finite |
404 |
difference representation of topography. The method is known to the |
difference representation of topography. The method is known to the |
405 |
engineering community as {\em intersecting boundary method}. It |
engineering community as {\em intersecting boundary method}. It |
406 |
involves allowing the boundary to intersect a grid of cells thereby |
involves allowing the boundary to intersect a grid of cells thereby |
407 |
modifying the shape of those cells intersected. We suggested allowing |
modifying the shape of those cells intersected. We suggested allowing |
408 |
the topgoraphy to take on a peice-wise linear representation (shaved |
the topography to take on a piece-wise linear representation (shaved |
409 |
cells) or a simpler piecewise constant representaion (partial step). |
cells) or a simpler piecewise constant representation (partial step). |
410 |
Both show dramatic improvements in solution compared to the |
Both show dramatic improvements in solution compared to the |
411 |
traditional full step representation, the piece-wise linear being the |
traditional full step representation, the piece-wise linear being the |
412 |
best. However, the storage requirements are excessive so the simpler |
best. However, the storage requirements are excessive so the simpler |
420 |
\marginpar{$h_s$: {\bf hFacS}} |
\marginpar{$h_s$: {\bf hFacS}} |
421 |
The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta |
The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta |
422 |
r_f(k)$ and the physical thickness of the open side is given by |
r_f(k)$ and the physical thickness of the open side is given by |
423 |
$h_w(i,j,k) \Delta r_f(k)$. Three 3-D discriptors $h_c$, $h_w$ and |
$h_w(i,j,k) \Delta r_f(k)$. Three 3-D descriptors $h_c$, $h_w$ and |
424 |
$h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and |
$h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and |
425 |
{\bf hFacS} respectively. These are calculated in subroutine {\em |
{\bf hFacS} respectively. These are calculated in subroutine {\em |
426 |
INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf |
INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf |
461 |
The core algorithm is based on the ``C grid'' discretization of the |
The core algorithm is based on the ``C grid'' discretization of the |
462 |
continuity equation which can be summarized as: |
continuity equation which can be summarized as: |
463 |
\begin{eqnarray} |
\begin{eqnarray} |
464 |
\partial_t u + \frac{1}{\Delta x_c} \delta_i \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta x_c} \delta_i \Phi_{nh}' & = & G_u - \frac{1}{\Delta x_c} \delta_i \Phi_h' \\ |
\partial_t u + \frac{1}{\Delta x_c} \delta_i \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta x_c} \delta_i \Phi_{nh}' & = & G_u - \frac{1}{\Delta x_c} \delta_i \Phi_h' \label{eq:discrete-momu} \\ |
465 |
\partial_t v + \frac{1}{\Delta y_c} \delta_j \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta y_c} \delta_j \Phi_{nh}' & = & G_v - \frac{1}{\Delta y_c} \delta_j \Phi_h' \\ |
\partial_t v + \frac{1}{\Delta y_c} \delta_j \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta y_c} \delta_j \Phi_{nh}' & = & G_v - \frac{1}{\Delta y_c} \delta_j \Phi_h' \label{eq:discrete-momv} \\ |
466 |
\epsilon_{nh} \left( \partial_t w + \frac{1}{\Delta r_c} \delta_k \Phi_{nh}' \right) & = & \epsilon_{nh} G_w + \overline{b}^k - \frac{1}{\Delta r_c} \delta_k \Phi_{h}' \\ |
\epsilon_{nh} \left( \partial_t w + \frac{1}{\Delta r_c} \delta_k \Phi_{nh}' \right) & = & \epsilon_{nh} G_w + \overline{b}^k - \frac{1}{\Delta r_c} \delta_k \Phi_{h}' \label{eq:discrete-momw} \\ |
467 |
\delta_i \Delta y_g \Delta r_f h_w u + |
\delta_i \Delta y_g \Delta r_f h_w u + |
468 |
\delta_j \Delta x_g \Delta r_f h_s v + |
\delta_j \Delta x_g \Delta r_f h_s v + |
469 |
\delta_k {\cal A}_c w & = & {\cal A}_c \delta_k (P-E)_{r=0} |
\delta_k {\cal A}_c w & = & {\cal A}_c \delta_k (P-E)_{r=0} |
471 |
\end{eqnarray} |
\end{eqnarray} |
472 |
where the continuity equation has been most naturally discretized by |
where the continuity equation has been most naturally discretized by |
473 |
staggering the three components of velocity as shown in |
staggering the three components of velocity as shown in |
474 |
Fig.~\ref{fig-cgrid3d}. The grid lengths $\Delta x_c$ and $\Delta y_c$ |
Fig.~\ref{fig:cgrid3d}. The grid lengths $\Delta x_c$ and $\Delta y_c$ |
475 |
are the lengths between tracer points (cell centers). The grid lengths |
are the lengths between tracer points (cell centers). The grid lengths |
476 |
$\Delta x_g$, $\Delta y_g$ are the grid lengths between cell |
$\Delta x_g$, $\Delta y_g$ are the grid lengths between cell |
477 |
corners. $\Delta r_f$ and $\Delta r_c$ are the distance (in units of |
corners. $\Delta r_f$ and $\Delta r_c$ are the distance (in units of |
484 |
\marginpar{$h_s$: {\bf hFacS}} |
\marginpar{$h_s$: {\bf hFacS}} |
485 |
|
|
486 |
The last equation, the discrete continuity equation, can be summed in |
The last equation, the discrete continuity equation, can be summed in |
487 |
the vertical to yeild the free-surface equation: |
the vertical to yield the free-surface equation: |
488 |
\begin{equation} |
\begin{equation} |
489 |
{\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = |
{\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w |
490 |
{\cal A}_c(P-E)_{r=0} |
u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal |
491 |
|
A}_c(P-E)_{r=0} \label{eq:discrete-freesurface} |
492 |
\end{equation} |
\end{equation} |
493 |
The source term $P-E$ on the rhs of continuity accounts for the local |
The source term $P-E$ on the rhs of continuity accounts for the local |
494 |
addition of volume due to excess precipitation and run-off over |
addition of volume due to excess precipitation and run-off over |
503 |
from the pressure gradient terms when forming the kinetic energy |
from the pressure gradient terms when forming the kinetic energy |
504 |
equation. |
equation. |
505 |
|
|
506 |
In the ocean, using z-ccordinates, the hydrostatic balance terms are |
In the ocean, using z-coordinates, the hydrostatic balance terms are |
507 |
discretized: |
discretized: |
508 |
\begin{equation} |
\begin{equation} |
509 |
\epsilon_{nh} \partial_t w |
\epsilon_{nh} \partial_t w |
510 |
+ g \overline{\rho'}^k + \frac{1}{\Delta z} \delta_k \Phi_h' = \ldots |
+ g \overline{\rho'}^k + \frac{1}{\Delta z} \delta_k \Phi_h' = \ldots |
511 |
|
\label{eq:discrete_hydro_ocean} |
512 |
\end{equation} |
\end{equation} |
513 |
|
|
514 |
In the atmosphere, using p-coordinates, hydrostatic balance is |
In the atmosphere, using p-coordinates, hydrostatic balance is |
515 |
discretized: |
discretized: |
516 |
\begin{equation} |
\begin{equation} |
517 |
\overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0 |
\overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0 |
518 |
|
\label{eq:discrete_hydro_atmos} |
519 |
\end{equation} |
\end{equation} |
520 |
where $\Delta \Pi$ is the difference in Exner function between the |
where $\Delta \Pi$ is the difference in Exner function between the |
521 |
pressure points. The non-hydrostatic equations are not available in |
pressure points. The non-hydrostatic equations are not available in |
523 |
|
|
524 |
The difference in approach between ocean and atmosphere occurs because |
The difference in approach between ocean and atmosphere occurs because |
525 |
of the direct use of the ideal gas equation in forming the potential |
of the direct use of the ideal gas equation in forming the potential |
526 |
energy conversion term $\alpha \omega$. The form of these consversion |
energy conversion term $\alpha \omega$. The form of these conversion |
527 |
terms is discussed at length in \cite{Adcroft01}. |
terms is discussed at length in \cite{adcroft:02}. |
528 |
|
|
529 |
Because of the different representation of hydrostatic balance between |
Because of the different representation of hydrostatic balance between |
530 |
ocean and atmosphere there is no elegant way to represent both systems |
ocean and atmosphere there is no elegant way to represent both systems |