%\documentstyle[11pt]{article} 
%\begin{document}
%\setcounter{page}{24}

\newpage
 
\begin{center}
\Large{\bf IV. CONTINUITY}
\addcontentsline{toc}{section}{\protect\numberline{IV.}{\bf ~Continuity}}
\end{center}   

\begin{flushleft}
\vspace*{3mm}
\large{\bf A. Differential Equation}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.A.}{ ~Differential Equation}}

\vspace*{2mm}
\[ \hspace{30mm}\frac{d(\rho\Delta V)}{dt} + \Delta(\rho r^{\mu} u_{rel}) = \Delta(r^{\mu}\sigma_{\rho} \frac{\Delta\rho}{\Delta r}) \hspace{24mm}\rm (C1)\]\\
 
\vspace*{2mm}
\large{\bf B. Difference Equation}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.B.}{ ~Difference Equation}}

\vspace*{2mm}
\rm For $(k = 2, \ldots, N - 1),$\\

\vspace*{2mm}
\[ \frac{\rho_{k}^{n+1} [(r^{n+1}_{k+1})^{\mu+1} - (r^{n+1}_{k})^{\mu+1}] - 
\rho_{k}^{n} [(r^{n}_{k+1})^{\mu+1} - (r^{n}_{k})^{\mu+1}]} {(\mu+1)\ dt} \hspace{30mm} \]\\

\vspace*{2mm}
\[ +\ (r^{n+\theta}_{k+1})^{\mu} \left[ u^{n+\theta}_{k+1} - \left( \frac {r^{n+1}_{k+1} - r^{n}_{k+1}}{dt} \right) \right] \overline{\rho}_{k+1} \hspace{55mm} \]\\

\[ -\ (r^{n+\theta}_{k})^{\mu} \left[ u^{n+\theta}_{k} - \left( \frac {r^{n+1}_{k}- \ \ r^{n}_{k}}{dt} \right) \right] \overline{\rho}_{k} \hspace{60mm} \]\\

\vspace*{2mm}
\[ -2\sigma_{\rho} \left[ (r_{k+1}^{n+\theta})^{\mu} \left(\frac{\rho^{n+\theta}_{k+1} - \rho^{n+\theta}_{k}}{r^{n+\theta}_{k+2} - r^{n+\theta}_{k}} \right) 
- (r_{k}^{n+\theta})^{\mu} \left(\frac{\rho^{n+\theta}_{k} - \rho^{n+\theta}_{k-1}}{r^{n+\theta}_{k+1} - r^{n+\theta}_{k-1}} \right) \right] = 0 \hspace{9mm} \rm (C2)\]\\

\vspace*{2mm}
\rm where $\mu =$ 0, 1, or 2.\\
\end{flushleft}

\noindent\rm Equation (C2) provides $N-2$ relations connecting the densities
at N gridpoints. We thus require two boundary conditions. In equation (C2) we
have written time-centered values of the physical variables as:\\

\vspace*{2mm}
\hspace{35mm}$x^{n+\theta} \equiv \theta x^{n+1} + (1 - \theta) x^{n}$ \hfill \rm (C3)\\

\begin{flushleft}
\vspace*{2mm}
\rm Values for the density advected at cell interfaces are defined as\\
\newpage

\vspace*{2mm}
$\overline{\rho}_{k} \equiv (0.5 + s_{k})[\theta(\rho^{n+1}_{k-1} + 0.5 D\rho^{n+1}_{k-1}) + (1-\theta)(\rho_{k-1}^{n} + 0.5D\rho_{k-1}^{n})]$\\
\hspace{5mm}$ + \ (0.5 - s_{k})[\theta(\rho^{n+1}_{k} - 0.5 D\rho^{n+1}_{k}) + (1-\theta)(\rho_{k}^{n}\ \ \  - 0.5D\rho_{k}^{n}\ \ \ )]$ \hfill \rm (C4)\\

\vspace*{2mm}
\rm The switch $s_{k}$\\

\vspace*{2mm}
\hspace{35mm}$s_{k} \equiv {\tt cvmgp}[0.5, -0.5, u_{rel,k}]$ \hspace{35mm} \rm (C5)\\
\end{flushleft}

\vspace*{2mm}
\noindent\rm chooses the upstream direction, and the slope $D\rho_{k}^{n+1}$ inside the cell is given by

\begin{flushleft}
\vspace*{2mm}
\[ \hspace{20mm} D\rho_{k}^{n+1} \equiv {\tt cvmgp} \left[ \frac{C \Delta\rho_{k-1}^{n+1}\Delta\rho_{k}^{n+1}}{\Delta\rho^{n+1}_{k-1} + \Delta\rho^{n+1}_{k}},\ 0,\ \Delta\rho^{n+1}_{k} \right] \hspace{20mm} \rm (C6) \]\\
\end{flushleft}

\vspace*{2mm}
\noindent\rm which yields monotonized van Leer advection if $C = 2$, and donor cell advection if $C = 0$. Here

\begin{flushleft}
\vspace*{2mm}
\hspace{45mm} $\Delta\rho_{k+1}^{n+1} \equiv \rho^{n+1}_{k+1} - \rho^{n+1}_{k}$ \hfill \rm (C7)

\vspace*{5mm}
\large{\bf C. Linearization}
\addcontentsline{toc}{subsection}{\protect\numberline{IV.C.}{ ~Linearization}}
\vspace*{2mm}

\vspace*{2mm}
\rm For $(k = 1, \ldots, N + 1)$ define\\

\vspace*{2mm}
\tt rmuo(k) $\equiv (r_{k}^{n})^{\mu}$ \hfill \rm (C8) \\

\vspace*{2mm}
\tt rmu\ (k) $\equiv (r_{k}^{n+\theta})^{\mu}$ \hfill \rm (C9) \\

\vspace*{2mm}
\tt rmun(k) $\equiv (r_{k}^{n+1})^{\mu}$ \hfill \rm (C10) \\

\vspace*{2mm}
\tt urel(k) $\equiv u_{k}^{n+\theta} - (r^{n+1}_{k} -r^{n}_{k}) \, / \, dt $ \hfill \rm (C11) \\

\vspace*{2mm}
$ a_{s,k}^{n+\theta} \hspace{7mm}\equiv \sqrt{\gamma p_{k}^{n+\theta}/\rho_{k}^{n+\theta}} $ \hfill \rm (C12) \\

\vspace*{2mm}
Calculating matrix elements we find:\\
\end{flushleft}
 
\begin{center}
\vspace*{2mm}
\it(1) Time Derivative
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.C.1.}{\sl ~Time Derivative}}
\end{center}   

\begin{flushleft}
\vspace*{2mm}
\tt e00(id,jr): $-\rho_{k}^{n+1}{\tt rmup1n(k\hspace{6mm})}\,/ \,dt \hfill \rm (C13) $\\

\vspace*{2mm}
\tt ep1(id,jr): $\hspace{3mm}\rho_{k}^{n+1}{\tt rmup1n(k+1)}\,/ \,dt \hfill \rm (C14) $\\

\vspace*{2mm}
\tt e00(id,jd): $\hspace{3mm}\rho_{k}^{n+1}{\tt dvoln\ (k\hspace{6mm})}\,/ \,dt \hfill \rm (C15) $\\
\end{flushleft}

\newpage
\vspace*{2mm}
\begin{center}
\it(2) Advection
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.C.2.}{\sl ~Advection}}
\end{center}   
 
\begin{flushleft}
\vspace*{2mm}
$${\tt e00(id,jr)}:{\tt rmu(k \hspace{8mm} )} \overline{\rho}_{k\ \ } \left[\ \ \left( \frac{r_{k}^{n+1}}{dt} \right) - \mu\theta \left( \frac{r^{n+1}_{k}}{r^{n+\theta}_{k}} \right) {\tt urel(k\hspace{8mm})}\right]\hspace{1mm}\rm (C16)$$\\
 
\vspace*{2mm}
$${\tt e00(id,jr)}:{\tt rmu(k+1)} \overline{\rho}_{k+1} \left[- \left( \frac{r_{k+1}^{n+1}}{dt} \right) + \mu\theta \left( \frac{r^{n+1}_{k+1}}{r^{n+\theta}_{k+1}} \right) {\tt urel(k+1)} \right] \hspace{1mm} \rm (C17) $$\\

\vspace*{2mm}
\tt e00(id,ju):$-\theta$ \tt rmu(k\hspace{4mm}) $\overline{\rho}_{k\hspace{4mm}}$ \tt unom(k\hspace{4mm}) \hfill \rm (C18)\\

\vspace*{2mm}
\tt ep1(id,ju):$\hspace{3mm} \theta$ \tt rmu(k+1) $\overline{\rho}_{k+1}$ \tt unom(k+1) \hfill \rm (C19)\\
\end{flushleft}

\noindent \rm The quantities needed to calculate the derivatives of the advection term with respect to density are generated in a separate subroutine, named {\tt advectc} to indicate that it treats the advection of cell-centered variables such as $\rho,$ $\ e,$ and $e+(E/\rho)$. For generality, and to avoid repetition, we shall denote the advected quantity as $q$. The inputs required by the subroutine are, for $(k = $ $-1,$ $\ldots,$ $N+2),$

\begin{flushleft}
\vspace*{2mm}
\hspace{10mm} $q\,\ \ \hspace{9mm}\equiv \ q_{k}^{n+1}$\hspace{5mm}  = advected quantity at advanced time\\

\vspace*{2mm}
\hspace{10mm} $qo\,\ \hspace{9mm}\equiv \ q_{k}^{n}$\hspace{9mm} = advected quantity at old time\\

\vspace*{2mm}
\hspace{10mm} $qso \hspace{9mm}\equiv  Dq_{k}^{n}$\hspace{7mm} = monotonized slope at old time\\

\vspace*{2mm}
\hspace{10mm} {\tt flow(k)} $\equiv$ {\tt urel(k)} $=$ direction of flow at interface $k$\\

\vspace*{2mm}
\rm Then\\

\vspace*{2mm}
\tt s(k) $\hspace{8mm}\equiv s_{k} = {\tt cvmgp}[0.5, - 0.5,{\tt flow(k)}\,] \hfill (k = 1,\ldots, N+1) \hspace{4mm} \rm (C20)$\\

\vspace*{2mm}
\tt dq(k) $\hspace{6mm}\equiv \Delta q_{k}^{n+1} = q_{k+1}^{n+1} - q^{n+1}_{k} \hfill (k = 0,\ldots, N+1) \hspace{4mm} \rm (C21)$\\

\vspace*{2mm}
\tt denom(k) $\equiv {\tt cvmgm[\,dq(k-1) + dq(k)} - \epsilon,$\\ 

\vspace*{1mm}
\hspace{31mm} ${\tt \,dq(k-1) + dq(k)} + \epsilon,$\\

\vspace*{1mm}
\hspace{30mm} ${\tt \ dq(k-1) + dq(k)\,] } \hspace{17mm} (k = 1, \ldots, N+1)$ \hfill \rm (C22)\\

\vspace*{2mm}
{\tt qr(k)} $\equiv R^{n+1}_{k} = C_{adv}{\tt dq(k-1)\,dq(k)\,/\,denom(k)} \hspace{1mm} (k = 1,\ldots, N+1) \hfill \rm (C23)$\\

\vspace*{2mm}
{\tt qsn(k)} $\equiv Dq^{n+1}_{k} = {\tt cvmgp[ qr(k), 0, dq(k-1)dq(k)]}(k = 0, \ldots, N+1) \hfill \rm (C24)$\\

\vspace*{2mm}
\rm and\\

\vspace*{2mm}
${\tt qb(k)} \equiv \overline{q}_{k}^{n+1} =$\\

\vspace*{1mm}
$\,\ (0.5 + s_{k})\{ \theta [{\tt q(k-1) + 0.5qsn(k-1)}] + (1 - \theta)[{\tt qo(k-1) + 0.5 qso(k-1)}]\}$\\

\vspace*{1mm}
$+ (0.5 - s_{k})\{ \theta [{\tt q(k\hspace{6mm}) + 0.5qs(k\hspace{7mm})}] + (1 - \theta)[{\tt qo(k\hspace{6mm}) + 0.5 qso(k\hspace{7mm})}] \}$
\end{flushleft}

\begin{flushright}
$(k = 1,\ldots, N+1)\ $ \rm (C25)\\
\end{flushright}

\begin{flushleft}
\rm Hence\\

\vspace*{2mm}
{\tt dqrdlqm1(k)}$ \equiv \partial R^{n+1}_{k}/\partial ln q^{n+1}_{k-1} = - C_{adv} q^{n+1}_{k-1}[\tt{dq(k)/denom(k)}]^{2} \hfill \rm(C26)$\\

\vspace*{2mm}
{\tt dqrdlq00(k)}$ \equiv \partial R^{n+1}_{k}/\partial ln q^{n+1}_{k} =$\\ 
\end{flushleft}

\begin{flushright}
\hfill $- C_{adv} q^{n+1}_{k} {\tt[dq(k)-dq(k-1)]/denom(k)} \hspace{5mm} \rm(C27)$\\
\end{flushright}

\begin{flushleft}
\vspace*{2mm}
{\tt dqrdlqp1(k)}$ \equiv \partial R^{n+1}_{k}/\partial ln q^{n+1}_{k+1} = + C_{adv} q^{n+1}_{k+1}[\tt{dq(k-1)/denom(k)}]^{2} \hfill \rm(C28)$\\

\vspace*{2mm}
\rm and\\

\vspace*{2mm}
{\tt dqsdlqm1(k)}$\equiv \partial(Dq^{n+1}_{k})/\partial ln q^{n+1}_{k-1} =$\\
\end{flushleft}

\begin{flushright}
$ {\tt cvmgp[ dqrdlqm1(k), 0, dq(k-1)dq(k)]} \hspace{5mm} \rm (C29)$\\
\end{flushright}

\begin{flushleft}
\vspace*{2mm}
{\tt dqsdlq00(k)}$\equiv \partial(Dq^{n+1}_{k})/\partial ln q^{n+1}_{k} =$\\
\end{flushleft}

\begin{flushright}
${\tt cvmgp[ dqrdlq00(k), 0, dq(k-1)dq(k)]} \hspace{5mm} \rm (C30)$\\
\end{flushright}

\begin{flushleft}
\vspace*{2mm}
{\tt dqsdlqp1(k)}$\equiv \partial(Dq^{n+1}_{k})/\partial ln q^{n+1}_{k+1} =$\\
\end{flushleft}

\begin{flushright}
$ {\tt cvmgp[ dqrdlqp1(k), 0, dq(k-1)dq(k)]} \hspace{5mm} \rm (C31)$\\
\end{flushright}

\begin{flushleft}
\vspace*{2mm}
\rm Then\\

$\delta\overline{q}^{n+1}_{k} = \theta[(0.5+s_{k})(\delta q^{n+1}_{k-1} +0.5 \delta Dq^{n+1}_{k-1})$\\

\hfill $+\  (0.5-s_{k})(\delta q^{n+1}_{k} -0.5 \delta Dq^{n+1}_{k})] \hspace{5mm} \rm (C32) $\\

\vspace*{2mm}
\rm Therefore\\

${\tt dqbdlqm2(k)} \equiv \partial(\delta\overline{q}^{n+1}_{k})/\partial ln q_{k-2}^{n+1} = \ \frac{1}{2}\theta (0.5+s_{k}){\tt dqsdlqm1(k-1)} \hfill \rm (C33) $\\

\newpage
\vspace*{2mm}
${\tt dqbdlqm1(k)} \equiv \partial(\delta\overline{q}^{n+1}_{k})/\partial ln q_{k-1}^{n+1} =$\\
\vspace*{1mm}
$\hfill \ \theta \{ (0.5+s_{k}){\tt [q(k-1) + \frac{1}{2} dqsdlq00(k-1)]} - \frac{1}{2}(0.5-s_{k}){\tt dqsdlqm1(k)} \} \hspace{1mm} \rm (C34)$\\

\vspace*{2mm}
${\tt dqbdlq00(k)} \equiv \partial(\delta\overline{q}^{n+1}_{k})/\partial ln q_{k}^{n+1} =$\\
\vspace*{1mm}
$\hfill \ \theta \{ \frac{1}{2} (0.5+s_{k}){\tt dqsdlqp1(k-1)} + (0.5-s_{k}){\tt [q(k) - \frac{1}{2}dqsdlq00(k)]} \}\hspace{1mm}\rm (C35)$\\

\vspace*{2mm}
${\tt dqbdlqp1(k)} \equiv \partial(\delta\overline{q}^{n+1}_{k})/\partial ln q_{k+1}^{n+1} = - \frac{1}{2}\theta (0.5-s_{k}){\tt dqsdlqp1(k)} \hfill \rm (C36)$\\
\end{flushleft}

\vspace*{2mm}
\noindent \rm For the continuity equation, the advected quantity $q_{k}$ $\equiv$ $\rho_{k},$ hence $\partial/\partial q$ $\equiv$ $\partial/\partial\rho.$ Thus

\begin{flushleft}
\vspace*{2mm}
\tt em2(id,jd): \hfill $-$ \tt rmu(k) urel(k) dqbdlqm2(k) \rm (C37)\\

\vspace*{1mm}
\tt em1(id,jd): rmu(k+1) urel(k+1) dqbdlqm2(k+1)\\
\end{flushleft}
\begin{flushright}
$-$ \tt rmu(k) urel(k) dqbdlqm1(k) \rm (C38)\\
\end{flushright}

\begin{flushleft}
\vspace*{1mm}
\tt e00(id,jd): rmu(k+1) urel(k+1) dqbdlqm1(k+1)\\
\end{flushleft}
\begin{flushright}
$-$ \tt rmu(k) urel(k) dqbdlq00(k) \rm (C39)\\
\end{flushright}

\begin{flushleft}
\vspace*{1mm}
\tt ep1(id,jd): rmu(k+1) urel(k+1) dqbdlq00(k+1)\\
\end{flushleft}
\begin{flushright}
$-$ \tt rmu(k) urel(k) dqbdlqp1(k) \rm (C40)\\
\end{flushright}

\begin{flushleft}
\vspace*{1mm}
\tt ep2(id,jd): rmu(k+1) urel(k+1) dqbdlqp1(k+1) \hfill \rm (C41)\\
\end{flushleft}

\begin{center}
\it(3) Diffusion
\addcontentsline{toc}{subsubsection}{\protect\numberline{IV.C.3.}{\sl ~Diffusion}}
\end{center}   

\noindent \rm The quantities needed to calculate the derivatives of the diffusion term are generated in a separate subroutine, named {\tt diffuse}. It treats the diffusion of both the density $\rho$ in the continuity equation, and the gas energy densi- ty $e$ in the radiating fluid energy equation. For generality, and to avoid re- petition, we shall denote the quantity diffused as $q$. The inputs required by the subroutine are, for $(k = 1,$ $\ldots,$ $N+1):$

\begin{center}
\begin{tabular}{l l l|c|c} \hline
 & & & \rm Mass & Energy\\
 & & & \rm diffusion & diffusion\\ \hline
\tt qd\ (k) $\equiv$ & $q^{n+\theta}_{k} =$ & \rm quantity diffused at $t^{n+\theta} =$ & $\rho_{k}^{n+\theta}$ & $e_{k}^{n+\theta}$ \\
\tt qdn(k) $\equiv$ & $q^{n+1}_{k} =$ & \rm quantity diffused at $t^{n+1} =$ & $\rho_{k}^{n+1}$ & $e_{k}^{n+1}$ \\
 & $\sigma\ \ \ \ =$ & diffusion coefficient \hspace{8mm} $=$ & $\sigma_{\rho}$ & $\sigma_{e}$ \\
 & {\tt zet}\ \ \ \ $=$ & density exponent \hspace{12mm} $=$ & 0 & 1 \\ \hline  
\end{tabular}
\end{center}

\begin{flushleft}
\vspace*{2mm}
\rm Then the diffusion flux {\tt df(k)} and its derivatives for $(k = 1,\ldots, N+1)$ are \\

\tt df(k) $\equiv Q^{n+\theta}_{k} =$\\
\end{flushleft}
\begin{flushright}
\hfill \tt{$2^{(1-{\tt zet})}\,\sigma\,$rmu(k)$(\rho_{k}^{n+\theta}+\rho_{k-1}^{n+\theta})^{\tt zet}($qd(k)$-$qd(k-1)$)/(r^{n+\theta}_{k+1} - r^{n+\theta}_{k-1}$) \rm (C42)}\\
\end{flushright}

\begin{flushleft}
\vspace*{2mm}
\tt ddfdlrm1(k) $\equiv \partial Q^{n+\theta}_{k}/\partial ln r^{n+1}_{k-1} = \  \theta$\,\tt{df(k)}$r_{k-1}^{n+1} /(r^{n+\theta}_{k+1} - r^{n+\theta}_{k-1})$ \hfill \rm (C43) \\

\vspace*{2mm}
\tt ddfdlr00(k) $\equiv \partial Q^{n+\theta}_{k}/\partial ln r^{n+1}_{k} = \mu \theta$\,\tt{df(k)}$(r_{k}^{n+1} /r^{n+\theta}_{k})$ \hfill \rm (C44)\\

\vspace*{2mm}
\tt ddfdlrp1(k) $\equiv \partial Q^{n+\theta}_{k}/\partial ln r^{n+1}_{k+1} = -\theta$\,\tt{df(k)}$r_{k+1}^{n+1} /(r^{n+\theta}_{k+1} - r^{n+\theta}_{k-1})$ \hfill \rm (C45)\\

\tt ddfdlqm1(k) $\equiv \partial Q^{n+\theta}_{k}/\partial ln q^{n+1}_{k-1} =$ \hfill\\
\end{flushleft}

\begin{flushright}
\hfill $ -2^{(1-{\tt zet})}\,\sigma \theta$ \tt{qdn(k-1)}\,$(\rho^{n+\theta}_{k} + \rho^{n+\theta}_{k-1})^{\tt zet}\,(r_{k}^{n+\theta})^{\mu} /(r^{n+\theta}_{k+1} - r^{n+\theta}_{k-1})$ \rm (C46)\\
\end{flushright}

{\flushleft  \tt{ddfdlq00(k)} $\equiv \partial Q^{n+\theta}_{k}/\partial ln q^{n+1}_{k} =$ \hfill}\\
{\flushright \hfill $ \ 2^{(1-{\tt zet})}\,\sigma \theta$ \tt{qdn(k\ \ )}\,$(\rho^{n+\theta}_{k} + \rho^{n+\theta}_{k-1})^{\tt zet}\,(r_{k}^{n+\theta})^{\mu} /(r^{n+\theta}_{k+1} - r^{n+\theta}_{k-1})$ \rm (C47) }\\

{\flushleft  \tt{ddfdldm1(k)} $\equiv \partial Q^{n+\theta}_{k}/\partial ln\rho^{n+1}_{k-1} = {\tt zet} \theta$\,\tt{df(k)}$\rho_{k-1}^{n+1} /(\rho^{n+\theta}_{k} + \rho^{n+\theta}_{k-1})$ \hfill \rm (C48) }\\

{\flushleft  \tt{ddfdld00(k)} $\equiv \partial Q^{n+\theta}_{k}/\partial ln\rho^{n+1}_{k} = {\tt zet} \theta$\,\tt{df(k)}$\rho_{k}^{n+1} /(\rho^{n+\theta}_{k} + \rho^{n+\theta}_{k-1})$ \hfill \rm (C49)}\\

\vspace*{2mm}
\noindent\rm Thus\\
{\flushleft \tt em1(id,jr): ddfdlrm1(k) \hfill \rm (C50)} \\
{\flushleft \tt e00(id,jr): ddfdlr00(k) $-$ ddfdlrm1(k+1) \hfill \rm (C51)} \\
{\flushleft \tt ep1(id,jr): ddfdlrp1(k) $-$ ddfdlr00(k+1) \hfill \rm (C52)} \\
{\flushleft \tt ep2(id,jr): \hspace{20mm} $-$ ddfdlrp1(k+1) \hfill \rm (C53)} \\

\vspace*{2mm}
{\flushleft \tt em1(id,jd): ddfdlqm1(k) \hfill \rm (C54)} \\
{\flushleft \tt e00(id,jd): ddfdlq00(k) $-$ ddfdlqm1(k+1) \hfill \rm (C55)} \\
{\flushleft \tt ep1(id,jd): \hspace{20mm} $-$ ddfdlq00(k+1) \hfill \rm (C56)} \\

\vspace*{2mm}
\noindent\rm Notice that in equations (C55) and (C56) we have made explicit use of the
facts that for mass diffusion, {\tt zet} $\equiv 0$, and $(\partial ln q/\partial ln \rho) \equiv 1$. Finally,\\
\vspace*{2mm}
{\flushleft $ - {\tt rhs(id)} = [\,\rho^{n+1}_{k} {\tt dvoln(k)} - \rho^{n}_{k}\,{\tt dvolo(k)}\,]\ /\ dt $}\\
{\flushleft $ \hspace{17mm} +\ {\tt rmu(k+1)\,urel(k+1)} \,\overline\rho_{k+1} - {\tt rmu(k)\,urel(k)}\ \overline\rho_{k} $}\\
{\flushleft $\hspace{17mm} -\  {\tt df(k+1)} + {\tt df(k)}$ \hfill \rm (C57)}\\

%\end{document}   
.
