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

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

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

-revision 1.3 by adcroft,
Thu Aug  9 16:26:43 2001 UTC
+revision 1.14 by afe,
Tue Mar 23 16:47:04 2004 UTC
 Line 8 
 method. This amounts to a grid-point met
  centered finite difference) in the fluid interior but allows
  boundaries to intersect a regular grid allowing a more accurate
  representation of the position of the boundary. We treat the
- horizontal and veritical directions as seperable and thus slightly
+ horizontal and vertical directions as separable and differently.
- differently.
- Initialization of grid data is controlled by subroutine {\em
+ \input{part2/notation}
- 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
+ \subsection{The finite volume method: finite volumes versus finite difference}
- {\em GRID.h} and a generally declared as {\em \_RS}.
+ \begin{rawhtml}
+ <!-- CMIREDIR:finite_volume: -->
- \fbox{ \begin{minipage}{4.75in}
+ \end{rawhtml}
- {\em S/R INI\_GRID} ({\em model/src/ini\_grid.F})
- {\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F})
- grid data: ({\em model/inc/GRID.h})
- \end{minipage} }
- \subsection{The finite volume method: finite volumes versus finite difference}
  The finite volume method is used to discretize the equations in
  space. The expression ``finite volume'' actually has two meanings; one
- is the method of cut or instecting boundaries (shaved or lopped cells
+ is the method of embedded or intersecting boundaries (shaved or lopped
- in our terminology) and the other is non-linear interpolation methods
+ cells in our terminology) and the other is non-linear interpolation
- that can deal with non-smooth solutions such as shocks (i.e. flux
+ methods that can deal with non-smooth solutions such as shocks
- limiters for advection). Both make use of the integral form of the
+ (i.e. flux limiters for advection). Both make use of the integral form
- conservation laws to which the {\it weak solution} is a solution on
+ of the conservation laws to which the {\it weak solution} is a
- each finite volume of (sub-domain). The weak solution can be
+ solution on each finite volume of (sub-domain). The weak solution can
- constructed outof piece-wise constant elements or be
+ be constructed out of piece-wise constant elements or be
  differentiable. The differentiable equations can not be satisfied by
  piece-wise constant functions.
-Line 58 
 the lengths $\Delta x_i$:
+Line 42 
 the lengths $\Delta x_i$:
  \begin{displaymath}
  \Delta x \partial_t \theta + \delta_i ( F ) = 0
  \end{displaymath}
- is exact if $\theta(x)$ is peice-wise constant over the interval
+ is exact if $\theta(x)$ is piece-wise constant over the interval
  $\Delta x_i$ or more generally if $\theta_i$ is defined as the average
  over the interval $\Delta x_i$.
-Line 78 
 recovers the same ODE's resulting from f
+Line 62 
 recovers the same ODE's resulting from f
  interior of a fluid. Differences arise at boundaries where a boundary
  is not positioned on a regular or smoothly varying grid. This method
  is used to represent the topography using lopped cell, see
- \cite{Adcroft98}. Subtle difference also appear in more than one
+ \cite{adcroft:97}. Subtle difference also appear in more than one
  dimension away from boundaries. This happens because the each
- direction is discretized independantly in the finite difference method
+ direction is discretized independently in the finite difference method
  while the integrating over finite volume implicitly treats all
  directions simultaneously. Illustration of this is given in
- \cite{Adcroft02}.
+ \cite{ac:02}.
  \subsection{C grid staggering of variables}
  \begin{figure}
- \centerline{ \resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} }
+ \begin{center}
+ \resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}}
+ \end{center}
  \caption{Three dimensional staggering of velocity components. This
  facilitates the natural discretization of the continuity and tracer
  equations. }
-Line 98 
 equations. }
+Line 84 
 equations. }
  The basic algorithm employed for stepping forward the momentum
  equations is based on retaining non-divergence of the flow at all
  times. This is most naturally done if the components of flow are
- 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}.
  Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$)
  staggered in space such that the zonal component falls on the
- interface between continiuty cells in the zonal direction. Similarly
+ interface between continuity cells in the zonal direction. Similarly
- for the meridional and vertical directions.  The continiuty cell is
+ for the meridional and vertical directions.  The continuity cell is
  synonymous with tracer cells (they are one and the same).
+ \subsection{Grid initialization and data}
+ 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})
+ {\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F})
+ grid data: ({\em model/inc/GRID.h})
+ \end{minipage} }
  \subsection{Horizontal grid}
+ \label{sec:spatial_discrete_horizontal_grid}
  \begin{figure}
- \centerline{ \begin{tabular}{cc}
+ \begin{center}
+ \begin{tabular}{cc}
    \raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}}
  & \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}}
  \\
    \raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}}
  & \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}}
- \end{tabular} }
+ \end{tabular}
+ \end{center}
  \caption{
  Staggering of horizontal grid descriptors (lengths and areas). The
  grid lines indicate the tracer cell boundaries and are the reference
-Line 133 
 is bordered by the lengths $\Delta x_f$
+Line 146 
 is bordered by the lengths $\Delta x_f$
  The model domain is decomposed into tiles and within each tile a
  quasi-regular grid is used. A tile is the basic unit of domain
- decomposition for parallelization but may be used whether parallized
+ decomposition for parallelization but may be used whether parallelized
  or not; see section \ref{sect:tiles} for more details. Although the
  tiles may be patched together in an unstructured manner
  (i.e. irregular or non-tessilating pattern), the interior of tiles is
- a structered grid of quadrilateral cells. The horizontal coordinate
+ a structured grid of quadrilateral cells. The horizontal coordinate
  system is orthogonal curvilinear meaning we can not necessarily treat
- the two horizontal directions as seperable. Instead, each cell in the
+ the two horizontal directions as separable. Instead, each cell in the
  horizontal grid is described by the length of it's sides and it's
  area.
  The grid information is quite general and describes any of the
  available coordinates systems, cartesian, spherical-polar or
  curvilinear. All that is necessary to distinguish between the
- coordinate systems is to initialize the grid data (discriptors)
+ coordinate systems is to initialize the grid data (descriptors)
  appropriately.
  In the following, we refer to the orientation of quantities on the
-Line 285 
 using\-Cartes\-ianGrid} in namelist {\em
+Line 298 
 using\-Cartes\-ianGrid} in namelist {\em
  spacing can be set to uniform via scalars {\bf dXspacing} and {\bf
  dYspacing} in namelist {\em PARM04} or to variable resolution by the
  vectors {\bf DELX} and {\bf DELY}. Units are normally
- meters. Non-dimensional coordinates can be used by interpretting the
+ meters. Non-dimensional coordinates can be used by interpreting the
  gravitational constant as the Rayleigh number.
  \subsubsection{Spherical-polar coordinates}
-Line 310 
 other grids, the horizontal grid descrip
+Line 323 
 other grids, the horizontal grid descrip
  \subsection{Vertical grid}
  \begin{figure}
- \centerline{ \begin{tabular}{cc}
+ \begin{center}
+   \begin{tabular}{cc}
    \raisebox{4in}{a)} \resizebox{!}{4in}{
    \includegraphics{part2/vgrid-cellcentered.eps}} & \raisebox{4in}{b)}
    \resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}}
- \end{tabular} }
+ \end{tabular}
+ \end{center}
  \caption{Two versions of the vertical grid. a) The cell centered
  approach where the interface depths are specified and the tracer
  points centered in between the interfaces. b) The interface centered
-Line 337 
 The vertical grid is calculated in subro
+Line 352 
 The vertical grid is calculated in subro
  INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in
  namelist {\em PARM04}. The units of ``r'' are either meters or Pascals
  depending on the isomorphism being used which in turn is dependent
- only on the choise of equation of state.
+ only on the choice of equation of state.
  There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which
  dictate whether z- or
-Line 379 
 $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e
+Line 394 
 $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e
  \subsection{Topography: partially filled cells}
  \begin{figure}
- \centerline{
+ \begin{center}
  \resizebox{4.5in}{!}{\includegraphics{part2/vgrid-xz.eps}}
- }
+ \end{center}
  \caption{
  A schematic of the x-r plane showing the location of the
  non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a
-Line 390 
 thickness of the open side is given by $
+Line 405 
 thickness of the open side is given by $
  \label{fig:hfacs}
  \end{figure}
- \cite{Adcroft97} presented two alternatives to the step-wise finite
+ \cite{adcroft:97} presented two alternatives to the step-wise finite
  difference representation of topography. The method is known to the
  engineering community as {\em intersecting boundary method}. It
  involves allowing the boundary to intersect a grid of cells thereby
  modifying the shape of those cells intersected. We suggested allowing
- the topgoraphy to take on a peice-wise linear representation (shaved
+ the topography to take on a piece-wise linear representation (shaved
- cells) or a simpler piecewise constant representaion (partial step).
+ cells) or a simpler piecewise constant representation (partial step).
  Both show dramatic improvements in solution compared to the
  traditional full step representation, the piece-wise linear being the
  best. However, the storage requirements are excessive so the simpler
-Line 410 
 the thickness of a level is determined a
+Line 425 
 the thickness of a level is determined a
  \marginpar{$h_s$: {\bf hFacS}}
  The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta
  r_f(k)$ and the physical thickness of the open side is given by
- $h_w(i,j,k) \Delta r_f(k)$. Three 3-D discriptors $h_c$, $h_w$ and
+ $h_w(i,j,k) \Delta r_f(k)$. Three 3-D descriptors $h_c$, $h_w$ and
  $h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and
  {\bf hFacS} respectively. These are calculated in subroutine {\em
  INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf
-Line 446 
 $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRI
+Line 461 
 $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRI
  \end{minipage} }
- \subsection{Continuity and horizontal pressure gradient terms}
+ \section{Continuity and horizontal pressure gradient terms}
  The core algorithm is based on the ``C grid'' discretization of the
  continuity equation which can be summarized as:
  \begin{eqnarray}
- \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} \\
- \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} \\
- \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} \\
  \delta_i \Delta y_g \Delta r_f h_w u +
  \delta_j \Delta x_g \Delta r_f h_s v +
  \delta_k {\cal A}_c w & = & {\cal A}_c \delta_k (P-E)_{r=0}
-Line 461 
 continuity equation which can be summari
+Line 476 
 continuity equation which can be summari
  \end{eqnarray}
  where the continuity equation has been most naturally discretized by
  staggering the three components of velocity as shown in
- 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$
  are the lengths between tracer points (cell centers). The grid lengths
  $\Delta x_g$, $\Delta y_g$ are the grid lengths between cell
  corners. $\Delta r_f$ and $\Delta r_c$ are the distance (in units of
-Line 474 
 A}_c$.  The factors $h_w$ and $h_s$ are
+Line 489 
 A}_c$.  The factors $h_w$ and $h_s$ are
  \marginpar{$h_s$: {\bf hFacS}}
  The last equation, the discrete continuity equation, can be summed in
- the vertical to yeild the free-surface equation:
+ the vertical to yield the free-surface equation:
  \begin{equation}
- {\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
- {\cal A}_c(P-E)_{r=0}
+ u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal
+ A}_c(P-E)_{r=0} \label{eq:discrete-freesurface}
  \end{equation}
  The source term $P-E$ on the rhs of continuity accounts for the local
  addition of volume due to excess precipitation and run-off over
  evaporation and only enters the top-level of the {\em ocean} model.
- \subsection{Hydrostatic balance}
+ \section{Hydrostatic balance}
  The vertical momentum equation has the hydrostatic or
  quasi-hydrostatic balance on the right hand side. This discretization
-Line 492 
 derived from the buoyancy equation exact
+Line 508 
 derived from the buoyancy equation exact
  from the pressure gradient terms when forming the kinetic energy
  equation.
- In the ocean, using z-ccordinates, the hydrostatic balance terms are
+ In the ocean, using z-coordinates, the hydrostatic balance terms are
  discretized:
  \begin{equation}
  \epsilon_{nh} \partial_t w
  + g \overline{\rho'}^k + \frac{1}{\Delta z} \delta_k \Phi_h' = \ldots
+ \label{eq:discrete_hydro_ocean}
  \end{equation}
  In the atmosphere, using p-coordinates, hydrostatic balance is
  discretized:
  \begin{equation}
  \overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0
+ \label{eq:discrete_hydro_atmos}
  \end{equation}
  where $\Delta \Pi$ is the difference in Exner function between the
  pressure points. The non-hydrostatic equations are not available in
-Line 510 
 the atmosphere.
+Line 528 
 the atmosphere.
  The difference in approach between ocean and atmosphere occurs because
  of the direct use of the ideal gas equation in forming the potential
- energy conversion term $\alpha \omega$. The form of these consversion
+ energy conversion term $\alpha \omega$. The form of these conversion
- terms is discussed at length in \cite{Adcroft01}.
+ terms is discussed at length in \cite{adcroft:02}.
  Because of the different representation of hydrostatic balance between
  ocean and atmosphere there is no elegant way to represent both systems
-Line 526 
 CALC\_PHI\_HYD}. Inside this routine, on
+Line 544 
 CALC\_PHI\_HYD}. Inside this routine, on
  atmospheric/oceanic form is selected based on the string variable {\bf
  buoyancyRelation}.
- \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(
-A_v \frac{1}{\Delta r_c}
- + r_b
- + C_d \sqrt{ \overline{2 KE}^i }
- \right) u \\
- \tau_{23}^{bottom-drag} & = &
- \left(
-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.3
 


changed lines


 
Added in v.1.14
 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.14
-Removed from v.1.3
+Added in v.1.14

	ViewVC Help
Powered by ViewVC 1.1.22