\subsection{A Very Brief Description of {\TITAN} }

We discretize the physical equations (1)-(5) with finite differences. 
The physical quantities are defined on a staggered mesh. This allows us to write down the 
spatial operators in a straightforward manner, although some attention must be paid to the 
advection terms. Here we have employed the second-order, upwind Van Leer advection scheme 
[8]. The temporal evolution is achieved with a simple formula to advance from the
old to the new time level. As a consequence our variables are functions of the integer 
pair $(k,n)$, \ie~of zone number and time index. In the equations all quantities $x$, 
aside from the temporal operator, are defined in between the old and the new time step with
the aid of the time centering parameter $\theta$,
\begin{flushleft}
$\hspace{60mm} x_k = \theta x_k^{n+1} + (1-\theta) x_k^n~.  \hfill (24)
$
\end{flushleft}
The algorithm becomes fully implicit for $\theta =1$ and assumes backward-Euler
differencing for $\theta = {1\over2}$. In general we adopt ${1\over2}< \theta\le 1$.

The physical equations are formulated in conservation form using a discrete volume scheme.
We provide the initial solution of the fundamental quantities for this set of algebraic 
equations. Their linearization results in a block pentadiagonal matrix. This system is inverted
by the Henyey scheme~[5] and then the solution is advanced in time by Newton's method. 
The data on the new time level are said to be converged when the maximum of the relative 
change of the variables is less than a specified tolerance, typically $10^{-5}$, and the 
number of iterations lies within a specified range, say 3 and 15. The new time step is 
evaluated by a nonlinear control which incorporates accuracy considerations such that the 
temporal change of the variables is less than a given limit, typically 10\%. 

\subsection{The Adaptive Grid Equation}

At the core of our code is the adaptive grid algorithm. This couples the physical 
equations with a grid equation nonlinearly in a twofold manner. One way is the 
transformation of the physical equations according to the adaptive mesh Reynolds transport 
theorem, \cf~Winkler, Norman, \& Mihalas~[13]. This changes the equations to an arbitrary 
frame of reference which can be quite different for each zone number. One discussion of 
how this transformation is carried out can be found in~[2]. The other way is the addition 
of the grid equation which has only numerical significance and which determines which frame
of reference for each zone at every time index is chosen. It needs input from physical
quantities in order to find the features one wishes to resolve.

We adopted the grid equation according to Dorfi \& Drury~[3]. Their basic idea is to 
balance the grid concentration
$$
       {\cal C}_k^n = {{\cal X}(r_{k+1}^n,r_k^n)\over(r_{k+1}^n-r_k^n)}       \eqno{(25)}
$$
with the resolution function
$$
       {\cal R}_k^n = \sqrt{1 + \Bigl[{\cal C}_k^n\Bigr]^2 \sum_j 
                               {\Bigl[{f_{j\;k+1}^n-f_{j\;k}^n\over 
                                 {\cal X}(f_{j\;k+1}^n,f_{j\;k}^n)}\Bigr]^2}} \eqno{(26)}
$$
of the problem under consideration.
The latter measures the arc length of various physical quantities $f_k^n$ which are 
typically the density, gas energy, radiative energy, or gas pressure.

Demanding the grid equation
$$
       {\cal C}_k^n \propto {\cal R}_k^n                               \eqno{(27)}
$$
under the constraint
$$
 {\alpha\over{\alpha+1}} \le {{\cal C}_{k+1}^n\over {\cal C}_k^n} 
                         \le {\alpha+1 \over\alpha}                    \eqno{(28)}
$$
with $\alpha > 1$ one forces the grid to distribute in such a way that steep features in
the fluid are resolved while maintaining approximately the same number of zones within
the unit arc length over the whole domain. The detailed derivation and formulation of the
grid equation can be found in~[3]. The function ${\cal X}$ renders the expressions in (14) \& 
(15) dimensionless and serves as a scaling factor. We define it in two different ways: linear
scaling with ${\cal X}(x_{k+1}^n,x_k^n) = x_{\scriptscriptstyle scale}$ and logarithmic 
scaling with ${\cal X}(x_{k+1}^n,x_k^n) = (x_{k+1}^n+x_k^n)$.

The grid equation itself contains further a time scale $\tau$ which sets the time step
below which the distribution of the zones remains frozen in.

In addition to the adaptive grid the user can choose a fixed grid in either the Lagrangean 
or the Eulerian frame of reference.

\subsection{The Artificial Stress}

An important computational tool in finite difference computations is the artificial 
stress. We apply a tensor formulation adapted from~[10] which is symmetric and 
trace free:
$$
 \PQ = - {4\over3} \rho \mq \left[{\p u\overp r} - {\mu\over2}{u\over r}\right]. \eqno{(29)}
$$
The artificial viscosity is evaluated according to:
$$
 \mq = C_q^1 \ell \cs 
     - C_q^2 \ell^2 \min\left[0,{\p(r^\mu u)\over r^\mu\p r}\right].             \eqno{(30)}
$$
It incorporates a linear term in the shock width $\ell$ proportional to the local sound 
speed and a quadratic term in $\ell$ proportional to the divergence of the velocity field. 
The latter is turned on only in compression regions. The shock width is given either by
a constant $\ell = \ell_0$ or set proportional to the radius $\ell = \ell_1\cdot r$. The 
artificial stress $\PQ$ appears in the momentum equation (3) in form of a shear force and 
in the energy equation (5) in form of stress work that dissipates some heat and functions 
as a source of entropy generation. 

In adaptive grid computation the artifical stress serves as a means to establish a well
defined ceiling for the discretization errors of the finite difference scheme. For this
reason the shock width $\ell$ can typically be taken considerably smaller than fixed grid 
algorithms demand. We often combine the constant scalings for the shock width, $\ell_0$, 
and the grid resolution, $R_{\scriptscriptstyle scale}$, and similarly choose the radius 
dependent expressions $\ell_1\cdot (r_k^n+r_{k+1}^n)$ and
${\cal X}(r_{k+1}^n,r_k^n) = (r_{k+1}^n+r_k^n)$ together. Because of this correspondence 
we can see that the shock width $\ell$ sets a scale for the concentration ${\cal C}$ that 
the adaptive grid must obtain to resolve the shock and hence defines the maximum 
resolution. This will be eluded further later on in the discussion of the test problems.

\subsection{Principles of Coding}

{\it Comprehensible}.
Each subroutine of has been given a header which contains a brief description
of the contents, a calling sequence for easier location within the code, and possibly
some references. In the body we make good use of comment cards to break the computation
into meaningful parts and to make references to corresponding sections in the code
documentation for comparison purposes.

{\it Versatile}.
One general guideline during the development phase was to keep functionality of the code 
operators separate from its arguments wherever is was meaningful. This applies in particular
to the linear solution precedure into which the Jacobian matrices of any set of equations 
can be fed. Hence they were formulated in a modular fashion so that they could be replaced
with any other set. In this regard our code becomes a general PDE solver for the 1-D 
case. Similarly we provide for a general interface for the constitutive relations so that
their specific evaluation is left to the user. The modular formulation of the physics
part also enables the user to adapt or extend the code to his/her particular problem
of interest.

{\it Portable}.
We have written {\TITAN} in standard FORTRAN 77. This allows for an immediate compilation and
execution on a host of workstations and mainframes.

\vfill
.