--- manual/s_algorithm/text/spatial-discrete.tex 2001/08/09 19:48:39 1.4 +++ manual/s_algorithm/text/spatial-discrete.tex 2001/10/24 15:21:27 1.9 @@ -1,4 +1,4 @@ -% $Header: /home/ubuntu/mnt/e9_copy/manual/s_algorithm/text/spatial-discrete.tex,v 1.4 2001/08/09 19:48:39 adcroft Exp $ +% $Header: /home/ubuntu/mnt/e9_copy/manual/s_algorithm/text/spatial-discrete.tex,v 1.9 2001/10/24 15:21:27 cnh Exp $ % $Name: $ \section{Spatial discretization of the dynamical equations} @@ -8,43 +8,22 @@ 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 -differently. +horizontal and veritical directions as seperable and differently. -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} } +\input{part2/notation} \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 -in our terminology) and the other is non-linear interpolation methods -that can deal with non-smooth solutions such as shocks (i.e. flux -limiters for advection). Both make use of the integral form of the -conservation laws to which the {\it weak solution} is a solution on -each finite volume of (sub-domain). The weak solution can be -constructed outof piece-wise constant elements or be +is the method of embedded or intersecting boundaries (shaved or lopped +cells in our terminology) and the other is non-linear interpolation +methods that can deal with non-smooth solutions such as shocks +(i.e. flux limiters for advection). Both make use of the integral form +of the conservation laws to which the {\it weak solution} is a +solution on each finite volume of (sub-domain). The weak solution can +be constructed out of piece-wise constant elements or be differentiable. The differentiable equations can not be satisfied by piece-wise constant functions. @@ -88,7 +67,9 @@ \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. } @@ -108,17 +89,44 @@ +\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 @@ -310,11 +318,13 @@ \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 @@ -379,9 +389,9 @@ \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 @@ -446,14 +456,14 @@ \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 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' \\ -\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}' \\ +\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' \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}' \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} @@ -476,14 +486,15 @@ The last equation, the discrete continuity equation, can be summed in the vertical to yeild 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(P-E)_{r=0} +{\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(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 @@ -497,12 +508,14 @@ \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