Our goal (for the purposes of this paper) is to construct self-consistent
models of expanding stellar atmospheres. The atmosphere itself is
parameterized by a number of parameters, e.g., the total energy emitted by the
object (luminosity *L*), the mass *M* of the star, the abundances of the
elements (in some cases as function of the location in the atmosphere). This
means that we have to find a set of physical variables such as temperatures,
densities, population number of each atomic energy level and the radiation
field, at *each* location in the atmosphere so that all constraint
equations are simultaneously fulfilled. In Fig. 1 we show this
requirement in a simplified graphical form where the arrows indicate *
direct* (usually non-linear) coupling between the different blocks. The number
of variables that need to be addressed is, formally, very large. A typical case
of a spherically symmetric shell model with 50 radial points requires a set of
50 temperatures and gas pressures (or matter densities). In addition to this we
include a set of about 6000 individual energy levels for atoms and ions
directly, the population of each must be known at every radial point, adding a
total of 300,000 variables. In order to calculate the population numbers, we
need a description of the radiation field at each radial point and on a set of
wavelength points (the rates that govern the transitions between atomic energy
levels are integrals of the mean intensity of the radiation field over
wavelength). We typically need 100,000 to 300,000 wavelength point to describe
the complete radiation field, which adds, in the worst case, 15 million
variables to the system.

Fortunately, most of these formal variables are tightly coupled to a much
smaller set of variables which we might, therefore, consider the
``fundamental'' variables of the model atmosphere problem. In our approach,
these fundamental variables are the temperatures *T*, the gas pressures
, and the population numbers *n*_{i} at each radial point *r*_{i}. The
radiation field is considered a ``derived'' quantity and the problem is thus
reduced to find a set of physical variables at each radial
point *i* so that the system outlined in Fig. 2 is self-consistent.
With this approach we have reduced the number of variables from several million
to a few 100 thousand, which is still a daunting number.

Although it is possible to analytically bring the system into a form so that it
could be solved by a Newton-Raphson approach [11], this idea is
computationally prohibitive because of its enormous memory and time
requirements (however, for smaller systems this approach has been used
successfully). Furthermore, this approach is complex to implement and it is
relatively hard to add more ``physics'' to the model atmosphere. We have thus
developed a scheme of nested iterative solutions that considers the direct (or
strong) couplings between important variables directly and iteratively accounts
for the indirect coupling between sets of variables. With this approach the
problem of constructing the model atmosphere can be separated into solving a
large number of smaller problems with only a few 100 variables. The global
requirement of a self-consistent solution is then reach by iteratively coupling
these sets of variables to each other until a prescribed accuracy has been
reached. This method works because the level of coupling between the variables
is very different. For example, the temperature structure of the atmosphere
depends mostly on the global constraint of energy conservation (represented by
wavelength integrals over the whole spectrum) and on the *ratios* of
several averaged opacities, but it does not depend strongly on the fine *
details* of the radiation field or the *individual* population of the vast
majority of the atomic levels. Therefore, *correction* to the temperature
structure can be calculated approximately. The current errors of, e.g., the
energy conservation equations, must be calculated exactly in order to this
scheme to function, however, this is relatively simple. The general idea of
calculating errors exactly but corrections to the variables approximately will
work if the approximations are good enough for the scheme to converge at all.
This method will require more iterations to reach convergence but this is more
than offset by faster individual iterations and (very often) by better
robustness. The latter is very important if many model atmospheres have to be
constructed or if no good initial guesses for the the variables are known.

In the following sections we will concentrate on a few key parts of the expanding atmosphere problem: the numerical solution of the radiative transfer equation for a single wavelength point (this will deliver the radiation field for a given set of variables at every wavelength), the solution of the NLTE statistical equilibrium equations (coupling the radiation field to the level populations), and an outline of the temperature correction procedure. The latter is important because it allows us to solve the NLTE statistical equilibrium equations separately for individual elements (and even ionization stages), which dramatically reduces the dimension of the sub-problems that have to be solved within the global nested iteration scheme. In this paper we will not discuss problems related to the hydrodynamics of the expanding medium or the details of the equation of state calculations, both of which are important topics.