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

\newpage
 
\begin{center}
\Large{\bf VI. RADIATING FLUID ENERGY}\\
\addcontentsline{toc}{section}{\protect\numberline{VI.}{\bf ~Radiating Fluid Energy}}
\end{center}   
 
\begin{flushleft}
\vspace*{3mm}
\large{\bf A. Differential Equation}\\
\addcontentsline{toc}{subsection}{\protect\numberline{VI.A.}{ ~Differential Equation}}

\vspace*{2mm}
\[ \frac{d}{dt}\left[\left(e +\frac{E}{\rho}\right)\rho\Delta V\right] - \Delta\left[\frac{dm}{dt}\left(e +\frac{E}{\rho}\right) -r^{\mu}F\right] \hspace{44mm} \]\\
\[ \hspace{12mm}+(p+P)\Delta(r^{\mu}u) + (E-3P)\frac{u}{r}\Delta V = \Delta\left(r^{\mu}\sigma_{e}\frac{\Delta e}{\Delta r}\right) + \epsilon_{Q}\Delta V \hspace{5mm} \rm (FE1) \]\\

\vspace*{2mm}
\rm where $\mu =$ 0 or 2, and\\

\vspace*{2mm}
\[ \epsilon_{Q} \equiv \frac{4}{3}\mu_{Q}\rho\left(\frac{\Delta u}{\Delta r} - \frac{\mu}{2} <\frac{u}{r}> \right)^{2} \hspace{60mm} \rm (FE2) \]\\

\vspace*{2mm}
\rm Note that $\epsilon_{Q}$ is cell centered, as is $\mu_{Q}$.\\

\vspace*{5mm}
\large{\bf B. Difference Equation}\\
\addcontentsline{toc}{subsection}{\protect\numberline{VI.B.}{ ~Difference Equation}}

\vspace*{2mm}
\rm For $(k = 2, \ldots, N-1),$
 
\vspace*{2mm}
$\hspace{3mm}\left[ (\rho_{k}^{n+1}e_{k}^{n+1}+E^{n+1}_{k})\Delta V^{n+1}_{k}-(\rho_{k}^{n}e_{k}^{n}+E^{n}_{k})\Delta V^{n}_{k} \right]/ dt \hfill $\\
  
\vspace*{4mm}
$ -[(m_{k+1}^{n+1}-m_{k+1}^{n})(\overline{e+\frac{E}{\rho}})_{k+1}
-(m_{k}^{n+1}-m_{k}^{n})(\overline{e+\frac{E}{\rho}})_{k}] /dt $\\
 
\vspace*{4mm}
$+(r_{k+1}^{n+\theta})^{\mu}F^{n+\theta}_{k+1} -(r_{k}^{n+\theta})^{\mu}F^{n+\theta}_{k}$\\
  
\vspace*{4mm}
$+(p^{n+\theta}_{k}+f^{n+\theta}_{k}E^{n+\theta}_{k})[(r_{k+1}^{n+\theta})^{\mu}u_{k+1}^{n+\theta}-(r_{k}^{n+\theta})^{\mu}u_{k}^{n+\theta}]$\\
  
\vspace*{2mm}
\[+\frac{\mu}{4}(1-3f_{k}^{n+\theta})E^{n+\theta}_{k}\left( \frac{u^{n+\theta}_{k+1}}{r^{n+\theta}_{k+1}} + \frac{u^{n+\theta}_{k}}{r^{n+\theta}_{k}} \right)\Delta V^{n+\theta}_{k} \hspace{50mm} \]\\

\vspace*{2mm}
\[-\sigma_{e}\left[(r^{n+\theta}_{k+1})^{\mu}(\rho^{n+\theta}_{k}+\rho^{n+\theta}_{k+1})\left(\frac{e^{n+\theta}_{k+1}-e^{n+\theta}_{k}}{r^{n+\theta}_{k+2}-r^{n+\theta}_{k}}\right)\right. \hspace{51mm}\]\\
\[\hspace{54mm} \left. -(r^{n+\theta}_{k})^{\mu}(\rho^{n+\theta}_{k-1}+\rho^{n+\theta}_{k})\left(\frac{e^{n+\theta}_{k}-e^{n+\theta}_{k-1}}{r^{n+\theta}_{k+1}-r^{n+\theta}_{k-1}}\right)\right] \]\\

\vspace*{2mm}
\[ -\frac{4}{3}(\mu_{Q})^{n+\theta}_{k}\rho^{n+\theta}_{k}
\left[ \frac{u^{n+\theta}_{k+1} - u^{n+\theta}_{k}}{r^{n+\theta}_{k+1} - r^{n+\theta}_{k}}-\frac{\mu}{2} \left(
\frac{u^{n+\theta}_{k+1}}{r^{n+\theta}_{k+1}} + \frac{u^{n+\theta}_{k}}{r^{n+\theta}_{k}} \right) \right]^{2}\Delta V^{n+\theta}_{k} \hspace{20mm} \]\\

\vspace*{2mm}
$=0$ \hfill \rm (FE3)\\
\end{flushleft}

\vspace*{2mm}
\noindent Equation (FE3) provides $N-2$ relations connecting the the energy density of the radiating fluid at $N$ gridpoints. We thus require two boundary conditions.
 
\begin{flushleft}
\vspace*{5mm}
\large{\bf C. Linearization}\\
\addcontentsline{toc}{subsection}{\protect\numberline{VI.C.}{ ~Linearization}}

\vspace*{2mm}
\rm Calculating matrix elements we find:\\
\end{flushleft}

\vspace*{2mm}
\begin{center}
\it(1) Time Derivative
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.1.}{\sl ~Time Derivative}}
\end{center}   

\begin{flushleft}
\vspace*{2mm}
{\tt e00(it,jr):}$-(\rho_{k}^{n+1}e_{k}^{n+1}+E_{k}^{n+1})\,${\tt rmup1n(k\hspace{4mm})}$/dt \hfill \rm (FE4)$\\

\vspace*{2mm}
{\tt ep1(it,jr):}$\hspace{3mm}(\rho_{k}^{n+1}e_{k}^{n+1}+E_{k}^{n+1})\,${\tt rmup1n(k+1)}$/dt \hfill \rm (FE5)$\\
 
\vspace*{2mm}
{\tt e00(it,jd):}$\hspace{5mm}\rho_{k}^{n+1}e_{k}^{n+1}[1+(\partial ln e/\partial ln \rho)_{k}^{n+1}]\,${\tt dvoln(k)}$/dt \hfill \rm (FE6)$\\

\vspace*{2mm}
{\tt e00(it,jt):}$\hspace{5mm}\rho_{k}^{n+1}e_{k}^{n+1}\hspace{8mm}(\partial ln e/\partial ln T)_{k}^{n+1}\hspace{2mm}${\tt dvoln(k)}$/dt \hfill \rm (FE7)$\\

\vspace*{2mm}
{\tt e00(it,je):}$\hspace{5mm}E_{k}^{n+1}${\tt dvoln(k)}$/dt \hfill \rm (FE8)$\\
\end{flushleft}
 
\vspace*{2mm}
\begin{center}
\it(2) Advection
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.2.}{\sl ~Advection}}
\end{center}   
 
\vspace*{2mm}
\noindent The quantities needed to calculate the derivatives of the advection
term are generated in the subroutine {\tt advectc}. As before, denote the
advected quantity as $q$. The inputs required by the subroutine for $(k=0,$ $\ldots,$ $N+2)$ are:

\[ \begin{array}{lll}
{\tt q\ \ \ (k)}&=\hspace{3mm} q_{k}^{n+1} &\equiv (e+\frac{E}{\rho})^{n+1}_{k} \\
{\tt qo\ \ (k)}&=\hspace{3mm} q_{k}^{n} &\equiv (e+\frac{E}{\rho})^{n}_{k} \\
{\tt qso\ (k)}&= Dq_{k}^{n} & = \rm monotonized\ slope\ at\ old\ time\\
{\tt flow(k)}&= -{\tt dmdt(k)}&\rm \hspace{4mm} direction\ of\ flow\ at\ interface\ k\\
\end{array} \] 

\begin{flushleft}
\rm Then, proceeding as in (C20) - (C36), and using the facts that\\ 

\[ \left( \frac{\partial ln q}{\partial lnE} \right)^{n+1}_{k} = \left(\frac{E}{\rho}\right)^{n+1}_{k}/\left(e+\frac{E}{\rho}\right)^{n+1}_{k} ={\tt dlqdle(k)} \hspace{28mm} \rm (FE9) \]

\vspace*{2mm}
\[ \left( \frac{\partial ln q}{\partial lnT} \right)^{n+1}_{k} = \left(e\frac{\partial ln e}{\partial ln T}\right)^{n+1}_{k}/\left(e+\frac{E}{\rho}\right)^{n+1}_{k} ={\tt dlqdlt(k)} \hspace{20mm}\rm (FE10)\]

\vspace*{2mm}
\[ \left( \frac{\partial ln q}{\partial ln\rho} \right)^{n+1}_{k} = \left(e\frac{\partial ln e}{\partial ln \rho}-\frac{E}{\rho}\right)^{n+1}_{k}/\left(e+\frac{E}{\rho}\right)^{n+1}_{k} ={\tt dlqdld(k)} \hspace{12mm} \rm (FE11) \]

\vspace*{2mm}
\rm we get\\

\vspace*{2mm}
\tt dqbdlem2(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln E^{n+1}_{k-2} =\ $\tt dqbdlqm2(k) dlqdle(k-2) \hfill \rm (FE12)\\

\vspace*{2mm}
\tt dqbdlem1(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln E^{n+1}_{k-1} =\ $\tt dqbdlqm1(k) dlqdle(k-1) \hfill \rm (FE13)\\

\vspace*{2mm}
\tt dqbdle00(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln E^{n+1}_{k} =\ $\tt dqbdlq00(k) dlqdle(k\hspace{4mm}) \hfill \rm (FE14)\\

\vspace*{2mm}
\tt dqbdlep1(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln E^{n+1}_{k+1} =\ $\tt dqbdlqp1(k) dlqdle(k+1) \hfill \rm (FE15)\

\vspace*{2mm}
\tt dqbdltm2(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln T^{n+1}_{k-2} =\ $\tt dqbdlqm2(k) dlqdlt(k-2) \hfill \rm (FE16)\\

\vspace*{2mm}
\tt dqbdltm1(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln T^{n+1}_{k-1} =\ $\tt dqbdlqm1(k) dlqdlt(k-1) \hfill \rm (FE17)\\

\vspace*{2mm}
\tt dqbdlt00(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln T^{n+1}_{k} =\ $\tt dqbdlq00(k) dlqdlt(k\hspace{4mm}) \hfill \rm (FE18)\\

\vspace*{2mm}
\tt dqbdltp1(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln T^{n+1}_{k+1} =\ $\tt dqbdlqp1(k) dlqdlt(k+1) \hfill \rm (FE19)\\

\vspace*{2mm}
\tt dqbdldm2(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln \rho^{n+1}_{k-2} =\ $\tt dqbdlqm2(k) dlqdld(k-2) \hfill \rm (FE20)\\

\vspace*{2mm}
\tt dqbdldm1(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln \rho^{n+1}_{k-1} =\ $\tt dqbdlqm1(k) dlqdld(k-1) \hfill \rm (FE21)\\

\vspace*{2mm}
\tt dqbdld00(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln \rho^{n+1}_{k} =\ $\tt dqbdlq00(k) dlqdld(k\hspace{4mm}) \hfill \rm (FE22)\\

\vspace*{2mm}
\tt dqbdldp1(k)$\equiv \partial(\delta\overline{q}_{k})/\partial ln \rho^{n+1}_{k+1} =\ $\tt dqbdlqp1(k) dlqdld(k+1) \hfill \rm (FE23)\\
 
\vspace*{2mm}
\rm Thus\\
 
\vspace*{2mm}
\tt e00(it,jm):$\hspace{3mm}\overline{q}_{k\hspace{4mm}}m_{k\hspace{4mm}}^{n+1}/dt$ \hfill \rm(FE24)
 
\vspace*{2mm}
\tt ep1(it,jm):$-\overline{q}_{k+1}m_{k+1}^{n+1}/dt$ \hfill \rm(FE25)
 
\vspace*{2mm}
\tt em2(it,jd):dmdt(k)dqbdldm2(k)\hfill\rm(FE26)\\
  
\vspace*{2mm}
\tt em1(it,jd):dmdt(k)dqbdldm1(k)$-$dmdt(k+1)dqbdldm2(k+1)\hfill\rm(FE27)\\
  
\vspace*{2mm}
\tt e00(it,jd):dmdt(k)dqbdld00(k)$-$dmdt(k+1)dqbdldm1(k+1)\hfill\rm(FE28)\\
  
\vspace*{2mm}
\tt ep1(it,jd):dmdt(k)dqbdldp1(k)$-$dmdt(k+1)dqbdld00(k+1)\hfill\rm(FE29)\\
  
\vspace*{2mm}
\tt ep2(it,jd):\hspace{36mm}$-$dmdt(k+1)dqbdldp1(k+1)\hfill\rm(FE30)\\
  
\vspace*{2mm}
\tt em2(it,jt):dmdt(k)dqbdltm2(k)\hfill\rm(FE31)\\
  
\vspace*{2mm}
\tt em1(it,jt):dmdt(k)dqbdltm1(k)$-$dmdt(k+1)dqbdltm2(k+1)\hfill\rm(FE32)\\
  
\vspace*{2mm}
\tt e00(it,jt):dmdt(k)dqbdlt00(k)$-$dmdt(k+1)dqbdltm1(k+1)\hfill\rm(FE33)\\
  
\vspace*{2mm}
\tt ep1(it,jt):dmdt(k)dqbdltp1(k)$-$dmdt(k+1)dqbdlt00(k+1)\hfill\rm(FE34)\\
  
\vspace*{2mm}
\tt ep2(it,jt):\hspace{36mm}$-$dmdt(k+1)dqbdltp1(k+1)\hfill\rm(FE35)\\
  
\vspace*{2mm}
\tt em2(it,je):dmdt(k)dqbdlem2(k)\hfill\rm(FE36)\\
  
\vspace*{2mm}
\tt em1(it,je):dmdt(k)dqbdlem1(k)$-$dmdt(k+1)dqbdetm2(k+1)\hfill\rm(FE37)\\
  
\vspace*{2mm}
\tt e00(it,je):dmdt(k)dqbdle00(k)$-$dmdt(k+1)dqbdlem1(k+1)\hfill\rm(FE38)\\
  
\vspace*{2mm}
\tt ep1(it,je):dmdt(k)dqbdlep1(k)$-$dmdt(k+1)dqbdle00(k+1)\hfill\rm(FE39)\\
  
\vspace*{2mm}
\tt ep2(it,je):\hspace{36mm}$-$dmdt(k+1)dqbdlep1(k+1)\hfill\rm(FE40)\\
\end{flushleft}

\begin{center}
\vspace*{2mm}
\it(3) Flux Divergence
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.3.}{\sl ~Flux Divergence}}
\end{center}   

\noindent \rm These matrix elements are the same as those given for the radiation energy equation, (RE6) -- (RE9), with the index {\tt ie} replaced by the index {\tt it}.

\vspace*{2mm}
\begin{center}
\it(4) Work
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.4.}{\sl ~Work}}
\end{center}   

\begin{flushleft}
\vspace*{2mm}
\tt e00(it,jr):$-\theta\mu(p_{k}^{n+\theta}+f_{k}^{n+\theta}E_{k}^{n+\theta})(r_{k}^{n+1}/r_{k}^{n+\theta})$\,\tt rmu(k\hspace{4mm})$u_{k}^{n+\theta}$ \hfill \rm(FE41)\\

\vspace*{2mm}
\tt ep1(it,jr):$\hspace{3mm}\theta\mu(p_{k}^{n+\theta}+f_{k}^{n+\theta}E_{k}^{n+\theta})(r_{k+1}^{n+1}/r_{k+1}^{n+\theta})$\,\tt rmu(k+1)$u_{k+1}^{n+\theta}$ \hfill \rm(FE42)\\

\vspace*{2mm}
\tt e00(it,jd):$\hspace{3mm}\theta p_{k}^{n+1}(\partial ln p/\partial ln\rho)_{k}^{n+1}\,[$\tt rmu(k+1)$u_{k+1}^{n+\theta}-$\tt rmu(k)$u_{k}^{n+\theta}$] \hfill \rm(FE43)\\

\vspace*{2mm}
\tt e00(it,ju):$-\theta(p_{k}^{n+\theta}+f_{k}^{n+\theta}E_{k}^{n+\theta})$\,rmu(k\hspace{4mm})unom(k\hspace{4mm}) \hfill \rm(FE44)\\

\vspace*{2mm}
\tt ep1(it,ju):$\hspace{3mm}\theta(p_{k}^{n+\theta}+f_{k}^{n+\theta}E_{k}^{n+\theta})$\,rmu(k+1)unom(k+1) \hfill \rm(FE45)\\

\vspace*{2mm}
\tt e00(it,jt):$\hspace{3mm}\theta p_{k}^{n+1}(\partial ln p/\partial ln T)_{k}^{n+1}[$\tt rmu(k+1)$u_{k+1}^{n+\theta}-$\tt rmu(k)$u_{k}^{n+\theta}$] \hfill \rm(FE46)\\

\vspace*{2mm}
\tt e00(it,je):$\hspace{3mm}\theta f_{k}^{n+\theta}E_{k}^{n+1}\,[$\tt rmu(k+1)$u_{k+1}^{n+\theta}-$\tt rmu(k)$u_{k}^{n+\theta}$] \hfill \rm(FE47)\\
\end{flushleft}

\newpage
\begin{center}
\vspace*{2mm}
\it(5) Anisotropy
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.5.}{\sl ~Anisotropy}}
\end{center}   

\noindent These matrix elements are the same as those given for the radiation
energy equation, (RE15) -- (RE19), with the index {\tt ie} replaced by the 
index {\tt it}.

\vspace*{2mm}
\begin{center}
\it(6) Diffusion 
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.6.}{\sl ~Diffusion}}
\end{center}   

\noindent \rm The quantities needed to calculate the derivatives of the energy
diffusion term are generated in subroutine {\tt diffuse} in the same way as
the mass diffusion in the continuity equation ({\it q.v.}). In this case
the inputs required are {\tt qd(k)} $=$ $e_{k}^{n+\theta}$ and {\tt qdn(k)}
$=$ $e_{k}^{n+1}$, for $(k = 1,$ $\ldots,$ $N+1)$, and also $\sigma$ $=$ $\sigma_{e}$, and $\zeta$ $=$ $0$. Then, as in the continuity equation,

\begin{flushleft}
\tt em1(it,jr): ddfdlrm1(k) \hfill \rm (FE48)\\
\tt e00(it,jr): ddfdlr00(k) $-$ ddfdlrm1(k+1) \hfill \rm (FE49)\\
\tt ep1(it,jr): ddfdlrp1(k) $-$ ddfdlr00(k+1) \hfill \rm (FE50)\\
\tt ep2(it,jr): \hspace{20mm} $-$ ddfdlrp1(k+1) \hfill \rm (FE51)\\

\vspace*{2mm}
\tt em1(it,jd): ddfdlqm1(k) \hfill \rm (FE52)\\
\tt e00(it,jd): ddfdlq00(k) $-$ ddfdlqm1(k+1) \hfill \rm (FE53)\\
\tt ep1(it,jd): \hspace{20mm} $-$ ddfdlq00(k+1) \hfill \rm (FE54)\\

\vspace*{2mm}
\tt em1(it,jd): ddfdltm1(k) \hfill \rm (FE55)\\
\tt e00(it,jd): ddfdlt00(k) $-$ ddfdltm1(k+1) \hfill \rm (FE56)\\
\tt ep1(it,jd): \hspace{20mm} $-$ ddfdlt00(k+1) \hfill \rm (FE57)\\

\vspace*{2mm}
\rm Here we have defined\\

\vspace*{2mm}
\tt ddfdltm1(k)$\equiv \partial Q_{k}^{n+\theta}/\partial ln T_{k-1}^{n+1}=\,${\tt ddfdlqm1(k)}$(\partial ln e/\partial ln T)_{k-1}^{n+1}$ \hfill \rm (FE58)\\

\vspace*{2mm}
\tt ddfdlt00(k)$\equiv \partial Q_{k}^{n+\theta}/\partial ln T_{k}^{n+1}=\,${\tt ddfdlq00(k)}$(\partial ln e/\partial ln T)_{k}^{n+1}$ \hfill \rm (FE59)\\
\end{flushleft}

\begin{center}
\vspace*{2mm}
\it(7) Viscous Heating
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.7.}{\sl ~Viscous Heating}}
\end{center}   

\noindent The quantities needed to calculate the viscous energy dissipation
rate {\tt qe(k)} $=$ $(\epsilon_{Q})_{k}$ and its derivatives are generated in
subroutine {\tt viscous}. It is a cell-centered quantity needed only within
the computational domain. Thus at $(k = 1,$ $\ldots,$ $N+1)$ we need the radius
$r$, velocity $u$, and some auxiliary vectors defined below, and we apply the
algorithm only at $(k = 2,$ $\ldots,$ $N-1).$

\begin{flushleft}
\rm Take $f1$ to be a discrete representation of $\frac{du}{dr} - \frac{\mu}{2}<\frac{u}{r}>$. Then define
$$f1\hspace{5mm}(r_{0},r_{+},u_{0},u_{+})=\frac{u_{+}-u_{0}}{r_{+}-r_{0}} -\frac{\mu}{4}(\frac{u_{+}}{r_{+}} + \frac{u_{0}}{r_{0}}) \hspace{38mm} \rm(FE60)$$

$$f1u0(r_{0},r_{+},u_{0},u_{+})\equiv \frac{\partial f_{1}}{\partial u_{0}}=-\frac{1}{r_{+}-r_{0}} -\frac{\mu}{4}\frac{1}{r_{0}} \hspace{38mm} \rm(FE61)$$

$$f1up(r_{0},r_{+},u_{0},u_{+})\equiv \frac{\partial f_{1}}{\partial u_{+}}=\hspace{3mm}\frac{1}{r_{+}-r_{0}} -\frac{\mu}{4}\frac{1}{r_{+}} \hspace{35mm} \rm(FE62)$$

$$f1r0(r_{0},r_{+},u_{0},u_{+})\equiv\frac{\partial f_{1}}{\partial r_{0}}=\hspace{3mm}\frac{(u_{+}-u_{0})}{(r_{+}-r_{0})^{2}} +\frac{\mu}{4}\frac{u_{0}}{r_{0}^{2}} \hspace{32mm} \rm(FE63)$$

$$f1rp(r_{0},r_{+},u_{0},u_{+})\equiv\frac{\partial f_{1}}{\partial r_{+}}=-\frac{(u_{+}-u_{0})}{(r_{+}-r_{0})^{2}} +\frac{\mu}{4}\frac{u_{+}}{r_{+}^{2}} \hspace{30mm} \rm(FE64)$$

\vspace*{2mm}
\rm Then
$${\tt dudr(k)}\equiv f1(r_{k}^{n+\theta}, r_{k+1}^{n+\theta}, u_{k}^{n+\theta}, u_{k+1}^{n+\theta}) \hspace{52mm} \rm(FE65)$$

$${\tt dudlr00(k)}\equiv \frac{\partial f_{1}}{\partial ln r_{0}}= \theta r_{k}^{n+1} f1r0(r_{k}^{n+\theta}, r_{k+1}^{n+\theta}, u_{k}^{n+\theta}, u_{k+1}^{n+\theta}) \hspace{17mm} \rm(FE66)$$

$${\tt dudlrp1(k)}\equiv \frac{\partial f_{1}}{\partial ln r_{+}}= \theta r_{k+1}^{n+1} f1rp(r_{k}^{n+\theta}, r_{k+1}^{n+\theta}, u_{k}^{n+\theta}, u_{k+1}^{n+\theta}) \hspace{17mm} \rm(FE67)$$

$${\tt dudlu00(k)}\equiv \frac{\partial f_{1}}{\partial ln u_{0}}= \theta\,{\tt unom} (k\hspace{4mm}) f1u0(r_{k}^{n+\theta}, r_{k+1}^{n+\theta}, u_{k}^{n+\theta}, u_{k+1}^{n+\theta}) \hspace{7mm} \rm(FE68)$$

$${\tt dudlup1(k)}\equiv \frac{\partial f_{1}}{\partial ln u_{+}}= \theta\,{\tt unom} (k+1) f1up(r_{k}^{n+\theta}, r_{k+1}^{n+\theta}, u_{k}^{n+\theta}, u_{k+1}^{n+\theta}) \hspace{5mm} \rm(FE69)$$
 
\vspace*{3mm}
\rm The dissipation length is

\vspace*{2mm}
\tt ql(k) $\equiv \, \ell_{0}+\frac{1}{2}(r_{k}^{n+\theta}+r_{k+1}^{n+\theta}) \, \ell_{1}$ \hfill \rm(FE70)\\

\vspace*{2mm}
\tt dqldr00(k) $\equiv (\partial ql / \partial ln r_{k}^{n+1})= \frac{1}{2}\theta r_{k}^{n+1}\, l_{1}$ \hfill \rm(FE71)\\

\vspace*{2mm}
\tt dqldrp1(k) $\equiv (\partial ql / \partial ln r_{k+1}^{n+1})= \frac{1}{2}\theta r_{k+1}^{n+1}\, l_{1}$ \hfill \rm(FE72)\\
 
\vspace*{3mm}
\rm The velocity divergence is

\vspace*{2mm}
\tt div(k) $\equiv$ [rmu(k+1)$u_{k+1}^{n+\theta}-$rmu(k)$u_{k}^{n+\theta}$] / dvol(k) \hfill \rm (FE73)\\

\vspace*{2mm}
\tt ddivdr00(k) $\equiv \partial$div(k)/$\partial r_{k}^{n+1}=$\\
\end{flushleft}
{\flushleft\hfill $[\mu\,$\tt rmum1(k\hspace{4mm})$u_{k}^{n+\theta}-$rmu(k\hspace{4mm})div(k)$]/$dvol(k) \rm (FE74)}\\

\vspace*{2mm}
{\flushleft \tt ddivdrp1(k) $\equiv \partial$div(k)/$\partial r_{k+1}^{n+1}=$}
{\flushleft\hfill $-[\mu\,$\tt rmum1(k+1)$u_{k+1}^{n+\theta}-$rmu(k+1)div(k)$]/$dvol(k) \rm (FE75)}\\

\vspace*{2mm}
{\flushleft \tt ddivdu00(k)$\equiv \partial{\tt div(k)}/\partial u_{k}^{n+1}=-$rmu(k\hspace{4mm})/dvol(k) \hfill \rm (FE76)}\\

\vspace*{2mm}
{\flushleft \tt ddivdup1(k)$\equiv \partial{\tt div(k)}/\partial u_{k+1}^{n+1}=$\hspace{4mm}rmu(k+1)/dvol(k) \hfill \rm (FE77)}\\

\vspace*{3mm}
\noindent\rm Define
{\flushleft {\tt qv(k)}$\equiv$ min[{\tt div(k)}, 0] \hfill \rm (FE78)}

\vspace*{2mm}
\noindent\rm Then
{\flushleft \tt dqvdlr00(k) $=$ cvmgm[$\theta r_{k}^{n+1}$\tt ddivdr00(k), 0, \tt div(k)] \hfill \rm(FE79)}

\vspace*{2mm}
{\flushleft \tt dqvdlrp1(k) $=$ cvmgm[$\theta r_{k+1}^{n+1}$\tt ddivdrp1(k), 0, \tt div(k)] \hfill \rm(FE80)}

\vspace*{2mm}
{\flushleft \tt dqvdlu00(k) $=$ cvmgm[$\theta\,unom$(k\hspace{4mm})ddivdu00(k), 0, \tt div(k)] \hfill \rm(FE81)}

\vspace*{2mm}
{\flushleft \tt dqvdlup1(k) $=$ cvmgm[$\theta\,unom$(k+1)ddivdru1(k), 0, \tt div(k)] \hfill \rm(FE82)}\\

\vspace*{2mm}
\noindent\rm The coefficient of viscosity is
{\flushleft \tt qm(k) $\equiv (\mu_{Q})^{n+\theta}_{k} = C_{1}$ql(k)$(a_{s})^{n+\theta}_{k} - C_{2}$[ql(k)]$^{2}$qv(k) \hfill \rm (FE83)}

\vspace*{2mm}
{\flushleft \tt dqmdlr00(k) $\equiv \partial(\mu_{Q})^{n+\theta}_{k} / \partial ln r^{n+1}_{k} =$}\\
{\flushleft\hfill \tt dqldr00(k)$\{C_{1}a_{s,k}^{n+\theta}-C_{2}$[2ql(k)qv(k)+ql(k)$^{2}$dqvdlr00(k)]\} \rm(FE84)}\\

\vspace*{2mm}
{\flushleft \tt dqmdlrp1(k) $\equiv \partial(\mu_{Q})^{n+\theta}_{k} / \partial ln r^{n+1}_{k+1} =$}\\
{\flushleft\hfill \tt dqldrp1(k)$\{C_{1}a_{s,k}^{n+\theta}-C_{2}$[2ql(k)qv(k)+ql(k)$^{2}$dqvdlrp1(k)]\} \rm(FE85)}\\

\vspace*{2mm}
{\flushleft \tt dqmdld00(k)$\equiv \partial(\mu_{Q})^{n+\theta}_{k}/\partial ln \rho^{n+1}_{k} =$}\\
{\flushleft\hfill $\frac{1}{2}\theta C_{1}${\tt ql(k)}$a_{s,k}^{n+\theta}\left[(p_{k}^{n+1}/p_{k}^{n+\theta})\left(\frac{\partial ln p}{\partial ln T}\right)_{k}^{n+1}-(\rho_{k}^{n+1}/\rho^{n+\theta}_{k}) \right]$\hspace{2mm}\rm(FE86)}\\

\vspace*{2mm}
{\flushleft \tt dqmdlt00(k)$\equiv \partial(\mu_{Q})^{n+\theta}_{k}/\partial ln T^{n+1}_{k} =$}\\
{\flushleft\hfill $\frac{1}{2}\theta C_{1}${\tt ql(k)}$a_{s,k}^{n+\theta}(p_{k}^{n+1}/p_{k}^{n+\theta})\left(\frac{\partial ln p}{\partial ln T}\right)_{k}^{n+1}$\hspace{2mm}\rm(FE87)}\\

\vspace*{2mm}
{\flushleft \tt dqmdlu00(k)$\equiv \partial(\mu_{Q})^{n+\theta}_{k} / \partial ln u^{n+1}_{k}=-C_{2}$\tt ql(k)$^{2}$\tt dqvdlu00(k) \hfill \rm(FE88)}\\

\vspace*{2mm}
{\flushleft \tt dqmdlup1(k)$\equiv \partial(\mu_{Q})^{n+\theta}_{k} / \partial ln u^{n+1}_{k+1}=-C_{2}$\tt ql(k)$^{2}$\tt dqvdlup1(k) \hfill \rm(FE89)}\\

\vspace*{2mm}
\noindent{The rate of viscous energy dissipation is (cf. FE2 and GM48)}
{\flushleft \tt qe(k)\,=\,-$\frac{4}{3}\rho_{k}^{n+\theta}$qm(k)[dudr(k)]$^{2}$dvol(k)$\,\equiv\,$\tt qf(k)dudr(k)dvol(k)\hfill\rm(FE90) }\\

\vspace*{2mm}
\noindent\rm Then
{\flushleft \tt e00(it,jr):dqfdlr00(k)dudr(k)dvol(k) + qf(k)durlu00(k)dvol(k)}\\
{\flushleft\hfill \tt $-\theta\,$qf(k)dudr(k)rmu(k\hspace{4mm})$r_{k}^{n+1}$ \rm(FE91)}\\

\vspace*{2mm}
{\flushleft \tt ep1(it,jr):dqfdlrp1(k)dudr(k)dvol(k) + qf(k)durlup1(k)dvol(k)}\\
{\flushleft\hfill \tt $-\theta\,$qf(k)dudr(k)rmu(k+1)$r_{k+1}^{n+1}$ \rm(FE92)}\\

\vspace*{2mm}
{\flushleft \tt e00(it,ju):}
{\flushleft\hfill \tt dqfdlu00(k)dudr(k)dvol(k) + qf(k)durlu00(k)dvol(k) \rm(FE93)}\\

\vspace*{2mm}
{\flushleft \tt ep1(it,ju):}
{\flushleft\hfill \tt dqfdlup1(k)dudr(k)dvol(k) + qf(k)durlup1(k)dvol(k) \rm(FE94)}\\

\vspace*{2mm}
{\flushleft \tt e00(it,jt):dqfdlt00(k)dudr(k)dvol(k) \hfill \rm(FE95)}\\

\vspace*{2mm}
{\flushleft \tt e00(it,jd):dqfdld00(k)dudr(k)dvol(k) \hfill \rm(FE96)}\\

\newpage
\begin{center}
\it(8) Right Hand Side
\addcontentsline{toc}{subsubsection}{\protect\numberline{VI.C.8.}{\sl ~Right Hand Side}}
\end{center}   

\vspace*{2mm}
{\flushleft $-$\ \tt rhs(it)$=$}
\vspace*{2mm}
{\flushleft \hspace{3mm}$\left[(\rho e + E)_{k}^{n+1}{\tt dvoln(k)} - (\rho e+E)_{k}^{n}{\tt dvolo(k)}\right]/ dt \hfill $}\\
  
\vspace*{2mm}
{\flushleft $-$\ \tt dmdt(k$+$1)$(\overline{e+\frac{E}{\rho}})_{k+1}\ +$ dmdt(k)$(\overline{e+\frac{E}{\rho}})_{k}$ }\\
 
\vspace*{2mm}
{\flushleft $+\ ${\tt rmu(k$+$1)}$F^{n+\theta}_{k+1}\,-\,{\tt rmu(k)}F^{n+\theta}_{k}$ }\\
  
\vspace*{2mm}
{\flushleft $+\ (p + fE)^{n+\theta}_{k}[$\tt rmu(k$+$1)$u_{k+1}^{n+\theta}\,-\,$rmu(k)$u_{k}^{n+\theta}] $}\\
  
\vspace*{2mm}
{\flushleft $+\ \frac{\mu}{4}(1-3f_{k}^{n+\theta})E^{n+\theta}_{k}\left[ (u^{n+\theta}_{k+1} / r^{n+\theta}_{k+1}) + (u^{n+\theta}_{k} / r^{n+\theta}_{k}) \right]$\tt dvol(k) }\\

\vspace*{2mm}
{\flushleft $+$\ \tt df(k)$\,-\,$df(k+1)$\,+\,$qe(k) \hfill \rm (FE97) }\\
 
%\end{document}   
.