1 |
% $Header: /u/gcmpack/manual/part2/spatial-discrete.tex,v 1.14 2004/03/23 16:47:04 afe Exp $ |
2 |
% $Name: $ |
3 |
|
4 |
\section{Spatial discretization of the dynamical equations} |
5 |
|
6 |
Spatial discretization is carried out using the finite volume |
7 |
method. This amounts to a grid-point method (namely second-order |
8 |
centered finite difference) in the fluid interior but allows |
9 |
boundaries to intersect a regular grid allowing a more accurate |
10 |
representation of the position of the boundary. We treat the |
11 |
horizontal and vertical directions as separable and differently. |
12 |
|
13 |
\input{part2/notation} |
14 |
|
15 |
|
16 |
\subsection{The finite volume method: finite volumes versus finite difference} |
17 |
\begin{rawhtml} |
18 |
<!-- CMIREDIR:finite_volume: --> |
19 |
\end{rawhtml} |
20 |
|
21 |
|
22 |
|
23 |
The finite volume method is used to discretize the equations in |
24 |
space. The expression ``finite volume'' actually has two meanings; one |
25 |
is the method of embedded or intersecting boundaries (shaved or lopped |
26 |
cells in our terminology) and the other is non-linear interpolation |
27 |
methods that can deal with non-smooth solutions such as shocks |
28 |
(i.e. flux limiters for advection). Both make use of the integral form |
29 |
of the conservation laws to which the {\it weak solution} is a |
30 |
solution on each finite volume of (sub-domain). The weak solution can |
31 |
be constructed out of piece-wise constant elements or be |
32 |
differentiable. The differentiable equations can not be satisfied by |
33 |
piece-wise constant functions. |
34 |
|
35 |
As an example, the 1-D constant coefficient advection-diffusion |
36 |
equation: |
37 |
\begin{displaymath} |
38 |
\partial_t \theta + \partial_x ( u \theta - \kappa \partial_x \theta ) = 0 |
39 |
\end{displaymath} |
40 |
can be discretized by integrating over finite sub-domains, i.e. |
41 |
the lengths $\Delta x_i$: |
42 |
\begin{displaymath} |
43 |
\Delta x \partial_t \theta + \delta_i ( F ) = 0 |
44 |
\end{displaymath} |
45 |
is exact if $\theta(x)$ is piece-wise constant over the interval |
46 |
$\Delta x_i$ or more generally if $\theta_i$ is defined as the average |
47 |
over the interval $\Delta x_i$. |
48 |
|
49 |
The flux, $F_{i-1/2}$, must be approximated: |
50 |
\begin{displaymath} |
51 |
F = u \overline{\theta} - \frac{\kappa}{\Delta x_c} \partial_i \theta |
52 |
\end{displaymath} |
53 |
and this is where truncation errors can enter the solution. The |
54 |
method for obtaining $\overline{\theta}$ is unspecified and a wide |
55 |
range of possibilities exist including centered and upwind |
56 |
interpolation, polynomial fits based on the the volume average |
57 |
definitions of quantities and non-linear interpolation such as |
58 |
flux-limiters. |
59 |
|
60 |
Choosing simple centered second-order interpolation and differencing |
61 |
recovers the same ODE's resulting from finite differencing for the |
62 |
interior of a fluid. Differences arise at boundaries where a boundary |
63 |
is not positioned on a regular or smoothly varying grid. This method |
64 |
is used to represent the topography using lopped cell, see |
65 |
\cite{adcroft:97}. Subtle difference also appear in more than one |
66 |
dimension away from boundaries. This happens because the each |
67 |
direction is discretized independently in the finite difference method |
68 |
while the integrating over finite volume implicitly treats all |
69 |
directions simultaneously. Illustration of this is given in |
70 |
\cite{ac:02}. |
71 |
|
72 |
\subsection{C grid staggering of variables} |
73 |
|
74 |
\begin{figure} |
75 |
\begin{center} |
76 |
\resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} |
77 |
\end{center} |
78 |
\caption{Three dimensional staggering of velocity components. This |
79 |
facilitates the natural discretization of the continuity and tracer |
80 |
equations. } |
81 |
\label{fig:cgrid3d} |
82 |
\end{figure} |
83 |
|
84 |
The basic algorithm employed for stepping forward the momentum |
85 |
equations is based on retaining non-divergence of the flow at all |
86 |
times. This is most naturally done if the components of flow are |
87 |
staggered in space in the form of an Arakawa C grid \cite{arakawa:77}. |
88 |
|
89 |
Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$) |
90 |
staggered in space such that the zonal component falls on the |
91 |
interface between continuity cells in the zonal direction. Similarly |
92 |
for the meridional and vertical directions. The continuity cell is |
93 |
synonymous with tracer cells (they are one and the same). |
94 |
|
95 |
|
96 |
|
97 |
\subsection{Grid initialization and data} |
98 |
|
99 |
Initialization of grid data is controlled by subroutine {\em |
100 |
INI\_GRID} which in calls {\em INI\_VERTICAL\_GRID} to initialize the |
101 |
vertical grid, and then either of {\em INI\_CARTESIAN\_GRID}, {\em |
102 |
INI\_SPHERICAL\_POLAR\_GRID} or {\em INI\_CURV\-ILINEAR\_GRID} to |
103 |
initialize the horizontal grid for cartesian, spherical-polar or |
104 |
curvilinear coordinates respectively. |
105 |
|
106 |
The reciprocals of all grid quantities are pre-calculated and this is |
107 |
done in subroutine {\em INI\_MASKS\_ETC} which is called later by |
108 |
subroutine {\em INITIALIZE\_FIXED}. |
109 |
|
110 |
All grid descriptors are global arrays and stored in common blocks in |
111 |
{\em GRID.h} and a generally declared as {\em \_RS}. |
112 |
|
113 |
\fbox{ \begin{minipage}{4.75in} |
114 |
{\em S/R INI\_GRID} ({\em model/src/ini\_grid.F}) |
115 |
|
116 |
{\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F}) |
117 |
|
118 |
grid data: ({\em model/inc/GRID.h}) |
119 |
\end{minipage} } |
120 |
|
121 |
|
122 |
\subsection{Horizontal grid} |
123 |
\label{sec:spatial_discrete_horizontal_grid} |
124 |
|
125 |
\begin{figure} |
126 |
\begin{center} |
127 |
\begin{tabular}{cc} |
128 |
\raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}} |
129 |
& \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}} |
130 |
\\ |
131 |
\raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}} |
132 |
& \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}} |
133 |
\end{tabular} |
134 |
\end{center} |
135 |
\caption{ |
136 |
Staggering of horizontal grid descriptors (lengths and areas). The |
137 |
grid lines indicate the tracer cell boundaries and are the reference |
138 |
grid for all panels. a) The area of a tracer cell, $A_c$, is bordered |
139 |
by the lengths $\Delta x_g$ and $\Delta y_g$. b) The area of a |
140 |
vorticity cell, $A_\zeta$, is bordered by the lengths $\Delta x_c$ and |
141 |
$\Delta y_c$. c) The area of a u cell, $A_c$, is bordered by the |
142 |
lengths $\Delta x_v$ and $\Delta y_f$. d) The area of a v cell, $A_c$, |
143 |
is bordered by the lengths $\Delta x_f$ and $\Delta y_u$.} |
144 |
\label{fig:hgrid} |
145 |
\end{figure} |
146 |
|
147 |
The model domain is decomposed into tiles and within each tile a |
148 |
quasi-regular grid is used. A tile is the basic unit of domain |
149 |
decomposition for parallelization but may be used whether parallelized |
150 |
or not; see section \ref{sect:tiles} for more details. Although the |
151 |
tiles may be patched together in an unstructured manner |
152 |
(i.e. irregular or non-tessilating pattern), the interior of tiles is |
153 |
a structured grid of quadrilateral cells. The horizontal coordinate |
154 |
system is orthogonal curvilinear meaning we can not necessarily treat |
155 |
the two horizontal directions as separable. Instead, each cell in the |
156 |
horizontal grid is described by the length of it's sides and it's |
157 |
area. |
158 |
|
159 |
The grid information is quite general and describes any of the |
160 |
available coordinates systems, cartesian, spherical-polar or |
161 |
curvilinear. All that is necessary to distinguish between the |
162 |
coordinate systems is to initialize the grid data (descriptors) |
163 |
appropriately. |
164 |
|
165 |
In the following, we refer to the orientation of quantities on the |
166 |
computational grid using geographic terminology such as points of the |
167 |
compass. |
168 |
\marginpar{Caution!} |
169 |
This is purely for convenience but should note be confused |
170 |
with the actual geographic orientation of model quantities. |
171 |
|
172 |
Fig.~\ref{fig:hgrid}a shows the tracer cell (synonymous with the |
173 |
continuity cell). The length of the southern edge, $\Delta x_g$, |
174 |
western edge, $\Delta y_g$ and surface area, $A_c$, presented in the |
175 |
vertical are stored in arrays {\bf DXg}, {\bf DYg} and {\bf rAc}. |
176 |
\marginpar{$A_c$: {\bf rAc}} |
177 |
\marginpar{$\Delta x_g$: {\bf DXg}} |
178 |
\marginpar{$\Delta y_g$: {\bf DYg}} |
179 |
The ``g'' suffix indicates that the lengths are along the defining |
180 |
grid boundaries. The ``c'' suffix associates the quantity with the |
181 |
cell centers. The quantities are staggered in space and the indexing |
182 |
is such that {\bf DXg(i,j)} is positioned to the south of {\bf |
183 |
rAc(i,j)} and {\bf DYg(i,j)} positioned to the west. |
184 |
|
185 |
Fig.~\ref{fig:hgrid}b shows the vorticity cell. The length of the |
186 |
southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface |
187 |
area, $A_\zeta$, presented in the vertical are stored in arrays {\bf |
188 |
DXg}, {\bf DYg} and {\bf rAz}. |
189 |
\marginpar{$A_\zeta$: {\bf rAz}} |
190 |
\marginpar{$\Delta x_c$: {\bf DXc}} |
191 |
\marginpar{$\Delta y_c$: {\bf DYc}} |
192 |
The ``z'' suffix indicates that the lengths are measured between the |
193 |
cell centers and the ``$\zeta$'' suffix associates points with the |
194 |
vorticity points. The quantities are staggered in space and the |
195 |
indexing is such that {\bf DXc(i,j)} is positioned to the north of |
196 |
{\bf rAc(i,j)} and {\bf DYc(i,j)} positioned to the east. |
197 |
|
198 |
Fig.~\ref{fig:hgrid}c shows the ``u'' or western (w) cell. The length of |
199 |
the southern edge, $\Delta x_v$, eastern edge, $\Delta y_f$ and |
200 |
surface area, $A_w$, presented in the vertical are stored in arrays |
201 |
{\bf DXv}, {\bf DYf} and {\bf rAw}. |
202 |
\marginpar{$A_w$: {\bf rAw}} |
203 |
\marginpar{$\Delta x_v$: {\bf DXv}} |
204 |
\marginpar{$\Delta y_f$: {\bf DYf}} |
205 |
The ``v'' suffix indicates that the length is measured between the |
206 |
v-points, the ``f'' suffix indicates that the length is measured |
207 |
between the (tracer) cell faces and the ``w'' suffix associates points |
208 |
with the u-points (w stands for west). The quantities are staggered in |
209 |
space and the indexing is such that {\bf DXv(i,j)} is positioned to |
210 |
the south of {\bf rAw(i,j)} and {\bf DYf(i,j)} positioned to the east. |
211 |
|
212 |
Fig.~\ref{fig:hgrid}d shows the ``v'' or southern (s) cell. The length of |
213 |
the northern edge, $\Delta x_f$, western edge, $\Delta y_u$ and |
214 |
surface area, $A_s$, presented in the vertical are stored in arrays |
215 |
{\bf DXf}, {\bf DYu} and {\bf rAs}. |
216 |
\marginpar{$A_s$: {\bf rAs}} |
217 |
\marginpar{$\Delta x_f$: {\bf DXf}} |
218 |
\marginpar{$\Delta y_u$: {\bf DYu}} |
219 |
The ``u'' suffix indicates that the length is measured between the |
220 |
u-points, the ``f'' suffix indicates that the length is measured |
221 |
between the (tracer) cell faces and the ``s'' suffix associates points |
222 |
with the v-points (s stands for south). The quantities are staggered |
223 |
in space and the indexing is such that {\bf DXf(i,j)} is positioned to |
224 |
the north of {\bf rAs(i,j)} and {\bf DYu(i,j)} positioned to the west. |
225 |
|
226 |
\fbox{ \begin{minipage}{4.75in} |
227 |
{\em S/R INI\_CARTESIAN\_GRID} ({\em |
228 |
model/src/ini\_cartesian\_grid.F}) |
229 |
|
230 |
{\em S/R INI\_SPHERICAL\_POLAR\_GRID} ({\em |
231 |
model/src/ini\_spherical\_polar\_grid.F}) |
232 |
|
233 |
{\em S/R INI\_CURVILINEAR\_GRID} ({\em |
234 |
model/src/ini\_curvilinear\_grid.F}) |
235 |
|
236 |
$A_c$, $A_\zeta$, $A_w$, $A_s$: {\bf rAc}, {\bf rAz}, {\bf rAw}, {\bf rAs} |
237 |
({\em GRID.h}) |
238 |
|
239 |
$\Delta x_g$, $\Delta y_g$: {\bf DXg}, {\bf DYg} ({\em GRID.h}) |
240 |
|
241 |
$\Delta x_c$, $\Delta y_c$: {\bf DXc}, {\bf DYc} ({\em GRID.h}) |
242 |
|
243 |
$\Delta x_f$, $\Delta y_f$: {\bf DXf}, {\bf DYf} ({\em GRID.h}) |
244 |
|
245 |
$\Delta x_v$, $\Delta y_u$: {\bf DXv}, {\bf DYu} ({\em GRID.h}) |
246 |
|
247 |
\end{minipage} } |
248 |
|
249 |
\subsubsection{Reciprocals of horizontal grid descriptors} |
250 |
|
251 |
%\marginpar{$A_c^{-1}$: {\bf RECIP\_rAc}} |
252 |
%\marginpar{$A_\zeta^{-1}$: {\bf RECIP\_rAz}} |
253 |
%\marginpar{$A_w^{-1}$: {\bf RECIP\_rAw}} |
254 |
%\marginpar{$A_s^{-1}$: {\bf RECIP\_rAs}} |
255 |
Lengths and areas appear in the denominator of expressions as much as |
256 |
in the numerator. For efficiency and portability, we pre-calculate the |
257 |
reciprocal of the horizontal grid quantities so that in-line divisions |
258 |
can be avoided. |
259 |
|
260 |
%\marginpar{$\Delta x_g^{-1}$: {\bf RECIP\_DXg}} |
261 |
%\marginpar{$\Delta y_g^{-1}$: {\bf RECIP\_DYg}} |
262 |
%\marginpar{$\Delta x_c^{-1}$: {\bf RECIP\_DXc}} |
263 |
%\marginpar{$\Delta y_c^{-1}$: {\bf RECIP\_DYc}} |
264 |
%\marginpar{$\Delta x_f^{-1}$: {\bf RECIP\_DXf}} |
265 |
%\marginpar{$\Delta y_f^{-1}$: {\bf RECIP\_DYf}} |
266 |
%\marginpar{$\Delta x_v^{-1}$: {\bf RECIP\_DXv}} |
267 |
%\marginpar{$\Delta y_u^{-1}$: {\bf RECIP\_DYu}} |
268 |
For each grid descriptor (array) there is a reciprocal named using the |
269 |
prefix {\bf RECIP\_}. This doubles the amount of storage in {\em |
270 |
GRID.h} but they are all only 2-D descriptors. |
271 |
|
272 |
\fbox{ \begin{minipage}{4.75in} |
273 |
{\em S/R INI\_MASKS\_ETC} ({\em |
274 |
model/src/ini\_masks\_etc.F}) |
275 |
|
276 |
$A_c^{-1}$: {\bf RECIP\_Ac} ({\em GRID.h}) |
277 |
|
278 |
$A_\zeta^{-1}$: {\bf RECIP\_Az} ({\em GRID.h}) |
279 |
|
280 |
$A_w^{-1}$: {\bf RECIP\_Aw} ({\em GRID.h}) |
281 |
|
282 |
$A_s^{-1}$: {\bf RECIP\_As} ({\em GRID.h}) |
283 |
|
284 |
$\Delta x_g^{-1}$, $\Delta y_g^{-1}$: {\bf RECIP\_DXg}, {\bf RECIP\_DYg} ({\em GRID.h}) |
285 |
|
286 |
$\Delta x_c^{-1}$, $\Delta y_c^{-1}$: {\bf RECIP\_DXc}, {\bf RECIP\_DYc} ({\em GRID.h}) |
287 |
|
288 |
$\Delta x_f^{-1}$, $\Delta y_f^{-1}$: {\bf RECIP\_DXf}, {\bf RECIP\_DYf} ({\em GRID.h}) |
289 |
|
290 |
$\Delta x_v^{-1}$, $\Delta y_u^{-1}$: {\bf RECIP\_DXv}, {\bf RECIP\_DYu} ({\em GRID.h}) |
291 |
|
292 |
\end{minipage} } |
293 |
|
294 |
\subsubsection{Cartesian coordinates} |
295 |
|
296 |
Cartesian coordinates are selected when the logical flag {\bf |
297 |
using\-Cartes\-ianGrid} in namelist {\em PARM04} is set to true. The grid |
298 |
spacing can be set to uniform via scalars {\bf dXspacing} and {\bf |
299 |
dYspacing} in namelist {\em PARM04} or to variable resolution by the |
300 |
vectors {\bf DELX} and {\bf DELY}. Units are normally |
301 |
meters. Non-dimensional coordinates can be used by interpreting the |
302 |
gravitational constant as the Rayleigh number. |
303 |
|
304 |
\subsubsection{Spherical-polar coordinates} |
305 |
|
306 |
Spherical coordinates are selected when the logical flag {\bf |
307 |
using\-Spherical\-PolarGrid} in namelist {\em PARM04} is set to true. The |
308 |
grid spacing can be set to uniform via scalars {\bf dXspacing} and |
309 |
{\bf dYspacing} in namelist {\em PARM04} or to variable resolution by |
310 |
the vectors {\bf DELX} and {\bf DELY}. Units of these namelist |
311 |
variables are alway degrees. The horizontal grid descriptors are |
312 |
calculated from these namelist variables have units of meters. |
313 |
|
314 |
\subsubsection{Curvilinear coordinates} |
315 |
|
316 |
Curvilinear coordinates are selected when the logical flag {\bf |
317 |
using\-Curvil\-inear\-Grid} in namelist {\em PARM04} is set to true. The |
318 |
grid spacing can not be set via the namelist. Instead, the grid |
319 |
descriptors are read from data files, one for each descriptor. As for |
320 |
other grids, the horizontal grid descriptors have units of meters. |
321 |
|
322 |
|
323 |
\subsection{Vertical grid} |
324 |
|
325 |
\begin{figure} |
326 |
\begin{center} |
327 |
\begin{tabular}{cc} |
328 |
\raisebox{4in}{a)} \resizebox{!}{4in}{ |
329 |
\includegraphics{part2/vgrid-cellcentered.eps}} & \raisebox{4in}{b)} |
330 |
\resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}} |
331 |
\end{tabular} |
332 |
\end{center} |
333 |
\caption{Two versions of the vertical grid. a) The cell centered |
334 |
approach where the interface depths are specified and the tracer |
335 |
points centered in between the interfaces. b) The interface centered |
336 |
approach where tracer levels are specified and the w-interfaces are |
337 |
centered in between.} |
338 |
\label{fig:vgrid} |
339 |
\end{figure} |
340 |
|
341 |
As for the horizontal grid, we use the suffixes ``c'' and ``f'' to |
342 |
indicates faces and centers. Fig.~\ref{fig:vgrid}a shows the default |
343 |
vertical grid used by the model. |
344 |
\marginpar{$\Delta r_f$: {\bf DRf}} |
345 |
\marginpar{$\Delta r_c$: {\bf DRc}} |
346 |
$\Delta r_f$ is the difference in $r$ |
347 |
(vertical coordinate) between the faces (i.e. $\Delta r_f \equiv - |
348 |
\delta_k r$ where the minus sign appears due to the convention that the |
349 |
surface layer has index $k=1$.). |
350 |
|
351 |
The vertical grid is calculated in subroutine {\em |
352 |
INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in |
353 |
namelist {\em PARM04}. The units of ``r'' are either meters or Pascals |
354 |
depending on the isomorphism being used which in turn is dependent |
355 |
only on the choice of equation of state. |
356 |
|
357 |
There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which |
358 |
dictate whether z- or |
359 |
\marginpar{Caution!} |
360 |
p- coordinates are to be used but we intend to |
361 |
phase this out since they are redundant. |
362 |
|
363 |
The reciprocals $\Delta r_f^{-1}$ and $\Delta r_c^{-1}$ are |
364 |
pre-calculated (also in subroutine {\em INI\_VERTICAL\_GRID}). All |
365 |
vertical grid descriptors are stored in common blocks in {\em GRID.h}. |
366 |
|
367 |
The above grid (Fig.~\ref{fig:vgrid}a) is known as the cell centered |
368 |
approach because the tracer points are at cell centers; the cell |
369 |
centers are mid-way between the cell interfaces. An alternative, the |
370 |
vertex or interface centered approach, is shown in |
371 |
Fig.~\ref{fig:vgrid}b. Here, the interior interfaces are positioned |
372 |
mid-way between the tracer nodes (no longer cell centers). This |
373 |
approach is formally more accurate for evaluation of hydrostatic |
374 |
pressure and vertical advection but historically the cell centered |
375 |
approach has been used. An alternative form of subroutine {\em |
376 |
INI\_VERTICAL\_GRID} is used to select the interface centered approach |
377 |
but no run time option is currently available. |
378 |
|
379 |
\fbox{ \begin{minipage}{4.75in} |
380 |
{\em S/R INI\_VERTICAL\_GRID} ({\em |
381 |
model/src/ini\_vertical\_grid.F}) |
382 |
|
383 |
$\Delta r_f$: {\bf DRf} ({\em GRID.h}) |
384 |
|
385 |
$\Delta r_c$: {\bf DRc} ({\em GRID.h}) |
386 |
|
387 |
$\Delta r_f^{-1}$: {\bf RECIP\_DRf} ({\em GRID.h}) |
388 |
|
389 |
$\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\em GRID.h}) |
390 |
|
391 |
\end{minipage} } |
392 |
|
393 |
|
394 |
\subsection{Topography: partially filled cells} |
395 |
\begin{rawhtml} |
396 |
<!-- CMIREDIR:topo_partial_cells: --> |
397 |
\end{rawhtml} |
398 |
|
399 |
\begin{figure} |
400 |
\begin{center} |
401 |
\resizebox{4.5in}{!}{\includegraphics{part2/vgrid-xz.eps}} |
402 |
\end{center} |
403 |
\caption{ |
404 |
A schematic of the x-r plane showing the location of the |
405 |
non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a |
406 |
tracer cell is given by $h_c(i,j,k) \Delta r_f(k)$ and the physical |
407 |
thickness of the open side is given by $h_w(i,j,k) \Delta r_f(k)$.} |
408 |
\label{fig:hfacs} |
409 |
\end{figure} |
410 |
|
411 |
\cite{adcroft:97} presented two alternatives to the step-wise finite |
412 |
difference representation of topography. The method is known to the |
413 |
engineering community as {\em intersecting boundary method}. It |
414 |
involves allowing the boundary to intersect a grid of cells thereby |
415 |
modifying the shape of those cells intersected. We suggested allowing |
416 |
the topography to take on a piece-wise linear representation (shaved |
417 |
cells) or a simpler piecewise constant representation (partial step). |
418 |
Both show dramatic improvements in solution compared to the |
419 |
traditional full step representation, the piece-wise linear being the |
420 |
best. However, the storage requirements are excessive so the simpler |
421 |
piece-wise constant or partial-step method is all that is currently |
422 |
supported. |
423 |
|
424 |
Fig.~\ref{fig:hfacs} shows a schematic of the x-r plane indicating how |
425 |
the thickness of a level is determined at tracer and u points. |
426 |
\marginpar{$h_c$: {\bf hFacC}} |
427 |
\marginpar{$h_w$: {\bf hFacW}} |
428 |
\marginpar{$h_s$: {\bf hFacS}} |
429 |
The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta |
430 |
r_f(k)$ and the physical thickness of the open side is given by |
431 |
$h_w(i,j,k) \Delta r_f(k)$. Three 3-D descriptors $h_c$, $h_w$ and |
432 |
$h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and |
433 |
{\bf hFacS} respectively. These are calculated in subroutine {\em |
434 |
INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf |
435 |
RECIP\_hFacW} and {\bf RECIP\_hFacS}. |
436 |
|
437 |
The non-dimensional fractions (or h-facs as we call them) are |
438 |
calculated from the model depth array and then processed to avoid tiny |
439 |
volumes. The rule is that if a fraction is less than {\bf hFacMin} |
440 |
then it is rounded to the nearer of $0$ or {\bf hFacMin} or if the |
441 |
physical thickness is less than {\bf hFacMinDr} then it is similarly |
442 |
rounded. The larger of the two methods is used when there is a |
443 |
conflict. By setting {\bf hFacMinDr} equal to or larger than the |
444 |
thinnest nominal layers, $\min{(\Delta z_f)}$, but setting {\bf |
445 |
hFacMin} to some small fraction then the model will only lop thick |
446 |
layers but retain stability based on the thinnest unlopped thickness; |
447 |
$\min{(\Delta z_f,\mbox{\bf hFacMinDr})}$. |
448 |
|
449 |
\fbox{ \begin{minipage}{4.75in} |
450 |
{\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F}) |
451 |
|
452 |
$h_c$: {\bf hFacC} ({\em GRID.h}) |
453 |
|
454 |
$h_w$: {\bf hFacW} ({\em GRID.h}) |
455 |
|
456 |
$h_s$: {\bf hFacS} ({\em GRID.h}) |
457 |
|
458 |
$h_c^{-1}$: {\bf RECIP\_hFacC} ({\em GRID.h}) |
459 |
|
460 |
$h_w^{-1}$: {\bf RECIP\_hFacW} ({\em GRID.h}) |
461 |
|
462 |
$h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRID.h}) |
463 |
|
464 |
\end{minipage} } |
465 |
|
466 |
|
467 |
\section{Continuity and horizontal pressure gradient terms} |
468 |
|
469 |
The core algorithm is based on the ``C grid'' discretization of the |
470 |
continuity equation which can be summarized as: |
471 |
\begin{eqnarray} |
472 |
\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} \\ |
473 |
\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} \\ |
474 |
\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} \\ |
475 |
\delta_i \Delta y_g \Delta r_f h_w u + |
476 |
\delta_j \Delta x_g \Delta r_f h_s v + |
477 |
\delta_k {\cal A}_c w & = & {\cal A}_c \delta_k (P-E)_{r=0} |
478 |
\label{eq:discrete-continuity} |
479 |
\end{eqnarray} |
480 |
where the continuity equation has been most naturally discretized by |
481 |
staggering the three components of velocity as shown in |
482 |
Fig.~\ref{fig:cgrid3d}. The grid lengths $\Delta x_c$ and $\Delta y_c$ |
483 |
are the lengths between tracer points (cell centers). The grid lengths |
484 |
$\Delta x_g$, $\Delta y_g$ are the grid lengths between cell |
485 |
corners. $\Delta r_f$ and $\Delta r_c$ are the distance (in units of |
486 |
$r$) between level interfaces (w-level) and level centers (tracer |
487 |
level). The surface area presented in the vertical is denoted ${\cal |
488 |
A}_c$. The factors $h_w$ and $h_s$ are non-dimensional fractions |
489 |
(between 0 and 1) that represent the fraction cell depth that is |
490 |
``open'' for fluid flow. |
491 |
\marginpar{$h_w$: {\bf hFacW}} |
492 |
\marginpar{$h_s$: {\bf hFacS}} |
493 |
|
494 |
The last equation, the discrete continuity equation, can be summed in |
495 |
the vertical to yield the free-surface equation: |
496 |
\begin{equation} |
497 |
{\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w |
498 |
u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal |
499 |
A}_c(P-E)_{r=0} \label{eq:discrete-freesurface} |
500 |
\end{equation} |
501 |
The source term $P-E$ on the rhs of continuity accounts for the local |
502 |
addition of volume due to excess precipitation and run-off over |
503 |
evaporation and only enters the top-level of the {\em ocean} model. |
504 |
|
505 |
\section{Hydrostatic balance} |
506 |
|
507 |
The vertical momentum equation has the hydrostatic or |
508 |
quasi-hydrostatic balance on the right hand side. This discretization |
509 |
guarantees that the conversion of potential to kinetic energy as |
510 |
derived from the buoyancy equation exactly matches the form derived |
511 |
from the pressure gradient terms when forming the kinetic energy |
512 |
equation. |
513 |
|
514 |
In the ocean, using z-coordinates, the hydrostatic balance terms are |
515 |
discretized: |
516 |
\begin{equation} |
517 |
\epsilon_{nh} \partial_t w |
518 |
+ g \overline{\rho'}^k + \frac{1}{\Delta z} \delta_k \Phi_h' = \ldots |
519 |
\label{eq:discrete_hydro_ocean} |
520 |
\end{equation} |
521 |
|
522 |
In the atmosphere, using p-coordinates, hydrostatic balance is |
523 |
discretized: |
524 |
\begin{equation} |
525 |
\overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0 |
526 |
\label{eq:discrete_hydro_atmos} |
527 |
\end{equation} |
528 |
where $\Delta \Pi$ is the difference in Exner function between the |
529 |
pressure points. The non-hydrostatic equations are not available in |
530 |
the atmosphere. |
531 |
|
532 |
The difference in approach between ocean and atmosphere occurs because |
533 |
of the direct use of the ideal gas equation in forming the potential |
534 |
energy conversion term $\alpha \omega$. The form of these conversion |
535 |
terms is discussed at length in \cite{adcroft:02}. |
536 |
|
537 |
Because of the different representation of hydrostatic balance between |
538 |
ocean and atmosphere there is no elegant way to represent both systems |
539 |
using an arbitrary coordinate. |
540 |
|
541 |
The integration for hydrostatic pressure is made in the positive $r$ |
542 |
direction (increasing k-index). For the ocean, this is from the |
543 |
free-surface down and for the atmosphere this is from the ground up. |
544 |
|
545 |
The calculations are made in the subroutine {\em |
546 |
CALC\_PHI\_HYD}. Inside this routine, one of other of the |
547 |
atmospheric/oceanic form is selected based on the string variable {\bf |
548 |
buoyancyRelation}. |
549 |
|