| 2 |
% $Name$ |
% $Name$ |
| 3 |
|
|
| 4 |
\section{Spatial discretization of the dynamical equations} |
\section{Spatial discretization of the dynamical equations} |
| 5 |
|
\begin{rawhtml} |
| 6 |
|
<!-- CMIREDIR:spatial_discretization_of_dyn_eq: --> |
| 7 |
|
\end{rawhtml} |
| 8 |
|
|
| 9 |
|
Spatial discretization is carried out using the finite volume |
| 10 |
|
method. This amounts to a grid-point method (namely second-order |
| 11 |
|
centered finite difference) in the fluid interior but allows |
| 12 |
|
boundaries to intersect a regular grid allowing a more accurate |
| 13 |
|
representation of the position of the boundary. We treat the |
| 14 |
|
horizontal and vertical directions as separable and differently. |
| 15 |
|
|
| 16 |
|
\input{part2/notation} |
| 17 |
|
|
| 18 |
|
|
| 19 |
|
\subsection{The finite volume method: finite volumes versus finite difference} |
| 20 |
|
\begin{rawhtml} |
| 21 |
|
<!-- CMIREDIR:finite_volume: --> |
| 22 |
|
\end{rawhtml} |
| 23 |
|
|
| 24 |
|
|
| 25 |
|
|
| 26 |
|
The finite volume method is used to discretize the equations in |
| 27 |
|
space. The expression ``finite volume'' actually has two meanings; one |
| 28 |
|
is the method of embedded or intersecting boundaries (shaved or lopped |
| 29 |
|
cells in our terminology) and the other is non-linear interpolation |
| 30 |
|
methods that can deal with non-smooth solutions such as shocks |
| 31 |
|
(i.e. flux limiters for advection). Both make use of the integral form |
| 32 |
|
of the conservation laws to which the {\it weak solution} is a |
| 33 |
|
solution on each finite volume of (sub-domain). The weak solution can |
| 34 |
|
be constructed out of piece-wise constant elements or be |
| 35 |
|
differentiable. The differentiable equations can not be satisfied by |
| 36 |
|
piece-wise constant functions. |
| 37 |
|
|
| 38 |
|
As an example, the 1-D constant coefficient advection-diffusion |
| 39 |
|
equation: |
| 40 |
|
\begin{displaymath} |
| 41 |
|
\partial_t \theta + \partial_x ( u \theta - \kappa \partial_x \theta ) = 0 |
| 42 |
|
\end{displaymath} |
| 43 |
|
can be discretized by integrating over finite sub-domains, i.e. |
| 44 |
|
the lengths $\Delta x_i$: |
| 45 |
|
\begin{displaymath} |
| 46 |
|
\Delta x \partial_t \theta + \delta_i ( F ) = 0 |
| 47 |
|
\end{displaymath} |
| 48 |
|
is exact if $\theta(x)$ is piece-wise constant over the interval |
| 49 |
|
$\Delta x_i$ or more generally if $\theta_i$ is defined as the average |
| 50 |
|
over the interval $\Delta x_i$. |
| 51 |
|
|
| 52 |
|
The flux, $F_{i-1/2}$, must be approximated: |
| 53 |
|
\begin{displaymath} |
| 54 |
|
F = u \overline{\theta} - \frac{\kappa}{\Delta x_c} \partial_i \theta |
| 55 |
|
\end{displaymath} |
| 56 |
|
and this is where truncation errors can enter the solution. The |
| 57 |
|
method for obtaining $\overline{\theta}$ is unspecified and a wide |
| 58 |
|
range of possibilities exist including centered and upwind |
| 59 |
|
interpolation, polynomial fits based on the the volume average |
| 60 |
|
definitions of quantities and non-linear interpolation such as |
| 61 |
|
flux-limiters. |
| 62 |
|
|
| 63 |
|
Choosing simple centered second-order interpolation and differencing |
| 64 |
|
recovers the same ODE's resulting from finite differencing for the |
| 65 |
|
interior of a fluid. Differences arise at boundaries where a boundary |
| 66 |
|
is not positioned on a regular or smoothly varying grid. This method |
| 67 |
|
is used to represent the topography using lopped cell, see |
| 68 |
|
\cite{adcroft:97}. Subtle difference also appear in more than one |
| 69 |
|
dimension away from boundaries. This happens because the each |
| 70 |
|
direction is discretized independently in the finite difference method |
| 71 |
|
while the integrating over finite volume implicitly treats all |
| 72 |
|
directions simultaneously. Illustration of this is given in |
| 73 |
|
\cite{ac:02}. |
| 74 |
|
|
| 75 |
\subsection{C grid staggering of variables} |
\subsection{C grid staggering of variables} |
| 76 |
|
|
| 77 |
\begin{figure} |
\begin{figure} |
| 78 |
\centerline{ \resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} } |
\begin{center} |
| 79 |
\label{fig-cgrid3d} |
\resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} |
| 80 |
|
\end{center} |
| 81 |
\caption{Three dimensional staggering of velocity components. This |
\caption{Three dimensional staggering of velocity components. This |
| 82 |
facilitates the natural discretization of the continuity and tracer |
facilitates the natural discretization of the continuity and tracer |
| 83 |
equations. } |
equations. } |
| 84 |
|
\label{fig:cgrid3d} |
| 85 |
\end{figure} |
\end{figure} |
| 86 |
|
|
| 87 |
|
The basic algorithm employed for stepping forward the momentum |
| 88 |
|
equations is based on retaining non-divergence of the flow at all |
| 89 |
|
times. This is most naturally done if the components of flow are |
| 90 |
|
staggered in space in the form of an Arakawa C grid \cite{arakawa:77}. |
| 91 |
|
|
| 92 |
|
Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$) |
| 93 |
|
staggered in space such that the zonal component falls on the |
| 94 |
|
interface between continuity cells in the zonal direction. Similarly |
| 95 |
|
for the meridional and vertical directions. The continuity cell is |
| 96 |
|
synonymous with tracer cells (they are one and the same). |
| 97 |
|
|
| 98 |
|
|
| 99 |
|
|
| 100 |
|
\subsection{Grid initialization and data} |
| 101 |
|
|
| 102 |
|
Initialization of grid data is controlled by subroutine {\em |
| 103 |
|
INI\_GRID} which in calls {\em INI\_VERTICAL\_GRID} to initialize the |
| 104 |
|
vertical grid, and then either of {\em INI\_CARTESIAN\_GRID}, {\em |
| 105 |
|
INI\_SPHERICAL\_POLAR\_GRID} or {\em INI\_CURV\-ILINEAR\_GRID} to |
| 106 |
|
initialize the horizontal grid for cartesian, spherical-polar or |
| 107 |
|
curvilinear coordinates respectively. |
| 108 |
|
|
| 109 |
|
The reciprocals of all grid quantities are pre-calculated and this is |
| 110 |
|
done in subroutine {\em INI\_MASKS\_ETC} which is called later by |
| 111 |
|
subroutine {\em INITIALIZE\_FIXED}. |
| 112 |
|
|
| 113 |
|
All grid descriptors are global arrays and stored in common blocks in |
| 114 |
|
{\em GRID.h} and a generally declared as {\em \_RS}. |
| 115 |
|
|
| 116 |
|
\fbox{ \begin{minipage}{4.75in} |
| 117 |
|
{\em S/R INI\_GRID} ({\em model/src/ini\_grid.F}) |
| 118 |
|
|
| 119 |
|
{\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F}) |
| 120 |
|
|
| 121 |
|
grid data: ({\em model/inc/GRID.h}) |
| 122 |
|
\end{minipage} } |
| 123 |
|
|
| 124 |
|
|
| 125 |
\subsection{Horizontal grid} |
\subsection{Horizontal grid} |
| 126 |
|
\label{sec:spatial_discrete_horizontal_grid} |
| 127 |
|
|
| 128 |
\begin{figure} |
\begin{figure} |
| 129 |
\centerline{ \begin{tabular}{cc} |
\begin{center} |
| 130 |
\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}} |
\begin{tabular}{cc} |
| 131 |
& \resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}} |
\raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}} |
| 132 |
|
& \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}} |
| 133 |
\\ |
\\ |
| 134 |
\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}} |
\raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}} |
| 135 |
& \resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}} |
& \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}} |
| 136 |
\end{tabular} } |
\end{tabular} |
| 137 |
\label{fig-hgrid} |
\end{center} |
| 138 |
\caption{Three dimensional staggering of velocity components. This |
\caption{ |
| 139 |
facilitates the natural discretization of the continuity and tracer |
Staggering of horizontal grid descriptors (lengths and areas). The |
| 140 |
equations. } |
grid lines indicate the tracer cell boundaries and are the reference |
| 141 |
|
grid for all panels. a) The area of a tracer cell, $A_c$, is bordered |
| 142 |
|
by the lengths $\Delta x_g$ and $\Delta y_g$. b) The area of a |
| 143 |
|
vorticity cell, $A_\zeta$, is bordered by the lengths $\Delta x_c$ and |
| 144 |
|
$\Delta y_c$. c) The area of a u cell, $A_c$, is bordered by the |
| 145 |
|
lengths $\Delta x_v$ and $\Delta y_f$. d) The area of a v cell, $A_c$, |
| 146 |
|
is bordered by the lengths $\Delta x_f$ and $\Delta y_u$.} |
| 147 |
|
\label{fig:hgrid} |
| 148 |
\end{figure} |
\end{figure} |
| 149 |
|
|
| 150 |
|
The model domain is decomposed into tiles and within each tile a |
| 151 |
|
quasi-regular grid is used. A tile is the basic unit of domain |
| 152 |
|
decomposition for parallelization but may be used whether parallelized |
| 153 |
|
or not; see section \ref{sect:domain_decomposition} for more details. |
| 154 |
|
Although the tiles may be patched together in an unstructured manner |
| 155 |
|
(i.e. irregular or non-tessilating pattern), the interior of tiles is |
| 156 |
|
a structured grid of quadrilateral cells. The horizontal coordinate |
| 157 |
|
system is orthogonal curvilinear meaning we can not necessarily treat |
| 158 |
|
the two horizontal directions as separable. Instead, each cell in the |
| 159 |
|
horizontal grid is described by the length of it's sides and it's |
| 160 |
|
area. |
| 161 |
|
|
| 162 |
|
The grid information is quite general and describes any of the |
| 163 |
|
available coordinates systems, cartesian, spherical-polar or |
| 164 |
|
curvilinear. All that is necessary to distinguish between the |
| 165 |
|
coordinate systems is to initialize the grid data (descriptors) |
| 166 |
|
appropriately. |
| 167 |
|
|
| 168 |
|
In the following, we refer to the orientation of quantities on the |
| 169 |
|
computational grid using geographic terminology such as points of the |
| 170 |
|
compass. |
| 171 |
|
\marginpar{Caution!} |
| 172 |
|
This is purely for convenience but should note be confused |
| 173 |
|
with the actual geographic orientation of model quantities. |
| 174 |
|
|
| 175 |
|
Fig.~\ref{fig:hgrid}a shows the tracer cell (synonymous with the |
| 176 |
|
continuity cell). The length of the southern edge, $\Delta x_g$, |
| 177 |
|
western edge, $\Delta y_g$ and surface area, $A_c$, presented in the |
| 178 |
|
vertical are stored in arrays {\bf DXg}, {\bf DYg} and {\bf rAc}. |
| 179 |
|
\marginpar{$A_c$: {\bf rAc}} |
| 180 |
|
\marginpar{$\Delta x_g$: {\bf DXg}} |
| 181 |
|
\marginpar{$\Delta y_g$: {\bf DYg}} |
| 182 |
|
The ``g'' suffix indicates that the lengths are along the defining |
| 183 |
|
grid boundaries. The ``c'' suffix associates the quantity with the |
| 184 |
|
cell centers. The quantities are staggered in space and the indexing |
| 185 |
|
is such that {\bf DXg(i,j)} is positioned to the south of {\bf |
| 186 |
|
rAc(i,j)} and {\bf DYg(i,j)} positioned to the west. |
| 187 |
|
|
| 188 |
|
Fig.~\ref{fig:hgrid}b shows the vorticity cell. The length of the |
| 189 |
|
southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface |
| 190 |
|
area, $A_\zeta$, presented in the vertical are stored in arrays {\bf |
| 191 |
|
DXg}, {\bf DYg} and {\bf rAz}. |
| 192 |
|
\marginpar{$A_\zeta$: {\bf rAz}} |
| 193 |
|
\marginpar{$\Delta x_c$: {\bf DXc}} |
| 194 |
|
\marginpar{$\Delta y_c$: {\bf DYc}} |
| 195 |
|
The ``z'' suffix indicates that the lengths are measured between the |
| 196 |
|
cell centers and the ``$\zeta$'' suffix associates points with the |
| 197 |
|
vorticity points. The quantities are staggered in space and the |
| 198 |
|
indexing is such that {\bf DXc(i,j)} is positioned to the north of |
| 199 |
|
{\bf rAc(i,j)} and {\bf DYc(i,j)} positioned to the east. |
| 200 |
|
|
| 201 |
|
Fig.~\ref{fig:hgrid}c shows the ``u'' or western (w) cell. The length of |
| 202 |
|
the southern edge, $\Delta x_v$, eastern edge, $\Delta y_f$ and |
| 203 |
|
surface area, $A_w$, presented in the vertical are stored in arrays |
| 204 |
|
{\bf DXv}, {\bf DYf} and {\bf rAw}. |
| 205 |
|
\marginpar{$A_w$: {\bf rAw}} |
| 206 |
|
\marginpar{$\Delta x_v$: {\bf DXv}} |
| 207 |
|
\marginpar{$\Delta y_f$: {\bf DYf}} |
| 208 |
|
The ``v'' suffix indicates that the length is measured between the |
| 209 |
|
v-points, the ``f'' suffix indicates that the length is measured |
| 210 |
|
between the (tracer) cell faces and the ``w'' suffix associates points |
| 211 |
|
with the u-points (w stands for west). The quantities are staggered in |
| 212 |
|
space and the indexing is such that {\bf DXv(i,j)} is positioned to |
| 213 |
|
the south of {\bf rAw(i,j)} and {\bf DYf(i,j)} positioned to the east. |
| 214 |
|
|
| 215 |
|
Fig.~\ref{fig:hgrid}d shows the ``v'' or southern (s) cell. The length of |
| 216 |
|
the northern edge, $\Delta x_f$, western edge, $\Delta y_u$ and |
| 217 |
|
surface area, $A_s$, presented in the vertical are stored in arrays |
| 218 |
|
{\bf DXf}, {\bf DYu} and {\bf rAs}. |
| 219 |
|
\marginpar{$A_s$: {\bf rAs}} |
| 220 |
|
\marginpar{$\Delta x_f$: {\bf DXf}} |
| 221 |
|
\marginpar{$\Delta y_u$: {\bf DYu}} |
| 222 |
|
The ``u'' suffix indicates that the length is measured between the |
| 223 |
|
u-points, the ``f'' suffix indicates that the length is measured |
| 224 |
|
between the (tracer) cell faces and the ``s'' suffix associates points |
| 225 |
|
with the v-points (s stands for south). The quantities are staggered |
| 226 |
|
in space and the indexing is such that {\bf DXf(i,j)} is positioned to |
| 227 |
|
the north of {\bf rAs(i,j)} and {\bf DYu(i,j)} positioned to the west. |
| 228 |
|
|
| 229 |
|
\fbox{ \begin{minipage}{4.75in} |
| 230 |
|
{\em S/R INI\_CARTESIAN\_GRID} ({\em |
| 231 |
|
model/src/ini\_cartesian\_grid.F}) |
| 232 |
|
|
| 233 |
|
{\em S/R INI\_SPHERICAL\_POLAR\_GRID} ({\em |
| 234 |
|
model/src/ini\_spherical\_polar\_grid.F}) |
| 235 |
|
|
| 236 |
|
{\em S/R INI\_CURVILINEAR\_GRID} ({\em |
| 237 |
|
model/src/ini\_curvilinear\_grid.F}) |
| 238 |
|
|
| 239 |
|
$A_c$, $A_\zeta$, $A_w$, $A_s$: {\bf rAc}, {\bf rAz}, {\bf rAw}, {\bf rAs} |
| 240 |
|
({\em GRID.h}) |
| 241 |
|
|
| 242 |
|
$\Delta x_g$, $\Delta y_g$: {\bf DXg}, {\bf DYg} ({\em GRID.h}) |
| 243 |
|
|
| 244 |
|
$\Delta x_c$, $\Delta y_c$: {\bf DXc}, {\bf DYc} ({\em GRID.h}) |
| 245 |
|
|
| 246 |
|
$\Delta x_f$, $\Delta y_f$: {\bf DXf}, {\bf DYf} ({\em GRID.h}) |
| 247 |
|
|
| 248 |
|
$\Delta x_v$, $\Delta y_u$: {\bf DXv}, {\bf DYu} ({\em GRID.h}) |
| 249 |
|
|
| 250 |
|
\end{minipage} } |
| 251 |
|
|
| 252 |
|
\subsubsection{Reciprocals of horizontal grid descriptors} |
| 253 |
|
|
| 254 |
|
%\marginpar{$A_c^{-1}$: {\bf RECIP\_rAc}} |
| 255 |
|
%\marginpar{$A_\zeta^{-1}$: {\bf RECIP\_rAz}} |
| 256 |
|
%\marginpar{$A_w^{-1}$: {\bf RECIP\_rAw}} |
| 257 |
|
%\marginpar{$A_s^{-1}$: {\bf RECIP\_rAs}} |
| 258 |
|
Lengths and areas appear in the denominator of expressions as much as |
| 259 |
|
in the numerator. For efficiency and portability, we pre-calculate the |
| 260 |
|
reciprocal of the horizontal grid quantities so that in-line divisions |
| 261 |
|
can be avoided. |
| 262 |
|
|
| 263 |
|
%\marginpar{$\Delta x_g^{-1}$: {\bf RECIP\_DXg}} |
| 264 |
|
%\marginpar{$\Delta y_g^{-1}$: {\bf RECIP\_DYg}} |
| 265 |
|
%\marginpar{$\Delta x_c^{-1}$: {\bf RECIP\_DXc}} |
| 266 |
|
%\marginpar{$\Delta y_c^{-1}$: {\bf RECIP\_DYc}} |
| 267 |
|
%\marginpar{$\Delta x_f^{-1}$: {\bf RECIP\_DXf}} |
| 268 |
|
%\marginpar{$\Delta y_f^{-1}$: {\bf RECIP\_DYf}} |
| 269 |
|
%\marginpar{$\Delta x_v^{-1}$: {\bf RECIP\_DXv}} |
| 270 |
|
%\marginpar{$\Delta y_u^{-1}$: {\bf RECIP\_DYu}} |
| 271 |
|
For each grid descriptor (array) there is a reciprocal named using the |
| 272 |
|
prefix {\bf RECIP\_}. This doubles the amount of storage in {\em |
| 273 |
|
GRID.h} but they are all only 2-D descriptors. |
| 274 |
|
|
| 275 |
|
\fbox{ \begin{minipage}{4.75in} |
| 276 |
|
{\em S/R INI\_MASKS\_ETC} ({\em |
| 277 |
|
model/src/ini\_masks\_etc.F}) |
| 278 |
|
|
| 279 |
|
$A_c^{-1}$: {\bf RECIP\_Ac} ({\em GRID.h}) |
| 280 |
|
|
| 281 |
|
$A_\zeta^{-1}$: {\bf RECIP\_Az} ({\em GRID.h}) |
| 282 |
|
|
| 283 |
|
$A_w^{-1}$: {\bf RECIP\_Aw} ({\em GRID.h}) |
| 284 |
|
|
| 285 |
|
$A_s^{-1}$: {\bf RECIP\_As} ({\em GRID.h}) |
| 286 |
|
|
| 287 |
|
$\Delta x_g^{-1}$, $\Delta y_g^{-1}$: {\bf RECIP\_DXg}, {\bf RECIP\_DYg} ({\em GRID.h}) |
| 288 |
|
|
| 289 |
|
$\Delta x_c^{-1}$, $\Delta y_c^{-1}$: {\bf RECIP\_DXc}, {\bf RECIP\_DYc} ({\em GRID.h}) |
| 290 |
|
|
| 291 |
|
$\Delta x_f^{-1}$, $\Delta y_f^{-1}$: {\bf RECIP\_DXf}, {\bf RECIP\_DYf} ({\em GRID.h}) |
| 292 |
|
|
| 293 |
|
$\Delta x_v^{-1}$, $\Delta y_u^{-1}$: {\bf RECIP\_DXv}, {\bf RECIP\_DYu} ({\em GRID.h}) |
| 294 |
|
|
| 295 |
|
\end{minipage} } |
| 296 |
|
|
| 297 |
|
\subsubsection{Cartesian coordinates} |
| 298 |
|
|
| 299 |
|
Cartesian coordinates are selected when the logical flag {\bf |
| 300 |
|
using\-Cartes\-ianGrid} in namelist {\em PARM04} is set to true. The grid |
| 301 |
|
spacing can be set to uniform via scalars {\bf dXspacing} and {\bf |
| 302 |
|
dYspacing} in namelist {\em PARM04} or to variable resolution by the |
| 303 |
|
vectors {\bf DELX} and {\bf DELY}. Units are normally |
| 304 |
|
meters. Non-dimensional coordinates can be used by interpreting the |
| 305 |
|
gravitational constant as the Rayleigh number. |
| 306 |
|
|
| 307 |
|
\subsubsection{Spherical-polar coordinates} |
| 308 |
|
|
| 309 |
|
Spherical coordinates are selected when the logical flag {\bf |
| 310 |
|
using\-Spherical\-PolarGrid} in namelist {\em PARM04} is set to true. The |
| 311 |
|
grid spacing can be set to uniform via scalars {\bf dXspacing} and |
| 312 |
|
{\bf dYspacing} in namelist {\em PARM04} or to variable resolution by |
| 313 |
|
the vectors {\bf DELX} and {\bf DELY}. Units of these namelist |
| 314 |
|
variables are alway degrees. The horizontal grid descriptors are |
| 315 |
|
calculated from these namelist variables have units of meters. |
| 316 |
|
|
| 317 |
|
\subsubsection{Curvilinear coordinates} |
| 318 |
|
|
| 319 |
|
Curvilinear coordinates are selected when the logical flag {\bf |
| 320 |
|
using\-Curvil\-inear\-Grid} in namelist {\em PARM04} is set to true. The |
| 321 |
|
grid spacing can not be set via the namelist. Instead, the grid |
| 322 |
|
descriptors are read from data files, one for each descriptor. As for |
| 323 |
|
other grids, the horizontal grid descriptors have units of meters. |
| 324 |
|
|
| 325 |
|
|
| 326 |
\subsection{Vertical grid} |
\subsection{Vertical grid} |
| 327 |
|
|
| 328 |
\begin{figure} |
\begin{figure} |
| 329 |
\centerline{ \begin{tabular}{cc} |
\begin{center} |
| 330 |
\raisebox{4in}{a)} |
\begin{tabular}{cc} |
| 331 |
\resizebox{!}{4in}{ \includegraphics{part2/vgrid-cellcentered.eps}} |
\raisebox{4in}{a)} \resizebox{!}{4in}{ |
| 332 |
& |
\includegraphics{part2/vgrid-cellcentered.eps}} & \raisebox{4in}{b)} |
| 333 |
\raisebox{4in}{b)} |
\resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}} |
| 334 |
\resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}} |
\end{tabular} |
| 335 |
\end{tabular} } |
\end{center} |
|
\label{fig-vgrid} |
|
| 336 |
\caption{Two versions of the vertical grid. a) The cell centered |
\caption{Two versions of the vertical grid. a) The cell centered |
| 337 |
approach where the interface depths are specified and the tracer |
approach where the interface depths are specified and the tracer |
| 338 |
points centered in between the interfaces. b) The interface centered |
points centered in between the interfaces. b) The interface centered |
| 339 |
approach where tracer levels are specified and the w-interfaces are |
approach where tracer levels are specified and the w-interfaces are |
| 340 |
centered in between.} |
centered in between.} |
| 341 |
|
\label{fig:vgrid} |
| 342 |
|
\end{figure} |
| 343 |
|
|
| 344 |
|
As for the horizontal grid, we use the suffixes ``c'' and ``f'' to |
| 345 |
|
indicates faces and centers. Fig.~\ref{fig:vgrid}a shows the default |
| 346 |
|
vertical grid used by the model. |
| 347 |
|
\marginpar{$\Delta r_f$: {\bf DRf}} |
| 348 |
|
\marginpar{$\Delta r_c$: {\bf DRc}} |
| 349 |
|
$\Delta r_f$ is the difference in $r$ |
| 350 |
|
(vertical coordinate) between the faces (i.e. $\Delta r_f \equiv - |
| 351 |
|
\delta_k r$ where the minus sign appears due to the convention that the |
| 352 |
|
surface layer has index $k=1$.). |
| 353 |
|
|
| 354 |
|
The vertical grid is calculated in subroutine {\em |
| 355 |
|
INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in |
| 356 |
|
namelist {\em PARM04}. The units of ``r'' are either meters or Pascals |
| 357 |
|
depending on the isomorphism being used which in turn is dependent |
| 358 |
|
only on the choice of equation of state. |
| 359 |
|
|
| 360 |
|
There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which |
| 361 |
|
dictate whether z- or |
| 362 |
|
\marginpar{Caution!} |
| 363 |
|
p- coordinates are to be used but we intend to |
| 364 |
|
phase this out since they are redundant. |
| 365 |
|
|
| 366 |
|
The reciprocals $\Delta r_f^{-1}$ and $\Delta r_c^{-1}$ are |
| 367 |
|
pre-calculated (also in subroutine {\em INI\_VERTICAL\_GRID}). All |
| 368 |
|
vertical grid descriptors are stored in common blocks in {\em GRID.h}. |
| 369 |
|
|
| 370 |
|
The above grid (Fig.~\ref{fig:vgrid}a) is known as the cell centered |
| 371 |
|
approach because the tracer points are at cell centers; the cell |
| 372 |
|
centers are mid-way between the cell interfaces. |
| 373 |
|
This discretization is selected when the thickness of the |
| 374 |
|
levels are provided ({\bf delR}, parameter file {\em data}, |
| 375 |
|
namelist {\em PARM04}) |
| 376 |
|
An alternative, the vertex or interface centered approach, is shown in |
| 377 |
|
Fig.~\ref{fig:vgrid}b. Here, the interior interfaces are positioned |
| 378 |
|
mid-way between the tracer nodes (no longer cell centers). This |
| 379 |
|
approach is formally more accurate for evaluation of hydrostatic |
| 380 |
|
pressure and vertical advection but historically the cell centered |
| 381 |
|
approach has been used. An alternative form of subroutine {\em |
| 382 |
|
INI\_VERTICAL\_GRID} is used to select the interface centered approach |
| 383 |
|
This form requires to specify $Nr+1$ vertical distances {\bf delRc} |
| 384 |
|
(parameter file {\em data}, namelist {\em PARM04}, e.g. |
| 385 |
|
{\em verification/ideal\_2D\_oce/input/data}) |
| 386 |
|
corresponding to surface to center, $Nr-1$ center to center, and center to |
| 387 |
|
bottom distances. |
| 388 |
|
%but no run time option is currently available. |
| 389 |
|
|
| 390 |
|
\fbox{ \begin{minipage}{4.75in} |
| 391 |
|
{\em S/R INI\_VERTICAL\_GRID} ({\em |
| 392 |
|
model/src/ini\_vertical\_grid.F}) |
| 393 |
|
|
| 394 |
|
$\Delta r_f$: {\bf DRf} ({\em GRID.h}) |
| 395 |
|
|
| 396 |
|
$\Delta r_c$: {\bf DRc} ({\em GRID.h}) |
| 397 |
|
|
| 398 |
|
$\Delta r_f^{-1}$: {\bf RECIP\_DRf} ({\em GRID.h}) |
| 399 |
|
|
| 400 |
|
$\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\em GRID.h}) |
| 401 |
|
|
| 402 |
|
\end{minipage} } |
| 403 |
|
|
| 404 |
|
|
| 405 |
|
\subsection{Topography: partially filled cells} |
| 406 |
|
\begin{rawhtml} |
| 407 |
|
<!-- CMIREDIR:topo_partial_cells: --> |
| 408 |
|
\end{rawhtml} |
| 409 |
|
|
| 410 |
|
\begin{figure} |
| 411 |
|
\begin{center} |
| 412 |
|
\resizebox{4.5in}{!}{\includegraphics{part2/vgrid-xz.eps}} |
| 413 |
|
\end{center} |
| 414 |
|
\caption{ |
| 415 |
|
A schematic of the x-r plane showing the location of the |
| 416 |
|
non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a |
| 417 |
|
tracer cell is given by $h_c(i,j,k) \Delta r_f(k)$ and the physical |
| 418 |
|
thickness of the open side is given by $h_w(i,j,k) \Delta r_f(k)$.} |
| 419 |
|
\label{fig:hfacs} |
| 420 |
\end{figure} |
\end{figure} |
| 421 |
|
|
| 422 |
|
\cite{adcroft:97} presented two alternatives to the step-wise finite |
| 423 |
|
difference representation of topography. The method is known to the |
| 424 |
|
engineering community as {\em intersecting boundary method}. It |
| 425 |
|
involves allowing the boundary to intersect a grid of cells thereby |
| 426 |
|
modifying the shape of those cells intersected. We suggested allowing |
| 427 |
|
the topography to take on a piece-wise linear representation (shaved |
| 428 |
|
cells) or a simpler piecewise constant representation (partial step). |
| 429 |
|
Both show dramatic improvements in solution compared to the |
| 430 |
|
traditional full step representation, the piece-wise linear being the |
| 431 |
|
best. However, the storage requirements are excessive so the simpler |
| 432 |
|
piece-wise constant or partial-step method is all that is currently |
| 433 |
|
supported. |
| 434 |
|
|
| 435 |
|
Fig.~\ref{fig:hfacs} shows a schematic of the x-r plane indicating how |
| 436 |
|
the thickness of a level is determined at tracer and u points. |
| 437 |
|
\marginpar{$h_c$: {\bf hFacC}} |
| 438 |
|
\marginpar{$h_w$: {\bf hFacW}} |
| 439 |
|
\marginpar{$h_s$: {\bf hFacS}} |
| 440 |
|
The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta |
| 441 |
|
r_f(k)$ and the physical thickness of the open side is given by |
| 442 |
|
$h_w(i,j,k) \Delta r_f(k)$. Three 3-D descriptors $h_c$, $h_w$ and |
| 443 |
|
$h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and |
| 444 |
|
{\bf hFacS} respectively. These are calculated in subroutine {\em |
| 445 |
|
INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf |
| 446 |
|
RECIP\_hFacW} and {\bf RECIP\_hFacS}. |
| 447 |
|
|
| 448 |
|
The non-dimensional fractions (or h-facs as we call them) are |
| 449 |
|
calculated from the model depth array and then processed to avoid tiny |
| 450 |
|
volumes. The rule is that if a fraction is less than {\bf hFacMin} |
| 451 |
|
then it is rounded to the nearer of $0$ or {\bf hFacMin} or if the |
| 452 |
|
physical thickness is less than {\bf hFacMinDr} then it is similarly |
| 453 |
|
rounded. The larger of the two methods is used when there is a |
| 454 |
|
conflict. By setting {\bf hFacMinDr} equal to or larger than the |
| 455 |
|
thinnest nominal layers, $\min{(\Delta z_f)}$, but setting {\bf |
| 456 |
|
hFacMin} to some small fraction then the model will only lop thick |
| 457 |
|
layers but retain stability based on the thinnest unlopped thickness; |
| 458 |
|
$\min{(\Delta z_f,\mbox{\bf hFacMinDr})}$. |
| 459 |
|
|
| 460 |
|
\fbox{ \begin{minipage}{4.75in} |
| 461 |
|
{\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F}) |
| 462 |
|
|
| 463 |
|
$h_c$: {\bf hFacC} ({\em GRID.h}) |
| 464 |
|
|
| 465 |
|
$h_w$: {\bf hFacW} ({\em GRID.h}) |
| 466 |
|
|
| 467 |
|
$h_s$: {\bf hFacS} ({\em GRID.h}) |
| 468 |
|
|
| 469 |
|
$h_c^{-1}$: {\bf RECIP\_hFacC} ({\em GRID.h}) |
| 470 |
|
|
| 471 |
|
$h_w^{-1}$: {\bf RECIP\_hFacW} ({\em GRID.h}) |
| 472 |
|
|
| 473 |
|
$h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRID.h}) |
| 474 |
|
|
| 475 |
|
\end{minipage} } |
| 476 |
|
|
| 477 |
|
|
| 478 |
|
\section{Continuity and horizontal pressure gradient terms} |
| 479 |
|
\begin{rawhtml} |
| 480 |
|
<!-- CMIREDIR:continuity_and_horizontal_pressure: --> |
| 481 |
|
\end{rawhtml} |
| 482 |
|
|
|
\subsection{Continuity and horizontal pressure gradient terms} |
|
| 483 |
|
|
| 484 |
The core algorithm is based on the ``C grid'' discretization of the |
The core algorithm is based on the ``C grid'' discretization of the |
| 485 |
continuity equation which can be summarized as: |
continuity equation which can be summarized as: |
| 486 |
\begin{eqnarray} |
\begin{eqnarray} |
| 487 |
\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} \\ |
| 488 |
\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} \\ |
| 489 |
\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} \\ |
| 490 |
\delta_i \Delta y_g \Delta r_f h_w u + |
\delta_i \Delta y_g \Delta r_f h_w u + |
| 491 |
\delta_j \Delta x_g \Delta r_f h_s v + |
\delta_j \Delta x_g \Delta r_f h_s v + |
| 492 |
\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} |
| 493 |
|
\label{eq:discrete-continuity} |
| 494 |
\end{eqnarray} |
\end{eqnarray} |
| 495 |
where the continuity equation has been most naturally discretized by |
where the continuity equation has been most naturally discretized by |
| 496 |
staggering the three components of velocity as shown in |
staggering the three components of velocity as shown in |
| 497 |
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$ |
| 498 |
are the lengths between tracer points (cell centers). The grid lengths |
are the lengths between tracer points (cell centers). The grid lengths |
| 499 |
$\Delta x_g$, $\Delta y_g$ are the grid lengths between cell |
$\Delta x_g$, $\Delta y_g$ are the grid lengths between cell |
| 500 |
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 |
| 507 |
\marginpar{$h_s$: {\bf hFacS}} |
\marginpar{$h_s$: {\bf hFacS}} |
| 508 |
|
|
| 509 |
The last equation, the discrete continuity equation, can be summed in |
The last equation, the discrete continuity equation, can be summed in |
| 510 |
the vertical to yeild the free-surface equation: |
the vertical to yield the free-surface equation: |
| 511 |
\begin{equation} |
\begin{equation} |
| 512 |
{\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 |
| 513 |
{\cal A}_c(P-E)_{r=0} |
u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal |
| 514 |
|
A}_c(P-E)_{r=0} \label{eq:discrete-freesurface} |
| 515 |
\end{equation} |
\end{equation} |
| 516 |
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 |
| 517 |
addition of volume due to excess precipitation and run-off over |
addition of volume due to excess precipitation and run-off over |
| 518 |
evaporation and only enters the top-level of the {\em ocean} model. |
evaporation and only enters the top-level of the {\em ocean} model. |
| 519 |
|
|
| 520 |
\subsection{Hydrostatic balance} |
\section{Hydrostatic balance} |
| 521 |
|
\begin{rawhtml} |
| 522 |
|
<!-- CMIREDIR:hydrostatic_balance: --> |
| 523 |
|
\end{rawhtml} |
| 524 |
|
|
| 525 |
The vertical momentum equation has the hydrostatic or |
The vertical momentum equation has the hydrostatic or |
| 526 |
quasi-hydrostatic balance on the right hand side. This discretization |
quasi-hydrostatic balance on the right hand side. This discretization |
| 529 |
from the pressure gradient terms when forming the kinetic energy |
from the pressure gradient terms when forming the kinetic energy |
| 530 |
equation. |
equation. |
| 531 |
|
|
| 532 |
In the ocean, using z-ccordinates, the hydrostatic balance terms are |
In the ocean, using z-coordinates, the hydrostatic balance terms are |
| 533 |
discretized: |
discretized: |
| 534 |
\begin{equation} |
\begin{equation} |
| 535 |
\epsilon_{nh} \partial_t w |
\epsilon_{nh} \partial_t w |
| 536 |
+ 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 |
| 537 |
|
\label{eq:discrete_hydro_ocean} |
| 538 |
\end{equation} |
\end{equation} |
| 539 |
|
|
| 540 |
In the atmosphere, using p-coordinates, hydrostatic balance is |
In the atmosphere, using p-coordinates, hydrostatic balance is |
| 541 |
discretized: |
discretized: |
| 542 |
\begin{equation} |
\begin{equation} |
| 543 |
\overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0 |
\overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0 |
| 544 |
|
\label{eq:discrete_hydro_atmos} |
| 545 |
\end{equation} |
\end{equation} |
| 546 |
where $\Delta \Pi$ is the difference in Exner function between the |
where $\Delta \Pi$ is the difference in Exner function between the |
| 547 |
pressure points. The non-hydrostatic equations are not available in |
pressure points. The non-hydrostatic equations are not available in |
| 549 |
|
|
| 550 |
The difference in approach between ocean and atmosphere occurs because |
The difference in approach between ocean and atmosphere occurs because |
| 551 |
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 |
| 552 |
energy conversion term $\alpha \omega$. The form of these consversion |
energy conversion term $\alpha \omega$. The form of these conversion |
| 553 |
terms is discussed at length in \cite{Adcroft01}. |
terms is discussed at length in \cite{adcroft:02}. |
| 554 |
|
|
| 555 |
Because of the different representation of hydrostatic balance between |
Because of the different representation of hydrostatic balance between |
| 556 |
ocean and atmosphere there is no elegant way to represent both systems |
ocean and atmosphere there is no elegant way to represent both systems |
| 565 |
atmospheric/oceanic form is selected based on the string variable {\bf |
atmospheric/oceanic form is selected based on the string variable {\bf |
| 566 |
buoyancyRelation}. |
buoyancyRelation}. |
| 567 |
|
|
|
\subsection{Flux-form momentum equations} |
|
|
|
|
|
The original finite volume model was based on the Eulerian flux form |
|
|
momentum equations. This is the default though the vector invariant |
|
|
form is optionally available (and recommended in some cases). |
|
|
|
|
|
The ``G's'' (our colloquial name for all terms on rhs!) are broken |
|
|
into the various advective, Coriolis, horizontal dissipation, vertical |
|
|
dissipation and metric forces: |
|
|
\marginpar{$G_u$: {\bf Gu} } |
|
|
\marginpar{$G_v$: {\bf Gv} } |
|
|
\marginpar{$G_w$: {\bf Gw} } |
|
|
\begin{eqnarray} |
|
|
G_u & = & G_u^{adv} + G_u^{cor} + G_u^{h-diss} + G_u^{v-diss} + |
|
|
G_u^{metric} + G_u^{nh-metric} \\ |
|
|
G_v & = & G_v^{adv} + G_v^{cor} + G_v^{h-diss} + G_v^{v-diss} + |
|
|
G_v^{metric} + G_v^{nh-metric} \\ |
|
|
G_w & = & G_w^{adv} + G_w^{cor} + G_w^{h-diss} + G_w^{v-diss} + |
|
|
G_w^{metric} + G_w^{nh-metric} |
|
|
\end{eqnarray} |
|
|
In the hydrostatic limit, $G_w=0$ and $\epsilon_{nh}=0$, reducing the |
|
|
vertical momentum to hydrostatic balance. |
|
|
|
|
|
These terms are calculated in routines called from subroutine {\em |
|
|
CALC\_MOM\_RHS} a collected into the global arrays {\bf Gu}, {\bf Gv}, |
|
|
and {\bf Gw}. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R CALC\_MOM\_RHS} ({\em pkg/mom\_fluxform/calc\_mom\_rhs.F}) |
|
|
|
|
|
$G_u$: {\bf Gu} ({\em DYNVARS.h}) |
|
|
|
|
|
$G_v$: {\bf Gv} ({\em DYNVARS.h}) |
|
|
|
|
|
$G_w$: {\bf Gw} ({\em DYNVARS.h}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Advection of momentum} |
|
|
|
|
|
The advective operator is second order accurate in space: |
|
|
\begin{eqnarray} |
|
|
{\cal A}_w \Delta r_f h_w G_u^{adv} & = & |
|
|
\delta_i \overline{ U }^i \overline{ u }^i |
|
|
+ \delta_j \overline{ V }^i \overline{ u }^j |
|
|
+ \delta_k \overline{ W }^i \overline{ u }^k \\ |
|
|
{\cal A}_s \Delta r_f h_s G_v^{adv} & = & |
|
|
\delta_i \overline{ U }^j \overline{ v }^i |
|
|
+ \delta_j \overline{ V }^j \overline{ v }^j |
|
|
+ \delta_k \overline{ W }^j \overline{ v }^k \\ |
|
|
{\cal A}_c \Delta r_c G_w^{adv} & = & |
|
|
\delta_i \overline{ U }^k \overline{ w }^i |
|
|
+ \delta_j \overline{ V }^k \overline{ w }^j |
|
|
+ \delta_k \overline{ W }^k \overline{ w }^k \\ |
|
|
\end{eqnarray} |
|
|
and because of the flux form does not contribute to the global budget |
|
|
of linear momentum. The quantities $U$, $V$ and $W$ are volume fluxes |
|
|
defined: |
|
|
\marginpar{$U$: {\bf uTrans} } |
|
|
\marginpar{$V$: {\bf vTrans} } |
|
|
\marginpar{$W$: {\bf rTrans} } |
|
|
\begin{eqnarray} |
|
|
U & = & \Delta y_g \Delta r_f h_w u \\ |
|
|
V & = & \Delta x_g \Delta r_f h_s v \\ |
|
|
W & = & {\cal A}_c w |
|
|
\end{eqnarray} |
|
|
The advection of momentum takes the same form as the advection of |
|
|
tracers but by a translated advective flow. Consequently, the |
|
|
conservation of second moments, derived for tracers later, applies to |
|
|
$u^2$ and $v^2$ and $w^2$ so that advection of momentum correctly |
|
|
conserves kinetic energy. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_ADV\_UU} ({\em mom\_u\_adv\_uu.F}) |
|
|
|
|
|
{\em S/R MOM\_U\_ADV\_VU} ({\em mom\_u\_adv\_vu.F}) |
|
|
|
|
|
{\em S/R MOM\_U\_ADV\_WU} ({\em mom\_u\_adv\_wu.F}) |
|
|
|
|
|
{\em S/R MOM\_U\_ADV\_UV} ({\em mom\_u\_adv\_uv.F}) |
|
|
|
|
|
{\em S/R MOM\_U\_ADV\_VV} ({\em mom\_u\_adv\_vv.F}) |
|
|
|
|
|
{\em S/R MOM\_U\_ADV\_WV} ({\em mom\_u\_adv\_wv.F}) |
|
|
|
|
|
$uu$, $uv$, $vu$, $vv$: {\bf aF} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
|
|
|
\subsubsection{Coriolis terms} |
|
|
|
|
|
The ``pure C grid'' Coriolis terms (i.e. in absence of C-D scheme) are |
|
|
discretized: |
|
|
\begin{eqnarray} |
|
|
{\cal A}_w \Delta r_f h_w G_u^{Cor} & = & |
|
|
\overline{ f {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i |
|
|
- \epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ w }^k }^i \\ |
|
|
{\cal A}_s \Delta r_f h_s G_v^{Cor} & = & |
|
|
- \overline{ f {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\ |
|
|
{\cal A}_c \Delta r_c G_w^{Cor} & = & |
|
|
\epsilon_{nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ u }^i }^k |
|
|
\end{eqnarray} |
|
|
where the Coriolis parameters $f$ and $f'$ are defined: |
|
|
\begin{eqnarray} |
|
|
f & = & 2 \Omega \sin{\phi} \\ |
|
|
f' & = & 2 \Omega \cos{\phi} |
|
|
\end{eqnarray} |
|
|
when using spherical geometry, otherwise the $\beta$-plane definition is used: |
|
|
\begin{eqnarray} |
|
|
f & = & f_o + \beta y \\ |
|
|
f' & = & 0 |
|
|
\end{eqnarray} |
|
|
|
|
|
This discretization globally conserves kinetic energy. It should be |
|
|
noted that despite the use of this discretization in former |
|
|
publications, all calculations to date have used the following |
|
|
different discretization: |
|
|
\begin{eqnarray} |
|
|
G_u^{Cor} & = & |
|
|
f_u \overline{ v }^{ji} |
|
|
- \epsilon_{nh} f_u' \overline{ w }^{ik} \\ |
|
|
G_v^{Cor} & = & |
|
|
- f_v \overline{ u }^{ij} \\ |
|
|
G_w^{Cor} & = & |
|
|
\epsilon_{nh} f_w' \overline{ u }^{ik} |
|
|
\end{eqnarray} |
|
|
\marginpar{Need to change the default in code to match this} |
|
|
where the subscripts on $f$ and $f'$ indicate evaluation of the |
|
|
Coriolis parameters at the appropriate points in space. The above |
|
|
discretization does {\em not} conserve anything, especially energy. An |
|
|
option to recover this discretization has been retained for backward |
|
|
compatibility testing (set run-time logical {\bf |
|
|
useNonconservingCoriolis} to {\em true} which otherwise defaults to |
|
|
{\em false}). |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_CDSCHEME} ({\em mom\_cdscheme.F}) |
|
|
|
|
|
{\em S/R MOM\_U\_CORIOLIS} ({\em mom\_u\_coriolis.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_CORIOLIS} ({\em mom\_v\_coriolis.F}) |
|
|
|
|
|
$G_u^{Cor}$, $G_v^{Cor}$: {\bf cF} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Curvature metric terms} |
|
|
|
|
|
The most commonly used coordinate system on the sphere is the |
|
|
geographic system $(\lambda,\phi)$. The curvilinear nature of these |
|
|
coordinates on the sphere lead to some ``metric'' terms in the |
|
|
component momentum equations. Under the thin-atmosphere and |
|
|
hydrostatic approximations these terms are discretized: |
|
|
\begin{eqnarray} |
|
|
{\cal A}_w \Delta r_f h_w G_u^{metric} & = & |
|
|
\overline{ \frac{ \overline{u}^i }{a} \tan{\phi} {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i \\ |
|
|
{\cal A}_s \Delta r_f h_s G_v^{metric} & = & |
|
|
- \overline{ \frac{ \overline{u}^i }{a} \tan{\phi} {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\ |
|
|
G_w^{metric} & = & 0 |
|
|
\end{eqnarray} |
|
|
where $a$ is the radius of the planet (sphericity is assumed) or the |
|
|
radial distance of the particle (i.e. a function of height). It is |
|
|
easy to see that this discretization satisfies all the properties of |
|
|
the discrete Coriolis terms since the metric factor $\frac{u}{a} |
|
|
\tan{\phi}$ can be viewed as a modification of the vertical Coriolis |
|
|
parameter: $f \rightarrow f+\frac{u}{a} \tan{\phi}$. |
|
|
|
|
|
However, as for the Coriolis terms, a non-energy conserving form has |
|
|
exclusively been used to date: |
|
|
\begin{eqnarray} |
|
|
G_u^{metric} & = & \frac{u \overline{v}^{ij} }{a} \tan{\phi} \\ |
|
|
G_v^{metric} & = & \frac{ \overline{u}^{ij} \overline{u}^{ij}}{a} \tan{\phi} |
|
|
\end{eqnarray} |
|
|
where $\tan{\phi}$ is evaluated at the $u$ and $v$ points |
|
|
respectively. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_METRIC\_SPHERE} ({\em mom\_u\_metric\_sphere.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_METRIC\_SPHERE} ({\em mom\_v\_metric\_sphere.F}) |
|
|
|
|
|
$G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
|
|
|
\subsubsection{Non-hydrostatic metric terms} |
|
|
|
|
|
For the non-hydrostatic equations, dropping the thin-atmosphere |
|
|
approximation re-introduces metric terms involving $w$ and are |
|
|
required to conserve anglular momentum: |
|
|
\begin{eqnarray} |
|
|
{\cal A}_w \Delta r_f h_w G_u^{metric} & = & |
|
|
- \overline{ \frac{ \overline{u}^i \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c }^i \\ |
|
|
{\cal A}_s \Delta r_f h_s G_v^{metric} & = & |
|
|
- \overline{ \frac{ \overline{v}^j \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c}^j \\ |
|
|
{\cal A}_c \Delta r_c G_w^{metric} & = & |
|
|
\overline{ \frac{ {\overline{u}^i}^2 + {\overline{v}^j}^2}{a} {\cal A}_c \Delta r_f h_c }^k |
|
|
\end{eqnarray} |
|
|
|
|
|
Because we are always consistent, even if consistently wrong, we have, |
|
|
in the past, used a different discretization in the model which is: |
|
|
\begin{eqnarray} |
|
|
G_u^{metric} & = & |
|
|
- \frac{u}{a} \overline{w}^{ik} \\ |
|
|
G_v^{metric} & = & |
|
|
- \frac{v}{a} \overline{w}^{jk} \\ |
|
|
G_w^{metric} & = & |
|
|
\frac{1}{a} ( {\overline{u}^{ik}}^2 + {\overline{v}^{jk}}^2 ) |
|
|
\end{eqnarray} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_METRIC\_NH} ({\em mom\_u\_metric\_nh.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_METRIC\_NH} ({\em mom\_v\_metric\_nh.F}) |
|
|
|
|
|
$G_u^{metric}$, $G_v^{metric}$: {\bf mT} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Lateral dissipation} |
|
|
|
|
|
Historically, we have represented the SGS Reynolds stresses as simply |
|
|
down gradient momentum fluxes, ignoring constraints on the stress |
|
|
tensor such as symmetry. |
|
|
\begin{eqnarray} |
|
|
{\cal A}_w \Delta r_f h_w G_u^{h-diss} & = & |
|
|
\delta_i \Delta y_f \Delta r_f h_c \tau_{11} |
|
|
+ \delta_j \Delta x_v \Delta r_f h_\zeta \tau_{12} \\ |
|
|
{\cal A}_s \Delta r_f h_s G_v^{h-diss} & = & |
|
|
\delta_i \Delta y_u \Delta r_f h_\zeta \tau_{21} |
|
|
+ \delta_j \Delta x_f \Delta r_f h_c \tau_{22} |
|
|
\end{eqnarray} |
|
|
\marginpar{Check signs of stress definitions} |
|
|
|
|
|
The lateral viscous stresses are discretized: |
|
|
\begin{eqnarray} |
|
|
\tau_{11} & = & A_h c_{11\Delta}(\phi) \frac{1}{\Delta x_f} \delta_i u |
|
|
-A_4 c_{11\Delta^2}(\phi) \frac{1}{\Delta x_f} \delta_i \nabla^2 u \\ |
|
|
\tau_{12} & = & A_h c_{12\Delta}(\phi) \frac{1}{\Delta y_u} \delta_j u |
|
|
-A_4 c_{12\Delta^2}(\phi)\frac{1}{\Delta y_u} \delta_j \nabla^2 u \\ |
|
|
\tau_{21} & = & A_h c_{21\Delta}(\phi) \frac{1}{\Delta x_v} \delta_i v |
|
|
-A_4 c_{21\Delta^2}(\phi) \frac{1}{\Delta x_v} \delta_i \nabla^2 v \\ |
|
|
\tau_{22} & = & A_h c_{22\Delta}(\phi) \frac{1}{\Delta y_f} \delta_j v |
|
|
-A_4 c_{22\Delta^2}(\phi) \frac{1}{\Delta y_f} \delta_j \nabla^2 v |
|
|
\end{eqnarray} |
|
|
where the non-dimensional factors $c_{lm\Delta^n}(\phi), \{l,m,n\} \in |
|
|
\{1,2\}$ define the ``cosine'' scaling with latitude which can be |
|
|
applied in various ad-hoc ways. For instance, $c_{11\Delta} = |
|
|
c_{21\Delta} = (\cos{\phi})^{3/2}$, $c_{12\Delta}=c_{22\Delta}=0$ would |
|
|
represent the an-isotropic cosine scaling typically used on the |
|
|
``lat-lon'' grid for Laplacian viscosity. |
|
|
\marginpar{Need to tidy up method for controlling this in code} |
|
|
|
|
|
It should be noted that dispite the ad-hoc nature of the scaling, some |
|
|
scaling must be done since on a lat-lon grid the converging meridians |
|
|
make it very unlikely that a stable viscosity parameter exists across |
|
|
the entire model domain. |
|
|
|
|
|
The Laplacian viscosity coefficient, $A_h$ ({\bf viscAh}), has units |
|
|
of $m^2 s^{-1}$. The bi-harmonic viscosity coefficient, $A_4$ ({\bf |
|
|
viscA4}), has units of $m^4 s^{-1}$. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_XVISCFLUX} ({\em mom\_u\_xviscflux.F}) |
|
|
|
|
|
{\em S/R MOM\_U\_YVISCFLUX} ({\em mom\_u\_yviscflux.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_XVISCFLUX} ({\em mom\_v\_xviscflux.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_YVISCFLUX} ({\em mom\_v\_yviscflux.F}) |
|
|
|
|
|
$\tau_{11}$, $\tau_{12}$, $\tau_{22}$, $\tau_{22}$: {\bf vF}, {\bf |
|
|
v4F} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
Two types of lateral boundary condition exist for the lateral viscous |
|
|
terms, no-slip and free-slip. |
|
|
|
|
|
The free-slip condition is most convenient to code since it is |
|
|
equivalent to zero-stress on boundaries. Simple masking of the stress |
|
|
components sets them to zero. The fractional open stress is properly |
|
|
handled using the lopped cells. |
|
|
|
|
|
The no-slip condition defines the normal gradient of a tangential flow |
|
|
such that the flow is zero on the boundary. Rather than modify the |
|
|
stresses by using complicated functions of the masks and ``ghost'' |
|
|
points (see \cite{Adcroft+Marshall98}) we add the boundary stresses as |
|
|
an additional source term in cells next to solid boundaries. This has |
|
|
the advantage of being able to cope with ``thin walls'' and also makes |
|
|
the interior stress calculation (code) independent of the boundary |
|
|
conditions. The ``body'' force takes the form: |
|
|
\begin{eqnarray} |
|
|
G_u^{side-drag} & = & |
|
|
\frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta x_v}{\Delta y_u} }^j |
|
|
\left( A_h c_{12\Delta}(\phi) u - A_4 c_{12\Delta^2}(\phi) \nabla^2 u \right) |
|
|
\\ |
|
|
G_v^{side-drag} & = & |
|
|
\frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta y_u}{\Delta x_v} }^i |
|
|
\left( A_h c_{21\Delta}(\phi) v - A_4 c_{21\Delta^2}(\phi) \nabla^2 v \right) |
|
|
\end{eqnarray} |
|
|
|
|
|
In fact, the above discretization is not quite complete because it |
|
|
assumes that the bathymetry at velocity points is deeper than at |
|
|
neighbouring vorticity points, e.g. $1-h_w < 1-h_\zeta$ |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_SIDEDRAG} ({\em mom\_u\_sidedrag.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_SIDEDRAG} ({\em mom\_v\_sidedrag.F}) |
|
|
|
|
|
$G_u^{side-drag}$, $G_v^{side-drag}$: {\bf vF} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Vertical dissipation} |
|
|
|
|
|
Vertical viscosity terms are discretized with only partial adherence |
|
|
to the variable grid lengths introduced by the finite volume |
|
|
formulation. This reduces the formal accuracy of these terms to just |
|
|
first order but only next to boundaries; exactly where other terms |
|
|
appear such as linar and quadratic bottom drag. |
|
|
\begin{eqnarray} |
|
|
G_u^{v-diss} & = & |
|
|
\frac{1}{\Delta r_f h_w} \delta_k \tau_{13} \\ |
|
|
G_v^{v-diss} & = & |
|
|
\frac{1}{\Delta r_f h_s} \delta_k \tau_{23} \\ |
|
|
G_w^{v-diss} & = & \epsilon_{nh} |
|
|
\frac{1}{\Delta r_f h_d} \delta_k \tau_{33} |
|
|
\end{eqnarray} |
|
|
represents the general discrete form of the vertical dissipation terms. |
|
|
|
|
|
In the interior the vertical stresses are discretized: |
|
|
\begin{eqnarray} |
|
|
\tau_{13} & = & A_v \frac{1}{\Delta r_c} \delta_k u \\ |
|
|
\tau_{23} & = & A_v \frac{1}{\Delta r_c} \delta_k v \\ |
|
|
\tau_{33} & = & A_v \frac{1}{\Delta r_f} \delta_k w |
|
|
\end{eqnarray} |
|
|
It should be noted that in the non-hydrostatic form, the stress tensor |
|
|
is even less consistent than for the hydrostatic (see Wazjowicz |
|
|
\cite{Waojz}). It is well known how to do this properly (see Griffies |
|
|
\cite{Griffies}) and is on the list of to-do's. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_RVISCLFUX} ({\em mom\_u\_rviscflux.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_RVISCLFUX} ({\em mom\_v\_rviscflux.F}) |
|
|
|
|
|
$\tau_{13}$: {\bf urf} (local to {\em calc\_mom\_rhs.F}) |
|
|
|
|
|
$\tau_{23}$: {\bf vrf} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
As for the lateral viscous terms, the free-slip condition is |
|
|
equivalent to simply setting the stress to zero on boundaries. The |
|
|
no-slip condition is implemented as an additional term acting on top |
|
|
of the interior and free-slip stresses. Bottom drag represents |
|
|
additional friction, in addition to that imposed by the no-slip |
|
|
condition at the bottom. The drag is cast as a stress expressed as a |
|
|
linear or quadratic function of the mean flow in the layer above the |
|
|
topography: |
|
|
\begin{eqnarray} |
|
|
\tau_{13}^{bottom-drag} & = & |
|
|
\left( |
|
|
2 A_v \frac{1}{\Delta r_c} |
|
|
+ r_b |
|
|
+ C_d \sqrt{ \overline{2 KE}^i } |
|
|
\right) u \\ |
|
|
\tau_{23}^{bottom-drag} & = & |
|
|
\left( |
|
|
2 A_v \frac{1}{\Delta r_c} |
|
|
+ r_b |
|
|
+ C_d \sqrt{ \overline{2 KE}^j } |
|
|
\right) v |
|
|
\end{eqnarray} |
|
|
where these terms are only evaluated immediately above topography. |
|
|
$r_b$ ({\bf bottomDragLinear}) has units of $m s^{-1}$ and a typical value |
|
|
of the order 0.0002 $m s^{-1}$. $C_d$ ({\bf bottomDragQuadratic}) is |
|
|
dimensionless with typical values in the range 0.001--0.003. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_BOTTOMDRAG} ({\em mom\_u\_bottomdrag.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_BOTTOMDRAG} ({\em mom\_v\_bottomdrag.F}) |
|
|
|
|
|
$\tau_{13}^{bottom-drag}$, $\tau_{23}^{bottom-drag}$: {\bf vf} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
\subsection{Tracer equations} |
|
|
|
|
|
The tracer equations are discretized consistantly with the continuity |
|
|
equation to facilitate conservation properties analogous to the |
|
|
continuum: |
|
|
\begin{equation} |
|
|
{\cal A}_c \Delta r_f h_c \partial_\theta |
|
|
+ \delta_i U \overline{ \theta }^i |
|
|
+ \delta_j V \overline{ \theta }^j |
|
|
+ \delta_k W \overline{ \theta }^k |
|
|
= {\cal A}_c \Delta r_f h_c {\cal S}_\theta + \theta {\cal A}_c \delta_k (P-E)_{r=0} |
|
|
\end{equation} |
|
|
The quantities $U$, $V$ and $W$ are volume fluxes defined: |
|
|
\marginpar{$U$: {\bf uTrans} } |
|
|
\marginpar{$V$: {\bf vTrans} } |
|
|
\marginpar{$W$: {\bf rTrans} } |
|
|
\begin{eqnarray} |
|
|
U & = & \Delta y_g \Delta r_f h_w u \\ |
|
|
V & = & \Delta x_g \Delta r_f h_s v \\ |
|
|
W & = & {\cal A}_c w |
|
|
\end{eqnarray} |
|
|
${\cal S}$ represents the ``parameterized'' SGS processes and |
|
|
physics associated with the tracer. For instance, potential |
|
|
temperature equation in the ocean has is forced by surface and |
|
|
partially penetrating heat fluxes: |
|
|
\begin{equation} |
|
|
{\cal A}_c \Delta r_f h_c {\cal S}_\theta = \frac{1}{c_p \rho_o} \delta_k {\cal A}_c {\cal Q} |
|
|
\end{equation} |
|
|
while the salt equation has no real sources, ${\cal S}=0$, which |
|
|
leaves just the $P-E$ term. |
|
|
|
|
|
The continuity equation can be recovered by setting ${\cal Q}=0$ and |
|
|
$\theta=1$. The term $\theta (P-E)_{r=0}$ is required to retain local |
|
|
conservation of $\theta$. Global conservation is not possible using |
|
|
the flux-form (as here) and a linearized free-surface |
|
|
(\cite{Griffies00,Campin02}). |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
\subsection{Derivation of discrete energy conservation} |
|
|
|
|
|
These discrete equations conserve kinetic plus potential energy using the |
|
|
following definitions: |
|
|
\begin{equation} |
|
|
KE = \frac{1}{2} \left( \overline{ u^2 }^i + \overline{ v^2 }^j + |
|
|
\epsilon_{nh} \overline{ w^2 }^k \right) |
|
|
\end{equation} |
|
|
|
|
|
|
|
|
\subsection{Vector invariant momentum equations} |
|
|
|
|
|
The finite volume method lends itself to describing the continuity and |
|
|
tracer equations in curvilinear coordinate systems but the appearance |
|
|
of new metric terms in the flux-form momentum equations makes |
|
|
generalizing them far from elegant. The vector invariant form of the |
|
|
momentum equations are exactly that; invariant under coordinate |
|
|
transformations. |
|
|
|
|
|
The non-hydrostatic vector invariant equations read: |
|
|
\begin{equation} |
|
|
\partial_t \vec{v} + ( 2\vec{\Omega} + \vec{\zeta}) \wedge \vec{v} |
|
|
- b \hat{r} |
|
|
+ \vec{\nabla} B = \vec{\nabla} \cdot \vec{\bf \tau} |
|
|
\end{equation} |
|
|
which describe motions in any orthogonal curvilinear coordinate |
|
|
system. Here, $B$ is the Bernoulli function and $\vec{\zeta}=\nabla |
|
|
\wedge \vec{v}$ is the vorticity vector. We can take advantage of the |
|
|
elegance of these equations when discretizing them and use the |
|
|
discrete definitions of the grad, curl and divergence operators to |
|
|
satisfy constraints. We can also consider the analogy to forming |
|
|
derived equations, such as the vorticity equation, and examine how the |
|
|
discretization can be adjusted to give suitable vorticity advection |
|
|
among other things. |
|
|
|
|
|
The underlying algorithm is the same as for the flux form |
|
|
equations. All that has changed is the contents of the ``G's''. For |
|
|
the time-being, only the hydrostatic terms have been coded but we will |
|
|
indicate the points where non-hydrostatic contributions will enter: |
|
|
\begin{eqnarray} |
|
|
G_u & = & G_u^{fv} + G_u^{\zeta_3 v} + G_u^{\zeta_2 w} + G_u^{\partial_x B} |
|
|
+ G_u^{\partial_z \tau^x} + G_u^{h-dissip} + G_u^{v-dissip} \\ |
|
|
G_v & = & G_v^{fu} + G_v^{\zeta_3 u} + G_v^{\zeta_1 w} + G_v^{\partial_y B} |
|
|
+ G_v^{\partial_z \tau^y} + G_v^{h-dissip} + G_v^{v-dissip} \\ |
|
|
G_w & = & G_w^{fu} + G_w^{\zeta_1 v} + G_w^{\zeta_2 u} + G_w^{\partial_z B} |
|
|
+ G_w^{h-dissip} + G_w^{v-dissip} |
|
|
\end{eqnarray} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R CALC\_MOM\_RHS} ({\em pkg/mom\_vecinv/calc\_mom\_rhs.F}) |
|
|
|
|
|
$G_u$: {\bf Gu} ({\em DYNVARS.h}) |
|
|
|
|
|
$G_v$: {\bf Gv} ({\em DYNVARS.h}) |
|
|
|
|
|
$G_w$: {\bf Gw} ({\em DYNVARS.h}) |
|
|
\end{minipage} } |
|
|
|
|
|
\subsubsection{Relative vorticity} |
|
|
|
|
|
The vertical component of relative vorticity is explicitly calculated |
|
|
and use in the discretization. The particular form is crucial for |
|
|
numerical stablility; alternative definitions break the conservation |
|
|
properties of the discrete equations. |
|
|
|
|
|
Relative vorticity is defined: |
|
|
\begin{equation} |
|
|
\zeta_3 = \frac{\Gamma}{A_\zeta} |
|
|
= \frac{1}{{\cal A}_\zeta} ( \delta_i \Delta y_c v - \delta_j \Delta x_c u ) |
|
|
\end{equation} |
|
|
where ${\cal A}_\zeta$ is the area of the vorticity cell presented in |
|
|
the vertical and $\Gamma$ is the circulation about that cell. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_VI\_CALC\_RELVORT3} ({\em mom\_vi\_calc\_relvort3.F}) |
|
|
|
|
|
$\zeta_3$: {\bf vort3} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Kinetic energy} |
|
|
|
|
|
The kinetic energy, denoted $KE$, is defined: |
|
|
\begin{equation} |
|
|
KE = \frac{1}{2} ( \overline{ u^2 }^i + \overline{ v^2 }^j |
|
|
+ \epsilon_{nh} \overline{ w^2 }^k ) |
|
|
\end{equation} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_VI\_CALC\_KE} ({\em mom\_vi\_calc\_ke.F}) |
|
|
|
|
|
$KE$: {\bf KE} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Coriolis terms} |
|
|
|
|
|
The potential enstrophy conserving form of the linear Coriolis terms |
|
|
are written: |
|
|
\begin{eqnarray} |
|
|
G_u^{fv} & = & |
|
|
\frac{1}{\Delta x_c} |
|
|
\overline{ \frac{f}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\ |
|
|
G_v^{fu} & = & - |
|
|
\frac{1}{\Delta y_c} |
|
|
\overline{ \frac{f}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j |
|
|
\end{eqnarray} |
|
|
Here, the Coriolis parameter $f$ is defined at vorticity (corner) |
|
|
points. |
|
|
\marginpar{$f$: {\bf fCoriG}} |
|
|
\marginpar{$h_\zeta$: {\bf hFacZ}} |
|
|
|
|
|
The potential enstrophy conserving form of the non-linear Coriolis |
|
|
terms are written: |
|
|
\begin{eqnarray} |
|
|
G_u^{\zeta_3 v} & = & |
|
|
\frac{1}{\Delta x_c} |
|
|
\overline{ \frac{\zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\ |
|
|
G_v^{\zeta_3 u} & = & - |
|
|
\frac{1}{\Delta y_c} |
|
|
\overline{ \frac{\zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j |
|
|
\end{eqnarray} |
|
|
\marginpar{$\zeta_3$: {\bf vort3}} |
|
|
|
|
|
The Coriolis terms can also be evaluated together and expressed in |
|
|
terms of absolute vorticity $f+\zeta_3$. The potential enstrophy |
|
|
conserving form using the absolute vorticity is written: |
|
|
\begin{eqnarray} |
|
|
G_u^{fv} + G_u^{\zeta_3 v} & = & |
|
|
\frac{1}{\Delta x_c} |
|
|
\overline{ \frac{f + \zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i \\ |
|
|
G_v^{fu} + G_v^{\zeta_3 u} & = & - |
|
|
\frac{1}{\Delta y_c} |
|
|
\overline{ \frac{f + \zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j |
|
|
\end{eqnarray} |
|
|
|
|
|
\marginpar{Run-time control needs to be added for these options} The |
|
|
disctinction between using absolute vorticity or relative vorticity is |
|
|
useful when constructing higher order advection schemes; monotone |
|
|
advection of relative vorticity behaves differently to monotone |
|
|
advection of absolute vorticity. Currently the choice of |
|
|
relative/absolute vorticity, centered/upwind/high order advection is |
|
|
available only through commented subroutine calls. |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_VI\_CORIOLIS} ({\em mom\_vi\_coriolis.F}) |
|
|
|
|
|
{\em S/R MOM\_VI\_U\_CORIOLIS} ({\em mom\_vi\_u\_coriolis.F}) |
|
|
|
|
|
{\em S/R MOM\_VI\_V\_CORIOLIS} ({\em mom\_vi\_v\_coriolis.F}) |
|
|
|
|
|
$G_u^{fv}$, $G_u^{\zeta_3 v}$: {\bf uCf} (local to {\em calc\_mom\_rhs.F}) |
|
|
|
|
|
$G_v^{fu}$, $G_v^{\zeta_3 u}$: {\bf vCf} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Shear terms} |
|
|
|
|
|
The shear terms ($\zeta_2w$ and $\zeta_1w$) are are discretized to |
|
|
guarantee that no spurious generation of kinetic energy is possible; |
|
|
the horizontal gradient of Bernoulli function has to be consistent |
|
|
with the vertical advection of shear: |
|
|
\marginpar{N-H terms have not been tried!} |
|
|
\begin{eqnarray} |
|
|
G_u^{\zeta_2 w} & = & |
|
|
\frac{1}{ {\cal A}_w \Delta r_f h_w } \overline{ |
|
|
\overline{ {\cal A}_c w }^i ( \delta_k u - \epsilon_{nh} \delta_j w ) |
|
|
}^k \\ |
|
|
G_v^{\zeta_1 w} & = & |
|
|
\frac{1}{ {\cal A}_s \Delta r_f h_s } \overline{ |
|
|
\overline{ {\cal A}_c w }^i ( \delta_k u - \epsilon_{nh} \delta_j w ) |
|
|
}^k |
|
|
\end{eqnarray} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_VI\_U\_VERTSHEAR} ({\em mom\_vi\_u\_vertshear.F}) |
|
|
|
|
|
{\em S/R MOM\_VI\_V\_VERTSHEAR} ({\em mom\_vi\_v\_vertshear.F}) |
|
|
|
|
|
$G_u^{\zeta_2 w}$: {\bf uCf} (local to {\em calc\_mom\_rhs.F}) |
|
|
|
|
|
$G_v^{\zeta_1 w}$: {\bf vCf} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
|
|
|
\subsubsection{Gradient of Bernoulli function} |
|
|
|
|
|
\begin{eqnarray} |
|
|
G_u^{\partial_x B} & = & |
|
|
\frac{1}{\Delta x_c} \delta_i ( \phi' + KE ) \\ |
|
|
G_v^{\partial_y B} & = & |
|
|
\frac{1}{\Delta x_y} \delta_j ( \phi' + KE ) |
|
|
%G_w^{\partial_z B} & = & |
|
|
%\frac{1}{\Delta r_c} h_c \delta_k ( \phi' + KE ) |
|
|
\end{eqnarray} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_VI\_U\_GRAD\_KE} ({\em mom\_vi\_u\_grad\_ke.F}) |
|
|
|
|
|
{\em S/R MOM\_VI\_V\_GRAD\_KE} ({\em mom\_vi\_v\_grad\_ke.F}) |
|
|
|
|
|
$G_u^{\partial_x KE}$: {\bf uCf} (local to {\em calc\_mom\_rhs.F}) |
|
|
|
|
|
$G_v^{\partial_y KE}$: {\bf vCf} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
|
|
|
\subsubsection{Horizontal dissipation} |
|
|
|
|
|
The horizontal divergence, a complimentary quantity to relative |
|
|
vorticity, is used in parameterizing the Reynolds stresses and is |
|
|
discretized: |
|
|
\begin{equation} |
|
|
D = \frac{1}{{\cal A}_c h_c} ( |
|
|
\delta_i \Delta y_g h_w u |
|
|
+ \delta_j \Delta x_g h_s v ) |
|
|
\end{equation} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_VI\_CALC\_HDIV} ({\em mom\_vi\_calc\_hdiv.F}) |
|
|
|
|
|
$D$: {\bf hDiv} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Horizontal dissipation} |
|
|
|
|
|
The following discretization of horizontal dissipation conserves |
|
|
potential vorticity (thickness weighted relative vorticity) and |
|
|
divergence and dissipates energy, enstrophy and divergence squared: |
|
|
\begin{eqnarray} |
|
|
G_u^{h-dissip} & = & |
|
|
\frac{1}{\Delta x_c} \delta_i ( A_D D - A_{D4} D^*) |
|
|
- \frac{1}{\Delta y_u h_w} \delta_j h_\zeta ( A_\zeta \zeta - A_{\zeta4} \zeta^* ) |
|
|
\\ |
|
|
G_v^{h-dissip} & = & |
|
|
\frac{1}{\Delta x_v h_s} \delta_i h_\zeta ( A_\zeta \zeta - A_\zeta \zeta^* ) |
|
|
+ \frac{1}{\Delta y_c} \delta_j ( A_D D - A_{D4} D^* ) |
|
|
\end{eqnarray} |
|
|
where |
|
|
\begin{eqnarray} |
|
|
D^* & = & \frac{1}{{\cal A}_c h_c} ( |
|
|
\delta_i \Delta y_g h_w \nabla^2 u |
|
|
+ \delta_j \Delta x_g h_s \nabla^2 v ) \\ |
|
|
\zeta^* & = & \frac{1}{{\cal A}_\zeta} ( |
|
|
\delta_i \Delta y_c \nabla^2 v |
|
|
- \delta_j \Delta x_c \nabla^2 u ) |
|
|
\end{eqnarray} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_VI\_HDISSIP} ({\em mom\_vi\_hdissip.F}) |
|
|
|
|
|
$G_u^{h-dissip}$: {\bf uDiss} (local to {\em calc\_mom\_rhs.F}) |
|
|
|
|
|
$G_v^{h-dissip}$: {\bf vDiss} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
|
|
|
|
|
|
|
\subsubsection{Vertical dissipation} |
|
|
|
|
|
Currently, this is exactly the same code as the flux form equations. |
|
|
\begin{eqnarray} |
|
|
G_u^{v-diss} & = & |
|
|
\frac{1}{\Delta r_f h_w} \delta_k \tau_{13} \\ |
|
|
G_v^{v-diss} & = & |
|
|
\frac{1}{\Delta r_f h_s} \delta_k \tau_{23} |
|
|
\end{eqnarray} |
|
|
represents the general discrete form of the vertical dissipation terms. |
|
|
|
|
|
In the interior the vertical stresses are discretized: |
|
|
\begin{eqnarray} |
|
|
\tau_{13} & = & A_v \frac{1}{\Delta r_c} \delta_k u \\ |
|
|
\tau_{23} & = & A_v \frac{1}{\Delta r_c} \delta_k v |
|
|
\end{eqnarray} |
|
|
|
|
|
\fbox{ \begin{minipage}{4.25in} |
|
|
{\em S/R MOM\_U\_RVISCLFUX} ({\em mom\_u\_rviscflux.F}) |
|
|
|
|
|
{\em S/R MOM\_V\_RVISCLFUX} ({\em mom\_v\_rviscflux.F}) |
|
|
|
|
|
$\tau_{13}$: {\bf urf} (local to {\em calc\_mom\_rhs.F}) |
|
|
|
|
|
$\tau_{23}$: {\bf vrf} (local to {\em calc\_mom\_rhs.F}) |
|
|
\end{minipage} } |
|
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch