<p><b>qchen3@fsu.edu</b> 2012-11-27 10:28:04 -0700 (Tue, 27 Nov 2012)</p><p>Some minor cleaning up and rearrangements on the GM design document.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/documents/ocean/current_design_doc/gm/gm.tex
===================================================================
--- trunk/documents/ocean/current_design_doc/gm/gm.tex 2012-11-27 17:07:19 UTC (rev 2322)
+++ trunk/documents/ocean/current_design_doc/gm/gm.tex 2012-11-27 17:28:04 UTC (rev 2323)
@@ -156,19 +156,11 @@
&{\bf U}_k = u_k + u^\ast_k,\label{eq:7}\\
& W_k = w_k + w^\ast_k,\label{eq:8},
\end{align}
-and the $3\times 3$ tensor $\Kb$ and the two dimensional slope $\Sb$
-are given by
-\begin{equation}
- \label{eq:9}
- \Kb = \left(
- \begin{aligned}
- &\mathbb{I} & \Sb\\
- &\Sb^T & \Sb\cdot\Sb
- \end{aligned}\right), \quad \Sb = -\dfrac{</font>
<font color="gray">abla_z \rho}{\rho_z},
-\end{equation}
-where $\mathbb{I}$ is the $2\times 2$ identity matrix.
+and $\Kb$ is a $3\times 3$ tensor and will be given later in Section
+\ref{sec:redi-isopycn-diff}.
-\section{The stream-function formulation of the GM parametrizations}
+\section{The stream-function formulation of the GM
+ parametrizations}\label{sec:stre-funct-form}
Let
\begin{equation}
\label{eq:10}
@@ -190,19 +182,27 @@
\label{eq:11}
\bs{\gamma}_\textrm{GM} = -\kappa\Sb,
\end{equation}
-where $\Sb$ is a 2D vector representing the slope of the isopycnal surfaces. Then
+where $\Sb$ is a 2D vector representing the slope of the isopycnal
+surfaces, and it is given by
+\begin{equation}
+ \label{eq:9}
+ \Sb = -\dfrac{</font>
<font color="blue">abla_z \rho}{\rho_z}.
+\end{equation}
+Then
\begin{align}
&\ub^\ast_\textrm{GM} = -\p_z(\kappa\Sb), \label{eq:12}\\
&w^\ast_\textrm{GM} = </font>
<font color="red">abla_z\cdot(\kappa \Sb).
\end{align}
-We note that $w^\ast$ can be determined via (\ref{eq:14}) or using (\ref{eq:13}) along with the condition that
+We note that $w^\ast$ can be determined via (\ref{eq:14}) or using
+(\ref{eq:13}) along with the incompressibility constraint
\begin{equation}
\label{eq:2}
- </font>
<font color="blue">abla \cdot \vb \equiv 0.
+ </font>
<font color="gray">abla \cdot \vb = 0.
\end{equation}
-\section{Tapering by solving a boundary value problem}
+\section{Tapering by solving a boundary value
+ problem}\label{sec:taper-solv-bound}
Following \cite{Ferrari10}, one solves a boundary value
problem for the horizontal 2D vector
$\bs{\gamma}$ in \eqref{eq:13} and \eqref{eq:14},
@@ -220,10 +220,6 @@
and bottom thin boundary layers where $\bs{\gamma}$ goes to zero
while $\bs{\gamma}_\textrm{GM}$ does not.
-\section{Normal components}
-Date last modified: 2012/9/18 \\
-Contributors: Qingshan Chen, Todd Ringler \\
-
MPAS-O utilizes a C-grid staggering technique. Hence it is necessary
to recast the full continuous GM formulation, given in the previous
section, in terms of the normal components of the variables. In what
@@ -250,115 +246,7 @@
Here $(\p \rho/\p\mathrm{n})|_z$ refers to the normal derivative of
the density field on constant height surfaces.
-%-----------------------------------------------------------------------
-
-\chapter{Design and Implementation}
-
-
-\section{Placement of the discrete variables}
-\begin{figure}[h]
- \centering
- \scalebox{0.8}{\includegraphics{variable-placement.pdf}}
- \caption{The placement of the discrete variables}
- \label{fig:1}
-\end{figure}
-We assume that a general vertical coordinate $r$ is in use. In the 3D
-MPAS-O model, the scalar variables $h,\, T,\, S,\, \rho$ etc.~are
-defined at the cell centers and in the middle of the layers. The
-horizontal normal velocity component $u$ is defined at the cell edges
-and in the middle of the layers. The vertical transport velocity $w$
-is defined at the cell centers and at the layer interfaces. The normal
-derivative of the density field is then most conveniently calculated
-on the cell edges and in the middle of the layers. The vertical
-derivative $\rho_z$ is most conveniently calculated at the cell
-centers and at the layer interfaces. The calculations of the normal
-component of the slope $\bs{s}$ of the constant density surfaces or
-the calculations of the normal component of the vector $\bs{\gamma}$
-involve both $\rho_z$ and $\p\rho/\p\mathrm{n}$, which suggests that
-$s$ and/or $\gamma$ should be specified on the cell edges and at the
-layer interfaces.
-
-\section{Formulae}
-Date last modified: 2012/09/21 \\
-Contributors: Qingshan Chen and Todd Ringler\\
-
-We label the layers from the top down by $k=1,2,3,\cdots$, and the
-layer interface between layers $k$ and $k+1$ by $k+1/2$. The height
-of the midpoint of layer $k$ is denoted as $z_k$, and the height of
-the layer interface $k+1/2$ is denoted as $z_{k+1/2}$. A first order
-approximation to $\rho_z$ at layer interface $k+1/2$ is given by
-\begin{equation}
- \left(\rho_z\right)_{k+\frac{1}{2}} =
- \dfrac{\rho_k - \rho_{k+1}}{z_k - z_{k+1}}.\label{eq:21}
-\end{equation}
-The calculation of $\p\rho/\p\mathrm{n}$ on constant height surfaces is
-a bit more complicated, because the general vertical coordinate usually
-does not align with the traditional $z$-level coordinate. We begin by
-deriving a formula for $</font>
<font color="black">abla_z\rho$ on the general coordinate. Projecting $</font>
<font color="black">abla_z\rho$ in the direction of $</font>
<font color="red">b$ will give $\p\rho/\p\mathrm{n}$ on
-constant height surfaces.
-
-We call the general vertical coordinate $r$. Thus it can be written that
-\begin{align}
-& r = r(x,y,z),\label{eq:22}\\
-& z = z(x,y,r).\label{eq:23}
-\end{align}
-The time variable $t$ has been dropped for now, since it is not relevant
-in the derivation. Applying the horizontal gradient operator $</font>
<font color="red">abla_z$
-to $\rho(x,y,r(x,y,z))$, and using the chain rule, we obtain
-\begin{equation}
-</font>
<font color="black">abla_z \rho = </font>
<font color="black">abla_r\rho + \dfrac{\p\rho}{\p r}</font>
<font color="red">abla_z r.
-\label{eq:24}
-\end{equation}
-Applying $</font>
<font color="red">abla_z$ to \eqref{eq:23}, we obtain
-\begin{equation*}
-0 = </font>
<font color="black">abla_r z + \dfrac{\p z}{\p r}</font>
<font color="red">abla_z r.
-\end{equation*}
-Hence
-\begin{equation}
-</font>
<font color="black">abla_z r = \frac{-</font>
<font color="red">abla_r z}{\frac{\p z}{\p r}}.\label{eq:25}
-\end{equation}
-Applying $\p/\p z$ to \eqref{eq:23}, we have
-\begin{equation}
-1 = \dfrac{\p z}{\p r}\dfrac{\p r}{\p z}.\label{eq:26}
-\end{equation}
-Combining \eqref{eq:24}--\eqref{eq:26} yields
-\begin{equation}
-</font>
<font color="black">abla_z \rho = </font>
<font color="red">abla_r \rho - \dfrac{\p p}{\p r}
-\dfrac{\p r}{\p z}</font>
<font color="red">abla_r z.\label{eq:27}
-\end{equation}
-We notice that
-\begin{equation}
-\dfrac{\p \rho}{\p z} = \dfrac{\p \rho}{\p r}
-\dfrac{\p r}{\p z}\label{eq:28}
-\end{equation}
-by taking the $z$-derivative of $\rho(x,y,r(x,y,z))$. Hence
-\eqref{eq:27} can be written as
-\begin{equation}
-</font>
<font color="black">abla_z \rho = </font>
<font color="black">abla_r \rho - \dfrac{\p\rho}{\p z}</font>
<font color="red">abla_r z.
-\label{eq:29}
-\end{equation}
-Taking dot-product of \eqref{eq:29} with the outer normal unit vector
-$</font>
<font color="gray">b$ yields
-\begin{equation}
-\dfrac{\p\rho}{\p \mathrm{n}}\rvert_z =
-\dfrac{\p\rho}{\p\mathrm{n}}\rvert_r - \dfrac{\p\rho}{\p z}
-\dfrac{\p z}{\p\mathrm{n}}\rvert_r.
-\label{eq:30}
-\end{equation}
-The discretization of the RHS of \eqref{eq:30} are as follow,
-\begin{align*}
-\dfrac{\p\rho}{\p n}\rvert_r &=
-\dfrac{\rho(\textrm{cell2}) - \rho(\textrm{cell1})}{\textrm{dcEdge}},\\
-\dfrac{\p z}{\p n}\rvert_r &=
-\dfrac{z(\textrm{cell2}) - z(\textrm{cell1})}{\textrm{dcEdge}}.
-\end{align*}
-$\p\rho/\p z$ has been calculated before at the cell centers and on
-the layer interfaces. The data first needs to be interpolated to the cell
-edges on the layer interfaces using area weighted averaging, and then
-be interpolated to the cell edges in the middle of the layer
-interfaces.
-
-\section{Redi isopycnal diffusion}
+\section{Redi isopycnal diffusion}\label{sec:redi-isopycn-diff}
Date last modified: 2012/11/19 \\
Contributors: Qingshan Chen, Todd Ringler, and Peter Gent\\
@@ -621,19 +509,20 @@
K^r = \cos^2\beta\left(
\begin{matrix}
1 + k_y^2 & { } & -k_xk_y & { } & { }\\
- { } & { } & { } & { } & \cos\gamma{\Large\mathbf{S}}\\
+ { } & { } & { } & { } & \cos\gamma{\tilde{\mathbf{S}}}\\
-k_xk_y & { } & 1 + k_x^2 & { } & { }\\
{ } & { } & { } & { } & { } \\
- { } & \cos\gamma\mathbf{S}^\textrm{T} & { } & { } &
+ { } & \cos\gamma\tilde{\mathbf{S}}^\textrm{T} & { } & { } &
\cos^2\gamma K_{33}
\end{matrix}\right),
\end{equation}
with
\begin{align*}
-& \mathbf{S} = \left(k_x-(1+k_y^2)l_x + k_xk_yl_y,\, k_y-(1+k_x^2)l_y
- + k_xk_yl_x\right)^\textrm{T},\\
-&K_{33} = l_x^2(1+k_y^2) + l_y^2(1+k_x^2) - 2(k_xl_x + k_yl_y) + k_x^2
-+ k_y^2 - 2k_xk_yl_xl_y.
+& \tilde{\Sb} = \left(k_x- l_x + k_y(k_xl_y - k_yl_x),\, k_y-l_y
+ + k_x(k_yl_x - k_xl_y)\right)^\textrm{T},\\
+&K_{33} = (k_x - l_x)^2 + (k_y - l_y)^2 + (k_yl_x - k_xl_y)^2.
+% l_x^2(1+k_y^2) + l_y^2(1+k_x^2) - 2(k_xl_x + k_yl_y) + k_x^2
+% + k_y^2 - 2k_xk_yl_xl_y.
\end{align*}
</font>
<font color="gray">oindent{\it Sanity checking}\\
@@ -690,26 +579,24 @@
\label{eq:46}
(l_x, l_y) = </font>
<font color="red">abla_r z
\end{equation}
-Using \eqref{eq:45} and \eqref{eq:46}, the vector $\mathbf{S}$ can
-rewritten and simplified as
+Using \eqref{eq:45} and \eqref{eq:46}, the vector $\tilde{\Sb}$ can
+written as
\begin{equation}
\label{eq:47}
- \mathbf{S} = -\dfrac{</font>
<font color="red">abla_r\rho}{\rho_z} +
- \left(-k_y^2l_x + k_xk_yl_y, -k_x^2l_y +
- k_xk_yl_x\right)^\textrm{T}.
+ \tilde{\Sb} = -\dfrac{</font>
<font color="red">abla_r\rho}{\rho_z} +
+ \left( k_y(k_xl_y - k_yl_x),\, k_x(k_yl_x - k_xl_y)\right)^\textrm{T}.
\end{equation}
-The term $K_{33}$ can be rewritten as simplified as
-\begin{align*}
- K_{33} &= (k_x-l_x)^2 + (k_y - l_y)^2 + l_x^2k_y^2 + l_y^2k_x^2 -
- 2k_xk_yl_xl_y\\
- &= \dfrac{|</font>
<font color="red">abla_r\rho|^2}{\rho_z^2}+ l_x^2k_y^2 + l_y^2k_x^2 -
- 2k_xk_yl_xl_y.
-\end{align*}
-Hence,
+The term $K_{33}$ can be written as
+% \begin{align*}
+% K_{33} &= (k_x-l_x)^2 + (k_y - l_y)^2 + l_x^2k_y^2 + l_y^2k_x^2 -
+% 2k_xk_yl_xl_y\\
+% &= \dfrac{|</font>
<font color="red">abla_r\rho|^2}{\rho_z^2}+ l_x^2k_y^2 + l_y^2k_x^2 -
+% 2k_xk_yl_xl_y.
+% \end{align*}
+% Hence,
\begin{equation}
\label{eq:49}
- K_{33} = \dfrac{|</font>
<font color="red">abla_r\rho|^2}{\rho_z^2}+ l_x^2k_y^2 +
- l_y^2k_x^2 - 2k_xk_yl_xl_y
+ K_{33} = \dfrac{|</font>
<font color="black">abla_r\rho|^2}{\rho_z^2}+ (k_yl_x - k_xl_y)^2.
\end{equation}
</font>
<font color="gray">oindent{\it Small angle approximation}\\
@@ -724,32 +611,33 @@
\begin{equation}
\label{eq:53}
K^r = \left(\begin{matrix}
- I_2 & { } & \cos\gamma\mathbf{S}\\
+ I_2 & { } & \cos\gamma\tilde{\Sb}\\
\\
- \cos\gamma\mathbf{S} & { } & \cos^2\gamma K_{33}\end{matrix}\right),
+ \cos\gamma\tilde{\mathbf{S}} & { } & \cos^2\gamma K_{33}\end{matrix}\right),
\end{equation}
with
\begin{align}
- &\mathbf{S} =
+ &\tilde{\mathbf{S}} =
-\dfrac{</font>
<font color="black">abla_r\rho}{\rho_z}, \label{eq:51}\\
&K_{33} = \dfrac{|</font>
<font color="red">abla_r\rho|^2}{\rho_z^2}\label{eq:52}
\end{align}
It can be easily checked that the diffusion tensor $K^r$ with
-$\mathbf{S}$ and $K_{33}$ taking the simplified forms \eqref{eq:51}
+$\tilde{\mathbf{S}}$ and $K_{33}$ taking the simplified forms \eqref{eq:51}
and \eqref{eq:52} has no effect on the density field, that is,
\begin{displaymath}
K^r</font>
<font color="blue">abla^r\rho = 0.
\end{displaymath}
+
</font>
<font color="black">oindent{\it Implementation}\\
With $K^r$ taking the form of \eqref{eq:53}, the Redi
diffusion can be rearranged as
\begin{multline}
\label{eq:54}
</font>
<font color="black">abla^r\cdot(K^r</font>
<font color="black">abla^r\varphi) = </font>
<font color="black">abla_r\cdot(</font>
<font color="red">abla_r\varphi)
- + </font>
<font color="blue">abla_r\cdot\left(\cos\gamma\mathbf{S}\dfrac{\p\varphi}{\p
+ + </font>
<font color="red">abla_r\cdot\left(\cos\gamma\tilde{\mathbf{S}}\dfrac{\p\varphi}{\p
r}\right) +\\ \dfrac{\p}{\p
- r}\left(\cos\gamma\mathbf{S}\cdot</font>
<font color="blue">abla_r\varphi\right) +
+ r}\left(\cos\gamma\tilde{\mathbf{S}}\cdot</font>
<font color="gray">abla_r\varphi\right) +
\dfrac{\p}{\p r}\left(\cos^2\gamma K_{33}\dfrac{\p}{\p
r}\varphi\right),
\end{multline}
@@ -767,9 +655,9 @@
\begin{multline}
\label{eq:55}
</font>
<font color="black">abla^r\cdot(K^r</font>
<font color="black">abla^r\varphi) = </font>
<font color="black">abla_r\cdot(</font>
<font color="red">abla_r\varphi)
- + </font>
<font color="blue">abla_r\cdot\left(\mathbf{S}\dfrac{\p\varphi}{\p
+ + </font>
<font color="red">abla_r\cdot\left(\tilde{\mathbf{S}}\dfrac{\p\varphi}{\p
z}\right) + \dfrac{\p}{\p
- z}\left(\mathbf{S}\cdot</font>
<font color="blue">abla_r\varphi\right) +
+ z}\left(\tilde{\mathbf{S}}\cdot</font>
<font color="gray">abla_r\varphi\right) +
\dfrac{\p}{\p z}\left( K_{33}\dfrac{\p}{\p
z}\varphi\right).
\end{multline}
@@ -780,17 +668,133 @@
$\mathbf{S}\cdot</font>
<font color="red">abla_r\varphi$ at the cell centers, which can be
approximated by
\begin{align*}
- \left[\mathbf{S}\cdot</font>
<font color="red">abla_r\varphi\right]_i &= 2\sum_{e\in
- EC(i)}\mathbf{S}\cdot\mathbf{n}_e \dfrac{\p\varphi}{\p n_e}\\
+ \left[\tilde{\mathbf{S}}\cdot</font>
<font color="blue">abla_r\varphi\right]_i &= 2\sum_{e\in
+ EC(i)}\tilde{\mathbf{S}}\cdot\mathbf{n}_e {</font>
<font color="red">abla_r\varphi}\cdot\mathbf{n}_e\\
&= -2\sum_{e\in
- EC(i)}\dfrac{1}{[\rho_z]_e}\dfrac{\p\rho}{\p n_e} \dfrac{\p\varphi}{\p n_e},
+ EC(i)}\dfrac{1}{[\rho_z]_e}\left.\dfrac{\p\rho}{\p
+ n_e}\right\vert_r \left.\dfrac{\p\varphi}{\p n_e}\right|_r,
\end{align*}
-with $\p/\p n_e$ denoting the derivative in the normal direction on
+with $(\p/\p n_e)|_r$ denoting the derivative in the normal direction
+along constant $r$ surfaces on
edge $e$.
+% \section{Normal components}\label{sec:normal-components}
+% Date last modified: 2012/9/18 \\
+% Contributors: Qingshan Chen, Todd Ringler \\
+
+
%-----------------------------------------------------------------------
-\chapter{Testing}
+\chapter{Design and Implementation}\label{cha:design-impl}
+
+
+\section{Placement of the discrete variables}\label{sec:plac-discr-vari}
+\begin{figure}[h]
+ \centering
+ \scalebox{0.8}{\includegraphics{variable-placement.pdf}}
+ \caption{The placement of the discrete variables}
+ \label{fig:1}
+\end{figure}
+We assume that a general vertical coordinate $r$ is in use. In the 3D
+MPAS-O model, the scalar variables $h,\, T,\, S,\, \rho$ etc.~are
+defined at the cell centers and in the middle of the layers. The
+horizontal normal velocity component $u$ is defined at the cell edges
+and in the middle of the layers. The vertical transport velocity $w$
+is defined at the cell centers and at the layer interfaces. The normal
+derivative of the density field is then most conveniently calculated
+on the cell edges and in the middle of the layers. The vertical
+derivative $\rho_z$ is most conveniently calculated at the cell
+centers and at the layer interfaces. The calculations of the normal
+component of the slope $\bs{s}$ of the constant density surfaces or
+the calculations of the normal component of the vector $\bs{\gamma}$
+involve both $\rho_z$ and $\p\rho/\p\mathrm{n}$, which suggests that
+$s$ and/or $\gamma$ should be specified on the cell edges and at the
+layer interfaces.
+
+\section{Formulae}\label{sec:formulae}
+Date last modified: 2012/09/21 \\
+Contributors: Qingshan Chen and Todd Ringler\\
+
+We label the layers from the top down by $k=1,2,3,\cdots$, and the
+layer interface between layers $k$ and $k+1$ by $k+1/2$. The height
+of the midpoint of layer $k$ is denoted as $z_k$, and the height of
+the layer interface $k+1/2$ is denoted as $z_{k+1/2}$. A first order
+approximation to $\rho_z$ at layer interface $k+1/2$ is given by
+\begin{equation}
+ \left(\rho_z\right)_{k+\frac{1}{2}} =
+ \dfrac{\rho_k - \rho_{k+1}}{z_k - z_{k+1}}.\label{eq:21}
+\end{equation}
+The calculation of $\p\rho/\p\mathrm{n}$ on constant height surfaces is
+a bit more complicated, because the general vertical coordinate usually
+does not align with the traditional $z$-level coordinate. We begin by
+deriving a formula for $</font>
<font color="black">abla_z\rho$ on the general coordinate. Projecting $</font>
<font color="black">abla_z\rho$ in the direction of $</font>
<font color="blue">b$ will give $\p\rho/\p\mathrm{n}$ on
+constant height surfaces.
+
+We call the general vertical coordinate $r$. Thus it can be written that
+\begin{align}
+& r = r(x,y,z),\label{eq:22}\\
+& z = z(x,y,r).\label{eq:23}
+\end{align}
+The time variable $t$ has been dropped for now, since it is not relevant
+in the derivation. Applying the horizontal gradient operator $</font>
<font color="blue">abla_z$
+to $\rho(x,y,r(x,y,z))$, and using the chain rule, we obtain
+\begin{equation}
+</font>
<font color="black">abla_z \rho = </font>
<font color="black">abla_r\rho + \dfrac{\p\rho}{\p r}</font>
<font color="blue">abla_z r.
+\label{eq:24}
+\end{equation}
+Applying $</font>
<font color="blue">abla_z$ to \eqref{eq:23}, we obtain
+\begin{equation*}
+0 = </font>
<font color="black">abla_r z + \dfrac{\p z}{\p r}</font>
<font color="blue">abla_z r.
+\end{equation*}
+Hence
+\begin{equation}
+</font>
<font color="black">abla_z r = \frac{-</font>
<font color="blue">abla_r z}{\frac{\p z}{\p r}}.\label{eq:25}
+\end{equation}
+Applying $\p/\p z$ to \eqref{eq:23}, we have
+\begin{equation}
+1 = \dfrac{\p z}{\p r}\dfrac{\p r}{\p z}.\label{eq:26}
+\end{equation}
+Combining \eqref{eq:24}--\eqref{eq:26} yields
+\begin{equation}
+</font>
<font color="black">abla_z \rho = </font>
<font color="blue">abla_r \rho - \dfrac{\p p}{\p r}
+\dfrac{\p r}{\p z}</font>
<font color="blue">abla_r z.\label{eq:27}
+\end{equation}
+We notice that
+\begin{equation}
+\dfrac{\p \rho}{\p z} = \dfrac{\p \rho}{\p r}
+\dfrac{\p r}{\p z}\label{eq:28}
+\end{equation}
+by taking the $z$-derivative of $\rho(x,y,r(x,y,z))$. Hence
+\eqref{eq:27} can be written as
+\begin{equation}
+</font>
<font color="black">abla_z \rho = </font>
<font color="black">abla_r \rho - \dfrac{\p\rho}{\p z}</font>
<font color="blue">abla_r z.
+\label{eq:29}
+\end{equation}
+Taking dot-product of \eqref{eq:29} with the outer normal unit vector
+$</font>
<font color="blue">b$ yields
+\begin{equation}
+\dfrac{\p\rho}{\p \mathrm{n}}\rvert_z =
+\dfrac{\p\rho}{\p\mathrm{n}}\rvert_r - \dfrac{\p\rho}{\p z}
+\dfrac{\p z}{\p\mathrm{n}}\rvert_r.
+\label{eq:30}
+\end{equation}
+The discretization of the RHS of \eqref{eq:30} are as follow,
+\begin{align*}
+\dfrac{\p\rho}{\p n}\rvert_r &=
+\dfrac{\rho(\textrm{cell2}) - \rho(\textrm{cell1})}{\textrm{dcEdge}},\\
+\dfrac{\p z}{\p n}\rvert_r &=
+\dfrac{z(\textrm{cell2}) - z(\textrm{cell1})}{\textrm{dcEdge}}.
+\end{align*}
+$\p\rho/\p z$ has been calculated before at the cell centers and on
+the layer interfaces. The data first needs to be interpolated to the cell
+edges on the layer interfaces using area weighted averaging, and then
+be interpolated to the cell edges in the middle of the layer
+interfaces.
+
+
+%-----------------------------------------------------------------------
+
+\chapter{Testing}\label{cha:testing}
The goal of the testing process is to ensure that the implementation
is working as it is intended by the design laid out in previous
sections. We particularly point out that the testing is not about
</font>
</pre>