/[MITgcm]/manual/s_algorithm/text/spatial-discrete.tex
ViewVC logotype

Diff of /manual/s_algorithm/text/spatial-discrete.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by adcroft, Wed Aug 8 22:19:02 2001 UTC revision 1.23 by jmc, Fri Aug 27 13:08:18 2010 UTC
# Line 2  Line 2 
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  Spatial discretization is carried out using the finite volume
10  method. This amounts to a grid-point method (namely second-order  method. This amounts to a grid-point method (namely second-order
11  centered finite difference) in the fluid interior but allows  centered finite difference) in the fluid interior but allows
12  boundaries to intersect a regular grid allowing a more accurate  boundaries to intersect a regular grid allowing a more accurate
13  representation of the position of the boundary. We treat the  representation of the position of the boundary. We treat the
14  horizontal and veritical directions as seperable and thus slightly  horizontal and vertical directions as separable and differently.
 differently.  
15    
 Initialization of grid data is controlled by subroutine {\em  
 INI\_GRID} which in calls {\em INI\_VERTICAL\_GRID} to initialize the  
 vertical grid, and then either of {\em INI\_CARTESIAN\_GRID}, {\em  
 INI\_SPHERICAL\_POLAR\_GRID} or {\em INI\_CURV\-ILINEAR\_GRID} to  
 initialize the horizontal grid for cartesian, spherical-polar or  
 curvilinear coordinates respectively.  
   
 The reciprocals of all grid quantities are pre-calculated and this is  
 done in subroutine {\em INI\_MASKS\_ETC} which is called later by  
 subroutine {\em INITIALIZE\_FIXED}.  
   
 All grid descriptors are global arrays and stored in common blocks in  
 {\em GRID.h} and a generally declared as {\em \_RS}.  
   
 \fbox{ \begin{minipage}{4.75in}  
 {\em S/R INI\_GRID} ({\em model/src/ini\_grid.F})  
16    
17  {\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F})  \subsection{The finite volume method: finite volumes versus finite difference}
18    \begin{rawhtml}
19  grid data: ({\em model/inc/GRID.h})  <!-- CMIREDIR:finite_volume: -->
20  \end{minipage} }  \end{rawhtml}
21    
22    
 \subsection{The finite volume method: finite volumes versus finite difference}  
23    
24  The finite volume method is used to discretize the equations in  The finite volume method is used to discretize the equations in
25  space. The expression ``finite volume'' actually has two meanings; one  space. The expression ``finite volume'' actually has two meanings; one
26  involves invocation of the weak formulation (e.g. integral  is the method of embedded or intersecting boundaries (shaved or lopped
27  formulation) and the other involves non-linear expressions for  cells in our terminology) and the other is non-linear interpolation
28  interpolation of data (e.g. flux limiters). We use both but they can  methods that can deal with non-smooth solutions such as shocks
29  and will be ddiscussed seperately. The finite volume method discretizes by invoking the weak formulation of the equations or integral form. For example, the 1-D advection-diffusion equation:  (i.e. flux limiters for advection). Both make use of the integral form
30    of the conservation laws to which the {\it weak solution} is a
31    solution on each finite volume of (sub-domain). The weak solution can
32    be constructed out of piece-wise constant elements or be
33    differentiable. The differentiable equations can not be satisfied by
34    piece-wise constant functions.
35    
36    As an example, the 1-D constant coefficient advection-diffusion
37    equation:
38  \begin{displaymath}  \begin{displaymath}
39  \partial_t \theta + \partial_x ( u \theta - \kappa \partial_x \theta ) = 0  \partial_t \theta + \partial_x ( u \theta - \kappa \partial_x \theta ) = 0
40  \end{displaymath}  \end{displaymath}
41  can be discretized by integrating of finite lengths $\Delta x$:  can be discretized by integrating over finite sub-domains, i.e.
42    the lengths $\Delta x_i$:
43  \begin{displaymath}  \begin{displaymath}
44  \Delta x \partial_t \theta + \delta_i ( F ) = 0  \Delta x \partial_t \theta + \delta_i ( F ) = 0
45  \end{displaymath}  \end{displaymath}
46  is exact but where the flux  is exact if $\theta(x)$ is piece-wise constant over the interval
47    $\Delta x_i$ or more generally if $\theta_i$ is defined as the average
48    over the interval $\Delta x_i$.
49    
50    The flux, $F_{i-1/2}$, must be approximated:
51  \begin{displaymath}  \begin{displaymath}
52  F = u \overline{\theta} - \frac{\kappa}{\Delta x_c} \partial_i \theta  F = u \overline{\theta} - \frac{\kappa}{\Delta x_c} \partial_i \theta
53  \end{displaymath}  \end{displaymath}
54  is approximate. The method for obtained $\overline{\theta}$ is  and this is where truncation errors can enter the solution. The
55  unspecified and non-lienar finite volume methods can be invoked.  method for obtaining $\overline{\theta}$ is unspecified and a wide
56    range of possibilities exist including centered and upwind
57  ....  INCOMPLETE  interpolation, polynomial fits based on the the volume average
58    definitions of quantities and non-linear interpolation such as
59    flux-limiters.
60    
61    Choosing simple centered second-order interpolation and differencing
62    recovers the same ODE's resulting from finite differencing for the
63    interior of a fluid. Differences arise at boundaries where a boundary
64    is not positioned on a regular or smoothly varying grid. This method
65    is used to represent the topography using lopped cell, see
66    \cite{adcroft:97}. Subtle difference also appear in more than one
67    dimension away from boundaries. This happens because the each
68    direction is discretized independently in the finite difference method
69    while the integrating over finite volume implicitly treats all
70    directions simultaneously. Illustration of this is given in
71    \cite{ac:02}.
72    
73  \subsection{C grid staggering of variables}  \subsection{C grid staggering of variables}
74    
75  \begin{figure}  \begin{figure}
76  \centerline{ \resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} }  \begin{center}
77    \resizebox{!}{2in}{ \includegraphics{s_algorithm/figs/cgrid3d.eps}}
78    \end{center}
79  \caption{Three dimensional staggering of velocity components. This  \caption{Three dimensional staggering of velocity components. This
80  facilitates the natural discretization of the continuity and tracer  facilitates the natural discretization of the continuity and tracer
81  equations. }  equations. }
# Line 71  equations. } Line 85  equations. }
85  The basic algorithm employed for stepping forward the momentum  The basic algorithm employed for stepping forward the momentum
86  equations is based on retaining non-divergence of the flow at all  equations is based on retaining non-divergence of the flow at all
87  times. This is most naturally done if the components of flow are  times. This is most naturally done if the components of flow are
88  staggered in space in the form of an Arakawa C grid \cite{Arakawa70}.  staggered in space in the form of an Arakawa C grid \cite{arakawa:77}.
89    
90  Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$)  Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$)
91  staggered in space such that the zonal component falls on the  staggered in space such that the zonal component falls on the
92  interface between continiuty cells in the zonal direction. Similarly  interface between continuity cells in the zonal direction. Similarly
93  for the meridional and vertical directions.  The continiuty cell is  for the meridional and vertical directions.  The continuity cell is
94  synonymous with tracer cells (they are one and the same).  synonymous with tracer cells (they are one and the same).
95    
96    
97    
98    \subsection{Grid initialization and data}
99    
100    Initialization of grid data is controlled by subroutine {\em
101    INI\_GRID} which in calls {\em INI\_VERTICAL\_GRID} to initialize the
102    vertical grid, and then either of {\em INI\_CARTESIAN\_GRID}, {\em
103    INI\_SPHERICAL\_POLAR\_GRID} or {\em INI\_CURV\-ILINEAR\_GRID} to
104    initialize the horizontal grid for cartesian, spherical-polar or
105    curvilinear coordinates respectively.
106    
107    The reciprocals of all grid quantities are pre-calculated and this is
108    done in subroutine {\em INI\_MASKS\_ETC} which is called later by
109    subroutine {\em INITIALIZE\_FIXED}.
110    
111    All grid descriptors are global arrays and stored in common blocks in
112    {\em GRID.h} and a generally declared as {\em \_RS}.
113    
114    \fbox{ \begin{minipage}{4.75in}
115    {\em S/R INI\_GRID} ({\em model/src/ini\_grid.F})
116    
117    {\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F})
118    
119    grid data: ({\em model/inc/GRID.h})
120    \end{minipage} }
121    
122    
123  \subsection{Horizontal grid}  \subsection{Horizontal grid}
124    \label{sec:spatial_discrete_horizontal_grid}
125    
126  \begin{figure}  \begin{figure}
127  \centerline{ \begin{tabular}{cc}  \begin{center}
128    \raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}}  \begin{tabular}{cc}
129  & \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}}    \raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{s_algorithm/figs/hgrid-Ac.eps}}
130    & \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{s_algorithm/figs/hgrid-Az.eps}}
131  \\  \\
132    \raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}}    \raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{s_algorithm/figs/hgrid-Au.eps}}
133  & \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}}  & \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{s_algorithm/figs/hgrid-Av.eps}}
134  \end{tabular} }  \end{tabular}
135    \end{center}
136  \caption{  \caption{
137  Staggering of horizontal grid descriptors (lengths and areas). The  Staggering of horizontal grid descriptors (lengths and areas). The
138  grid lines indicate the tracer cell boundaries and are the reference  grid lines indicate the tracer cell boundaries and are the reference
139  grid for all panels. a) The area of a tracer cell, $A_c$, is bordered  grid for all panels. a) The area of a tracer cell, $A_c$, is bordered
140  by the lengths $\Delta x_g$ and $\Delta y_g$. b) The area of a  by the lengths $\Delta x_g$ and $\Delta y_g$. b) The area of a
141  vorticity cell, $A_\zeta$, is bordered by the lengths $\Delta x_c$ and  vorticity cell, $A_\zeta$, is bordered by the lengths $\Delta x_c$ and
142  $\Delta y_c$. c) The area of a u cell, $A_c$, is bordered by the  $\Delta y_c$. c) The area of a u cell, $A_w$, is bordered by the
143  lengths $\Delta x_v$ and $\Delta y_f$. d) The area of a v cell, $A_c$,  lengths $\Delta x_v$ and $\Delta y_f$. d) The area of a v cell, $A_s$,
144  is bordered by the lengths $\Delta x_f$ and $\Delta y_u$.}  is bordered by the lengths $\Delta x_f$ and $\Delta y_u$.}
145  \label{fig:hgrid}  \label{fig:hgrid}
146  \end{figure}  \end{figure}
147    
148  The model domain is decomposed into tiles and within each tile a  The model domain is decomposed into tiles and within each tile a
149  quasi-regular grid is used. A tile is the basic unit of domain  quasi-regular grid is used. A tile is the basic unit of domain
150  decomposition for parallelization but may be used whether parallized  decomposition for parallelization but may be used whether parallelized
151  or not; see section \ref{sect:tiles} for more details. Although the  or not; see section \ref{sect:domain_decomposition} for more details.
152  tiles may be patched together in an unstructured manner  Although the tiles may be patched together in an unstructured manner
153  (i.e. irregular or non-tessilating pattern), the interior of tiles is  (i.e. irregular or non-tessilating pattern), the interior of tiles is
154  a structered grid of quadrilateral cells. The horizontal coordinate  a structured grid of quadrilateral cells. The horizontal coordinate
155  system is orthogonal curvilinear meaning we can not necessarily treat  system is orthogonal curvilinear meaning we can not necessarily treat
156  the two horizontal directions as seperable. Instead, each cell in the  the two horizontal directions as separable. Instead, each cell in the
157  horizontal grid is described by the length of it's sides and it's  horizontal grid is described by the length of it's sides and it's
158  area.  area.
159    
160  The grid information is quite general and describes any of the  The grid information is quite general and describes any of the
161  available coordinates systems, cartesian, spherical-polar or  available coordinates systems, cartesian, spherical-polar or
162  curvilinear. All that is necessary to distinguish between the  curvilinear. All that is necessary to distinguish between the
163  coordinate systems is to initialize the grid data (discriptors)  coordinate systems is to initialize the grid data (descriptors)
164  appropriately.  appropriately.
165    
 \marginpar{Caution!}  
166  In the following, we refer to the orientation of quantities on the  In the following, we refer to the orientation of quantities on the
167  computational grid using geographic terminology such as points of the  computational grid using geographic terminology such as points of the
168  compass. This is purely for convenience but should note be confused  compass.
169    \marginpar{Caution!}
170    This is purely for convenience but should note be confused
171  with the actual geographic orientation of model quantities.  with the actual geographic orientation of model quantities.
172    
 \marginpar{$A_c$: {\bf rAc}}  
 \marginpar{$\Delta x_g$: {\bf DXg}}  
 \marginpar{$\Delta y_g$: {\bf DYg}}  
173  Fig.~\ref{fig:hgrid}a shows the tracer cell (synonymous with the  Fig.~\ref{fig:hgrid}a shows the tracer cell (synonymous with the
174  continuity cell). The length of the southern edge, $\Delta x_g$,  continuity cell). The length of the southern edge, $\Delta x_g$,
175  western edge, $\Delta y_g$ and surface area, $A_c$, presented in the  western edge, $\Delta y_g$ and surface area, $A_c$, presented in the
176  vertical are stored in arrays {\bf DXg}, {\bf DYg} and {\bf rAc}. The  vertical are stored in arrays {\bf DXg}, {\bf DYg} and {\bf rAc}.
177  ``g'' suffix indicates that the lengths are along the defining grid  \marginpar{$A_c$: {\bf rAc}}
178  boundaries. The ``c'' suffix associates the quantity with the cell  \marginpar{$\Delta x_g$: {\bf DXg}}
179  centers. The quantities are staggered in space and the indexing is  \marginpar{$\Delta y_g$: {\bf DYg}}
180  such that {\bf DXg(i,j)} is positioned to the south of {\bf rAc(i,j)}  The ``g'' suffix indicates that the lengths are along the defining
181  and {\bf DYg(i,j)} positioned to the west.  grid boundaries. The ``c'' suffix associates the quantity with the
182    cell centers. The quantities are staggered in space and the indexing
183    is such that {\bf DXg(i,j)} is positioned to the south of {\bf
184    rAc(i,j)} and {\bf DYg(i,j)} positioned to the west.
185    
 \marginpar{$A_\zeta$: {\bf rAz}}  
 \marginpar{$\Delta x_c$: {\bf DXc}}  
 \marginpar{$\Delta y_c$: {\bf DYc}}  
186  Fig.~\ref{fig:hgrid}b shows the vorticity cell. The length of the  Fig.~\ref{fig:hgrid}b shows the vorticity cell. The length of the
187  southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface  southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface
188  area, $A_\zeta$, presented in the vertical are stored in arrays {\bf  area, $A_\zeta$, presented in the vertical are stored in arrays {\bf
189  DXg}, {\bf DYg} and {\bf rAz}.  The ``z'' suffix indicates that the  DXc}, {\bf DYc} and {\bf rAz}.
190  lengths are measured between the cell centers and the ``$\zeta$'' suffix  \marginpar{$A_\zeta$: {\bf rAz}}
191  associates points with the vorticity points. The quantities are  \marginpar{$\Delta x_c$: {\bf DXc}}
192  staggered in space and the indexing is such that {\bf DXc(i,j)} is  \marginpar{$\Delta y_c$: {\bf DYc}}
193  positioned to the north of {\bf rAc(i,j)} and {\bf DYc(i,j)} positioned  The ``z'' suffix indicates that the lengths are measured between the
194  to the east.  cell centers and the ``$\zeta$'' suffix associates points with the
195    vorticity points. The quantities are staggered in space and the
196    indexing is such that {\bf DXc(i,j)} is positioned to the north of
197    {\bf rAc(i,j)} and {\bf DYc(i,j)} positioned to the east.
198    
 \marginpar{$A_w$: {\bf rAw}}  
 \marginpar{$\Delta x_v$: {\bf DXv}}  
 \marginpar{$\Delta y_f$: {\bf DYf}}  
199  Fig.~\ref{fig:hgrid}c shows the ``u'' or western (w) cell. The length of  Fig.~\ref{fig:hgrid}c shows the ``u'' or western (w) cell. The length of
200  the southern edge, $\Delta x_v$, eastern edge, $\Delta y_f$ and  the southern edge, $\Delta x_v$, eastern edge, $\Delta y_f$ and
201  surface area, $A_w$, presented in the vertical are stored in arrays  surface area, $A_w$, presented in the vertical are stored in arrays
202  {\bf DXv}, {\bf DYf} and {\bf rAw}.  The ``v'' suffix indicates that  {\bf DXv}, {\bf DYf} and {\bf rAw}.
203  the length is measured between the v-points, the ``f'' suffix  \marginpar{$A_w$: {\bf rAw}}
204  indicates that the length is measured between the (tracer) cell faces  \marginpar{$\Delta x_v$: {\bf DXv}}
205  and the ``w'' suffix associates points with the u-points (w stands for  \marginpar{$\Delta y_f$: {\bf DYf}}
206  west). The quantities are staggered in space and the indexing is such  The ``v'' suffix indicates that the length is measured between the
207  that {\bf DXv(i,j)} is positioned to the south of {\bf rAw(i,j)} and  v-points, the ``f'' suffix indicates that the length is measured
208  {\bf DYf(i,j)} positioned to the east.  between the (tracer) cell faces and the ``w'' suffix associates points
209    with the u-points (w stands for west). The quantities are staggered in
210    space and the indexing is such that {\bf DXv(i,j)} is positioned to
211    the south of {\bf rAw(i,j)} and {\bf DYf(i,j)} positioned to the east.
212    
 \marginpar{$A_s$: {\bf rAs}}  
 \marginpar{$\Delta x_f$: {\bf DXf}}  
 \marginpar{$\Delta y_u$: {\bf DYu}}  
213  Fig.~\ref{fig:hgrid}d shows the ``v'' or southern (s) cell. The length of  Fig.~\ref{fig:hgrid}d shows the ``v'' or southern (s) cell. The length of
214  the northern edge, $\Delta x_f$, western edge, $\Delta y_u$ and  the northern edge, $\Delta x_f$, western edge, $\Delta y_u$ and
215  surface area, $A_s$, presented in the vertical are stored in arrays  surface area, $A_s$, presented in the vertical are stored in arrays
216  {\bf DXf}, {\bf DYu} and {\bf rAs}.  The ``u'' suffix indicates that  {\bf DXf}, {\bf DYu} and {\bf rAs}.
217  the length is measured between the u-points, the ``f'' suffix  \marginpar{$A_s$: {\bf rAs}}
218  indicates that the length is measured between the (tracer) cell faces  \marginpar{$\Delta x_f$: {\bf DXf}}
219  and the ``s'' suffix associates points with the v-points (s stands for  \marginpar{$\Delta y_u$: {\bf DYu}}
220  south). The quantities are staggered in space and the indexing is such  The ``u'' suffix indicates that the length is measured between the
221  that {\bf DXf(i,j)} is positioned to the north of {\bf rAs(i,j)} and  u-points, the ``f'' suffix indicates that the length is measured
222  {\bf DYu(i,j)} positioned to the west.  between the (tracer) cell faces and the ``s'' suffix associates points
223    with the v-points (s stands for south). The quantities are staggered
224    in space and the indexing is such that {\bf DXf(i,j)} is positioned to
225    the north of {\bf rAs(i,j)} and {\bf DYu(i,j)} positioned to the west.
226    
227  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
228  {\em S/R INI\_CARTESIAN\_GRID} ({\em  {\em S/R INI\_CARTESIAN\_GRID} ({\em
# Line 257  using\-Cartes\-ianGrid} in namelist {\em Line 299  using\-Cartes\-ianGrid} in namelist {\em
299  spacing can be set to uniform via scalars {\bf dXspacing} and {\bf  spacing can be set to uniform via scalars {\bf dXspacing} and {\bf
300  dYspacing} in namelist {\em PARM04} or to variable resolution by the  dYspacing} in namelist {\em PARM04} or to variable resolution by the
301  vectors {\bf DELX} and {\bf DELY}. Units are normally  vectors {\bf DELX} and {\bf DELY}. Units are normally
302  meters. Non-dimensional coordinates can be used by interpretting the  meters. Non-dimensional coordinates can be used by interpreting the
303  gravitational constant as the Rayleigh number.  gravitational constant as the Rayleigh number.
304    
305  \subsubsection{Spherical-polar coordinates}  \subsubsection{Spherical-polar coordinates}
# Line 282  other grids, the horizontal grid descrip Line 324  other grids, the horizontal grid descrip
324  \subsection{Vertical grid}  \subsection{Vertical grid}
325    
326  \begin{figure}  \begin{figure}
327  \centerline{ \begin{tabular}{cc}  \begin{center}
328      \begin{tabular}{cc}
329    \raisebox{4in}{a)} \resizebox{!}{4in}{    \raisebox{4in}{a)} \resizebox{!}{4in}{
330    \includegraphics{part2/vgrid-cellcentered.eps}} & \raisebox{4in}{b)}    \includegraphics{s_algorithm/figs/vgrid-cellcentered.eps}} & \raisebox{4in}{b)}
331    \resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}}    \resizebox{!}{4in}{ \includegraphics{s_algorithm/figs/vgrid-accurate.eps}}
332  \end{tabular} }  \end{tabular}
333    \end{center}
334  \caption{Two versions of the vertical grid. a) The cell centered  \caption{Two versions of the vertical grid. a) The cell centered
335  approach where the interface depths are specified and the tracer  approach where the interface depths are specified and the tracer
336  points centered in between the interfaces. b) The interface centered  points centered in between the interfaces. b) The interface centered
# Line 309  The vertical grid is calculated in subro Line 353  The vertical grid is calculated in subro
353  INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in  INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in
354  namelist {\em PARM04}. The units of ``r'' are either meters or Pascals  namelist {\em PARM04}. The units of ``r'' are either meters or Pascals
355  depending on the isomorphism being used which in turn is dependent  depending on the isomorphism being used which in turn is dependent
356  only on the choise of equation of state.  only on the choice of equation of state.
357    
358  There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which  There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which
359  dictate whether z- or  dictate whether z- or
# Line 323  vertical grid descriptors are stored in Line 367  vertical grid descriptors are stored in
367    
368  The above grid (Fig.~\ref{fig:vgrid}a) is known as the cell centered  The above grid (Fig.~\ref{fig:vgrid}a) is known as the cell centered
369  approach because the tracer points are at cell centers; the cell  approach because the tracer points are at cell centers; the cell
370  centers are mid-way between the cell interfaces. An alternative, the  centers are mid-way between the cell interfaces.
371  vertex or interface centered approach, is shown in  This discretization is selected when the thickness of the
372    levels are provided ({\bf delR}, parameter file {\em data},
373    namelist {\em PARM04})
374    An alternative, the vertex or interface centered approach, is shown in
375  Fig.~\ref{fig:vgrid}b. Here, the interior interfaces are positioned  Fig.~\ref{fig:vgrid}b. Here, the interior interfaces are positioned
376  mid-way between the tracer nodes (no longer cell centers). This  mid-way between the tracer nodes (no longer cell centers). This
377  approach is formally more accurate for evaluation of hydrostatic  approach is formally more accurate for evaluation of hydrostatic
378  pressure and vertical advection but historically the cell centered  pressure and vertical advection but historically the cell centered
379  approach has been used. An alternative form of subroutine {\em  approach has been used. An alternative form of subroutine {\em
380  INI\_VERTICAL\_GRID} is used to select the interface centered approach  INI\_VERTICAL\_GRID} is used to select the interface centered approach
381  but no run time option is currently available.  This form requires to specify $Nr+1$ vertical distances {\bf delRc}
382    (parameter file {\em data}, namelist {\em PARM04}, e.g.
383    {\em verification/ideal\_2D\_oce/input/data})
384    corresponding to surface to center, $Nr-1$ center to center, and center to
385    bottom distances.
386    %but no run time option is currently available.
387    
388  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
389  {\em S/R INI\_VERTICAL\_GRID} ({\em  {\em S/R INI\_VERTICAL\_GRID} ({\em
# Line 348  $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e Line 400  $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e
400  \end{minipage} }  \end{minipage} }
401    
402    
403  \subsection{Continuity and horizontal pressure gradient terms}  \subsection{Topography: partially filled cells}
404    \begin{rawhtml}
405    <!-- CMIREDIR:topo_partial_cells: -->
406    \end{rawhtml}
407    
408    \begin{figure}
409    \begin{center}
410    \resizebox{4.5in}{!}{\includegraphics{s_algorithm/figs/vgrid-xz.eps}}
411    \end{center}
412    \caption{
413    A schematic of the x-r plane showing the location of the
414    non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a
415    tracer cell is given by $h_c(i,j,k) \Delta r_f(k)$ and the physical
416    thickness of the open side is given by $h_w(i,j,k) \Delta r_f(k)$.}
417    \label{fig:hfacs}
418    \end{figure}
419    
420    \cite{adcroft:97} presented two alternatives to the step-wise finite
421    difference representation of topography. The method is known to the
422    engineering community as {\em intersecting boundary method}. It
423    involves allowing the boundary to intersect a grid of cells thereby
424    modifying the shape of those cells intersected. We suggested allowing
425    the topography to take on a piece-wise linear representation (shaved
426    cells) or a simpler piecewise constant representation (partial step).
427    Both show dramatic improvements in solution compared to the
428    traditional full step representation, the piece-wise linear being the
429    best. However, the storage requirements are excessive so the simpler
430    piece-wise constant or partial-step method is all that is currently
431    supported.
432    
433    Fig.~\ref{fig:hfacs} shows a schematic of the x-r plane indicating how
434    the thickness of a level is determined at tracer and u points.
435    \marginpar{$h_c$: {\bf hFacC}}
436    \marginpar{$h_w$: {\bf hFacW}}
437    \marginpar{$h_s$: {\bf hFacS}}
438    The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta
439    r_f(k)$ and the physical thickness of the open side is given by
440    $h_w(i,j,k) \Delta r_f(k)$. Three 3-D descriptors $h_c$, $h_w$ and
441    $h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and
442    {\bf hFacS} respectively. These are calculated in subroutine {\em
443    INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf
444    RECIP\_hFacW} and {\bf RECIP\_hFacS}.
445    
446    The non-dimensional fractions (or h-facs as we call them) are
447    calculated from the model depth array and then processed to avoid tiny
448    volumes. The rule is that if a fraction is less than {\bf hFacMin}
449    then it is rounded to the nearer of $0$ or {\bf hFacMin} or if the
450    physical thickness is less than {\bf hFacMinDr} then it is similarly
451    rounded. The larger of the two methods is used when there is a
452    conflict. By setting {\bf hFacMinDr} equal to or larger than the
453    thinnest nominal layers, $\min{(\Delta z_f)}$, but setting {\bf
454    hFacMin} to some small fraction then the model will only lop thick
455    layers but retain stability based on the thinnest unlopped thickness;
456    $\min{(\Delta z_f,\mbox{\bf hFacMinDr})}$.
457    
458    \fbox{ \begin{minipage}{4.75in}
459    {\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F})
460    
461    $h_c$: {\bf hFacC} ({\em GRID.h})
462    
463    $h_w$: {\bf hFacW} ({\em GRID.h})
464    
465    $h_s$: {\bf hFacS} ({\em GRID.h})
466    
467    $h_c^{-1}$: {\bf RECIP\_hFacC} ({\em GRID.h})
468    
469    $h_w^{-1}$: {\bf RECIP\_hFacW} ({\em GRID.h})
470    
471    $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRID.h})
472    
473    \end{minipage} }
474    
475    
476    \section{Continuity and horizontal pressure gradient terms}
477    \begin{rawhtml}
478    <!-- CMIREDIR:continuity_and_horizontal_pressure: -->
479    \end{rawhtml}
480    
481    
482  The core algorithm is based on the ``C grid'' discretization of the  The core algorithm is based on the ``C grid'' discretization of the
483  continuity equation which can be summarized as:  continuity equation which can be summarized as:
484  \begin{eqnarray}  \begin{eqnarray}
485  \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} \\
486  \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} \\
487  \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} \\
488  \delta_i \Delta y_g \Delta r_f h_w u +  \delta_i \Delta y_g \Delta r_f h_w u +
489  \delta_j \Delta x_g \Delta r_f h_s v +  \delta_j \Delta x_g \Delta r_f h_s v +
490  \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}
# Line 363  continuity equation which can be summari Line 492  continuity equation which can be summari
492  \end{eqnarray}  \end{eqnarray}
493  where the continuity equation has been most naturally discretized by  where the continuity equation has been most naturally discretized by
494  staggering the three components of velocity as shown in  staggering the three components of velocity as shown in
495  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$
496  are the lengths between tracer points (cell centers). The grid lengths  are the lengths between tracer points (cell centers). The grid lengths
497  $\Delta x_g$, $\Delta y_g$ are the grid lengths between cell  $\Delta x_g$, $\Delta y_g$ are the grid lengths between cell
498  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
# Line 376  A}_c$.  The factors $h_w$ and $h_s$ are Line 505  A}_c$.  The factors $h_w$ and $h_s$ are
505  \marginpar{$h_s$: {\bf hFacS}}  \marginpar{$h_s$: {\bf hFacS}}
506    
507  The last equation, the discrete continuity equation, can be summed in  The last equation, the discrete continuity equation, can be summed in
508  the vertical to yeild the free-surface equation:  the vertical to yield the free-surface equation:
509  \begin{equation}  \begin{equation}
510  {\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
511  {\cal A}_c(P-E)_{r=0}  u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal
512    A}_c(P-E)_{r=0} \label{eq:discrete-freesurface}
513  \end{equation}  \end{equation}
514  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
515  addition of volume due to excess precipitation and run-off over  addition of volume due to excess precipitation and run-off over
516  evaporation and only enters the top-level of the {\em ocean} model.  evaporation and only enters the top-level of the {\em ocean} model.
517    
518  \subsection{Hydrostatic balance}  \section{Hydrostatic balance}
519    \begin{rawhtml}
520    <!-- CMIREDIR:hydrostatic_balance: -->
521    \end{rawhtml}
522    
523  The vertical momentum equation has the hydrostatic or  The vertical momentum equation has the hydrostatic or
524  quasi-hydrostatic balance on the right hand side. This discretization  quasi-hydrostatic balance on the right hand side. This discretization
# Line 394  derived from the buoyancy equation exact Line 527  derived from the buoyancy equation exact
527  from the pressure gradient terms when forming the kinetic energy  from the pressure gradient terms when forming the kinetic energy
528  equation.  equation.
529    
530  In the ocean, using z-ccordinates, the hydrostatic balance terms are  In the ocean, using z-coordinates, the hydrostatic balance terms are
531  discretized:  discretized:
532  \begin{equation}  \begin{equation}
533  \epsilon_{nh} \partial_t w  \epsilon_{nh} \partial_t w
534  + 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
535    \label{eq:discrete_hydro_ocean}
536  \end{equation}  \end{equation}
537    
538  In the atmosphere, using p-coordinates, hydrostatic balance is  In the atmosphere, using p-coordinates, hydrostatic balance is
539  discretized:  discretized:
540  \begin{equation}  \begin{equation}
541  \overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0  \overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0
542    \label{eq:discrete_hydro_atmos}
543  \end{equation}  \end{equation}
544  where $\Delta \Pi$ is the difference in Exner function between the  where $\Delta \Pi$ is the difference in Exner function between the
545  pressure points. The non-hydrostatic equations are not available in  pressure points. The non-hydrostatic equations are not available in
# Line 412  the atmosphere. Line 547  the atmosphere.
547    
548  The difference in approach between ocean and atmosphere occurs because  The difference in approach between ocean and atmosphere occurs because
549  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
550  energy conversion term $\alpha \omega$. The form of these consversion  energy conversion term $\alpha \omega$. The form of these conversion
551  terms is discussed at length in \cite{Adcroft01}.  terms is discussed at length in \cite{adcroft:02}.
552    
553  Because of the different representation of hydrostatic balance between  Because of the different representation of hydrostatic balance between
554  ocean and atmosphere there is no elegant way to represent both systems  ocean and atmosphere there is no elegant way to represent both systems
# Line 428  CALC\_PHI\_HYD}. Inside this routine, on Line 563  CALC\_PHI\_HYD}. Inside this routine, on
563  atmospheric/oceanic form is selected based on the string variable {\bf  atmospheric/oceanic form is selected based on the string variable {\bf
564  buoyancyRelation}.  buoyancyRelation}.
565    
 \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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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.75in}  
 {\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} }  
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22