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 ni at each radial point ri. 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.