MPS-ATLAS: A fast all-in-one code for synthesising stellar spectra

Context. Stellar spectral synthesis is essential for various applications, ranging from determining stellar parameters to comprehensive stellar variability calculations. New observational resources as well as advanced stellar atmosphere modelling, taking three dimensional (3D) effects from radiative magnetohydrodynamics calculations into account, require a more efficient radiative transfer. Aims. For accurate, fast and flexible calculations of opacity distribution functions (ODFs), stellar atmospheres and stellar spectra we developed an efficient code building on the well-established ATLAS9 code. The new code also paves the way for an easy and fast access to different elemental compositions in stellar calculations. Methods. For the generation of ODF tables we further developed the well-established DFSYNTHE code by implementing additional functionality, and a speed-up by employing a parallel computation scheme. In addition, the line lists used can be changed from Kurucz's recent lists. In particular, we implemented the VALD3 line list. Results. A new code, the Merged Parallelised Simplified ATLAS is presented. It combines the efficient generation of ODF, atmosphere modelling and spectral synthesis in local thermodynamic equilibrium, therefore being an all-in-one code. This all-in-one code provides more numerical functionality and is substantially faster compared to other available codes. The fully portable MPS-ATLAS code is validated against previous ATLAS9 calculations, the PHOENIX code calculations, and high quality observations.


Introduction
Electromagnetic radiation emitted from the stellar photosphere is one of the key sources of information about a star. To advance our understanding of stars, accurate spectroscopic and photometric measurements, as well as accurate modelling that allows us to connect the measured stellar electromagnetic radiation to properties of the stellar atmospheres and stellar interiors are essential.
The data gathered by the numerous ground-based and spaceborne telescopes that started observing during the last decade have highlighted the need for accurate modelling of stellar atmospheres and their spectra. For example, the Echelle Spectrograph for Rocky Exoplanet-and Stable Spectroscopic Observations (ESPRESSO, see Pepe et al. 2010) and the High Accuracy Radial Velocity Planet Searcher (HARPS, see Mayor et al. 2003) made high resolution spectral data for thousands of stars available, while the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST, see Ai et al. 2016) provides low-resolution spectra for millions of stars. The interpretation of these spectral measurements requires modelling of stellar atmospheres on a fine grid of stellar fundamental parameters, such as effective temperature, surface gravity, and chemical composition. e-mail: witzke@mps.mpg.de Furthermore, the advent of planetary hunting missions, e.g. Kepler (Borucki et al. 2010), TESS (Ricker et al. 2014), CHEOPS (Benz et al. 2020), and WASP (Pollacco et al. 2006) brought measurements of the photometric variability for several hundred thousands of stars. Even more data are expected from the forthcoming PLATO mission (Rauer et al. 2014). The main source for photometric variability of cool stars, like the Sun, are surface magnetic fields that affect the local structure in stellar atmospheres. Consequently, interpreting stellar photometric data requires an assessment of the effect of magnetic field in stellar atmospheres on the emergent radiation.
There are different approaches to modelling stellar spectral and photometric fluxes. One of the simplest, and most widely used approach relies on 1D modelling of stellar atmospheres under the assumption of radiative-convective equilibrium (with a simple parameterisation for convective flux and overshooting). While such a 1D approach has a number of shortcomings (see, e.g. Koesterke et al. 2008;Uitenbroek & Criscuoli 2011) it proved itself to be an invaluable tool for various applications (see, e.g. Castelli & Kurucz 1994;Claret 2000;Mészáros et al. 2012;Marfil et al. 2020) and is extensively used in stellar physics.
A more comprehensive approach relies on 3D hydrodynamic and magnetohydrodynamic (HD and MHD, respectively) simulations of near-surface convection in stars (see, e.g. Nordlund Article Stein 2012;Freytag et al. 2012;Magic et al. 2013;Beeck et al. 2015b,a). Using the simulated 3D cubes the emitted radiation can be calculated following a 1.5D approach, i.e. along many rays passing through such a 3D cube (see, Asplund et al. 2006;Riethmüller et al. 2014;Norris et al. 2017, for a detailed description of the 1.5D approach). An important advantage of the 3D MHD simulations and 1.5D approach over 1D radiative equilibrium modelling is that it allows us to directly account for the effects of the magnetic field on the emergent radiation and, consequently, to model the stellar spectral and photometric variability.
Altogether, there are two separate but similar challenges that demand comprehensive and fast spectral synthesis on an adjustable frequency grid, and broad band spectral intervals. First, the 1.5D modelling needs a huge amount of fast radiative transfer (RT) computations on 1D structures. Second, accurate and fast 1D atmosphere modelling with subsequent RT calculations for any stellar parameter is needed. The aim of the present study is to develop a fast and easily applicable RT code to address both these challenges.
The accurate treatment of the line opacity poses the main challenge to spectral synthesis over broad spectral ranges because of the immense number of spectral lines: where current line lists contain up to hundreds of millions of lines (Kurucz 2005b) that need to be taking into account. A careful treatment of both atomic and molecular lines is imperative because atomic and molecular lines are interspersed in particular in the spectra of cool stars. Spectral lines not only dominate some spectral regions (e.g. the UV), but also affect the atmospheric structure by blocking photons. While a forward spectral synthesis on a high-resolution wavelength grid, i.e. typically with a resolving power, R = 500000, is computationally expensive, not all applications require high resolution spectra. To that end, different techniques were developed to correctly account for line opacity, but reduce the computational cost on coarse resolution grids. The most commonly used methods are the opacity distribution functions (ODF) method and the Opacity Sampling (OS) method (Carbon 1984;Castelli 2005a). Both methods approximate the opacity using different ways of sampling it. We preferred the ODF method for developing the RT code, as it results in significantly faster RT calculations compared to the OS method, thus making it more suitable for the 1.5D approach. Recently, the ODF approach was further optimised to make it more efficient. Cernetic et al. (2019) found the best configurations of the ODF to reach several times faster computations while maintaining accuracy. Moreover, the ODF approach was extended to calculate stellar fluxes as they are observed in various filters (Cernetic et al. 2019;Anusha et al. submitted).
Widely used RT codes used for spectral synthesis of cool stars include the local thermodynamic equilibrium (LTE) codes MARCS (Gustafsson et al. 2008), and MAFAGS-OS (Grupp 2004), and the non-LTE codes PHOENIX (Husser et al. 2013), and TLUSTY (Hubeny & Lanz 2017). The PHOENIX code was developed to account for expanding atmospheres and deviations from LTE, and is therefore more complex and CPU intensive. Another successful RT code in LTE, which can calculate both model atmospheres in radiative equilibrium and the emergent spectra, is the ATLAS code by Kurucz (Kurucz 1970). Since for most applications mentioned above the spectral range does not include the extreme and far UV, it is sufficient to consider LTE codes. Thus, the ATLAS code includes all crucial physics for radiative transfer in the atmospheres of main-sequence stars while keeping the setup simple. Consequently, we chose to build on the ATLAS code. While the ATLAS code comes in two ver-sions; ATLAS12 (Castelli 2005a) and ATLAS9, we prefer the ATLAS9 version as it uses ODF whereas the ATLAS12 works with the OS method. The advantage of the ATLAS9 code is the short time required to compute a single model (a few minutes on a single core using the standard ODF table with 328 bins).
When updating ATLAS9 it is essential to also consider the DFSYNTHE code (Kurucz 2005a;Castelli 2005b). This code computes ODF tables for the ATLAS9 code. The main disadvantage of the ODF so far was the limitation to the chemical composition and microturbulent velocities for which the ODF tables were pre-tabulated. The generation of ODFs using the DFSYN-THE code was not suitable for massive computations as several routines had to be successively executed, and the computation time was very long (Kurucz 2005a;Castelli 2005b). To achieve our goal of fast spectral synthesis for arbitrary abundances, we need to eliminate this bottleneck, and in addition make the ODF generation more user friendly. This is achieved by merging the DFSYNTHE code, which calculates high-resolution opacities and uses them to obtain ODF, with the ATLAS9 code that can calculate both the atmospheric structures as well as the emergent spectra. Additionally, a more flexible treatment of the ODF was implemented. The frequency resolution, and ODF configuration can be changed, and a spectral filter included. Furthermore, we parallelise the code to speed up the high resolution opacity calculations and model calculations. Finally, the line lists used can be exchanged and we give example calculations with both, Kurucz's line list as well as the most up to date VALD3 line list (Ryabchikova et al. 2015). In this paper we describe the structure, and improvements of the resulting code, which we call the Merged Parallelised Simplified ATLAS code (MPS-ATLAS). Moreover, we validate MPS-ATLAS against previous ATLAS9 calculations, PHOENIX code calculations, and observations. The MPS-ATLAS code is available on request, and it will be made publicly available soon.
The paper is structured as follows. A brief summary of the physics and its implementation is given in Section 2 and Section 3, respectively. All code improvements are listed in Section 4, where we also discuss limitations of the code. In Sect. 5 we comment on the code performance and test the resulting emergent spectra by a code-to-code comparison and a code to observations comparison. The outcome of the work is summarised in Sect. 6, where also conclusions are drawn.

Radiative transfer calculations
The energy transport in stars defines the structure of their interior and atmosphere, as well as the emitted electromagnetic radiation. Focusing on the upper layers around the optical surface of a star, the radiative transport of energy is dominant. Thus, the imperative problem for modelling stellar atmospheres is solving the radiation transfer equation (RTE) along a ray for a time-independent system, where the subscript λ denotes the wavelength and implies that all quantities are monochromatic. I λ is the intensity, S λ is the source function, and µ = cos θ, where θ is the angle between the viewing direction and the normal to the stellar surface. The optical depth τ λ is determined from dτ λ = −χ λ (s) ds, where s is the height in the atmosphere. The total extinction coefficient per unit volume is χ λ = α λ + σ λ , which contains the absorption coefficient, α λ , and the scattering coefficient, σ λ . The source function is generally defined as the ratio of the total emission coefficient, j λ , to Solving the RTE becomes straightforward once the emission coefficient and opacity are known along the atmosphere. In local thermodynamic equilibrium (LTE) and under the assumption of coherent isotropic scattering, the emissivity can be expressed as where B λ is the Planck-function, and J λ is the mean intensity (Mihalas & Mihalas 1984). Then, calculating the opacity becomes the keystone to solving the RTE (except scattering coefficients).

Calculating opacity
The opacity κ λ (hereafter, we refer to the opacity normalised per unit of mass κ λ ≡ χ λ /ρ, where ρ is the density) can be decomposed into continuum opacity, κ λ,c , associated with different processes that involve atomic and molecular transitions with non-discrete wavelength (i.e. bound-free and free-free transitions), and line opacity, κ λ,l , due to discrete transitions, i.e. atomic and molecular spectral lines. To determine the total opacity κ λ = κ λ,c + κ λ,l , at each point, the atomic and molecular level populations have to be determined. Under the assumption of LTE, the atomic and molecular level populations in equilibrium can be calculated using the Saha-Boltzmann (SB) equation. This implies that they depend only on elemental composition (i.e. abundances), local temperature and pressure. In addition, in stellar atmospheres the Doppler broadening of lines occurs due to turbulence, which in 1D models is taken into account using the micro-turbulence parameter, v turb . Hence, for a given composition, opacity is a function of wavelength, local temperature, pressure, and micro-turbulence κ ≡ κ(λ, T, P, v turb ). In addition to the temperature and pressure values, calculating the level populations using the SB equation requires knowledge of the electron number density, n e , which however is not known a priori. Thus under the assumption of particle conservation the SB equation has to be solved iteratively. For a detailed description of the solver implemented in MPS-ATLAS see Appendix A and Kurucz (1970).
Having determined the populations, the continuous opacity, κ λ,c , and the line opacity, κ λ,l , can be calculated. For the continuous opacity, the MPS-ATLAS code takes the following contributors into account: Free-free (ff) and bound-free (bf) transitions in H − , H 2 , He, He − , C, N, O, Ne, Mg, Al, Si, Ca, Fe, the molecules CH, OH and NH, and their ions. While for the calculation of the equilibrium number densities and the line opacity also C 2 and CN are included, their contribution to the continuum opacity via photodissociation is neglected. Moreover, electron scattering and Rayleigh scattering on H I, He I, and H 2 were considered.
The line opacity is more expensive to compute, simply because of the huge number of lines for which the line absorption coefficient needs to be calculated. For atomic and molecular transitions from an initial energy level, denoted by i to a final energy level j, the line absorption is defined as: where ν 0 is the frequency corresponding to the transition, e it the elementary charge, m e is the electron rest mass, c is the speed of light, g i is the statistical weight of the level i, f ij is the dimensionless oscillator strength of the transition i − → j, and ∆ν D is the Doppler width. The gas density is ρ, the temperature is T, and E i is the energy of the initial energy level. The number density over partition functions for the entire ionisation stage, is given by N k /U k , where the subscript k indicates the ionisation stage. The term in brackets in Equation (4) (Cowley & Castelli 2002).
2.2. Atmosphere models in radiative equilibrium and emergent spectra Modelling a stellar atmosphere for a particular set of stellar fundamental parameters and elemental composition involves an iterative process of recalculating the atmospheric structure until radiative equilibrium (RE) is reached (see for example Collins 1989;Hubeny & Mihalas 2015, and references therein). In this process, in each iteration the RTE for a currently estimated atmospheric structure has to be solved in order to determine by how much the model structure needs to be corrected to satisfy the RE (for convergence criteria see Appendix B). The iterations are usually initialised with a starting structure that represents a converged solution for some other set of stellar fundamental parameters and composition. Before the iterative procedure starts, using the effective temperature of this starting structure, the temperature value at each of its depth points is re-scaled by applying the ratio of this effective temperature to the required one (for more details see Appendix C).
Then the RE calculations start. At each of the iterations, the wavelength-and depth-dependent source function (consisting of the thermal and scattering parts) needs to be determined. This can be achieved by various methods. By default MPS-ATLAS uses a Feautrier method, but a second iterative method is also implemented (more details see Appendix D). Having found the source function, the moments of the intensity are calculated, namely the mean intensity J λ , and the flux, H λ . This process is repeated for each wavelength point, typically in an interval from 9 nm to 10 000 nm, in order to obtain the total wavelength integrated Eddington flux, H.
The aim is to match the total Eddington flux, H, to the flux that corresponds to the desired effective temperature. Considering the flux, H, and keeping in mind that the majority of it comes from the region where 0.1 ≤ τ ≤ 2.0, the temperature corrections δT can be found in the optically thick regions. For that a modified version of the Avrett-Krook procedure is applied to take deviations from the RE due to convection and overshooting into account (for more details see Avrett & Krook (1963); Kurucz (1970)). In the optically thin regions, the temperature does not greatly affect the overall flux, whereas the derivative of it, dH ν /dτ ν = 0, is sensitive to a temperature change. The boundary condition for the temperature corrections δT in the optically thin part of the atmosphere requires the flux derivative of the Eddington flux to be zero (for a discussion see Mihalas 1978). This is achieved using a Λ correction method (Böhm-Vitense 1964). Then, the part of the atmosphere where both approaches overlap is smoothed out to match.
After the atmospheric structure is either calculated by the method described above or taken from some other source, e.g. from a ray through a 3D cube, and the source function is found, the emergent intensity can be calculated by evaluating the inte-A&A proofs: manuscript no. main gral The integral in Eq. (5) has to be computed for each wavelength separately to get the whole emergent spectrum. Note, that the framework described here is restricted to a plane parallel setup. Thus, when using a computed 1D model the emergent intensities for different view angles µ can be obtained by setting a set of µ angles in the code. To calculate spectra emerging from 3D cubes for different view angles a 3D cube is rotated and for each view angle a different set of 1D rays is obtained for which the RT is solved along the ray, but keeping the view angle µ = 1.

The ODF approach
The synthesis of the line opacity is computationally demanding due to the tremendous number of spectral lines to be taking into account. For example, the default MPS-ATLAS line list contains more than 100 million atomic and molecular lines (Kurucz 2005b). These millions of lines lead to very complex and rich spectra for cool stars. Consequently, a very fine wavelength grid is needed to catch all the details in the spectra. For a number of applications such detailed spectra are not required, and low resolution calculations are sufficient. However, lines still affect even low-resolution spectra. A straightforward way to include the effect of lines on low-resolution spectra is to simply average high-resolution spectra. Let us consider a small wavelength interval (hereafter, bin), e.g. between 0.1 and 10 nm wide, between λ i and λ i+1 . This bin represents one low resolution grid point, in which the intensity, I bin , representing the whole bin, has to be calculated The simplest way of approximating I bin is by obtaining I λ on a fine wavelength grid using a large number of points, N, and taking the sum where ∆λ j is the j−th discretised wavelength step. While this method is straightforward, the main issue is that I λ has to be calculated N times for an accurate approximation, where N for example is of order O(10 4 ) for a one nanometre interval in the UV, if the resolving power R is 500 000. An alternative method to avoid solving the RT many times is the ODF approach (see e.g. Hubeny & Mihalas 2015, pp. 625-627), which approximates the complex structure of the opacity. The main goal of the ODF method is to reduce the number of points N in a bin to a minimum, and still approximate the flux accurately. However, the opacity in one bin can abruptly change by several orders of magnitude within very narrow spectral intervals (see Fig. 1a). Consequently, a fine spectral grid is required to accurately describe the opacity profile. One could think about averaging the opacity in the entire bin to avoid the high-resolution calculations. In most cases such an averaging would lead to a gross overestimation of the line opacity effect, and would essentially lead to trapping photons that otherwise would escape (see Fig. 2 and its detailed discussion in Cernetic et al. 2019  by grouping points in a certain opacity range together. A straightforward way is to sort the points in one bin in an ascending order by their opacity value. The sorted opacity profile can be described using significantly fewer points than that of the original opacity (compare Figs. 1a and 1b). The opacity profile in each bin is usually divided into several sub-bins, whereafter the opacity is averaged over each sub-bin using the geometric mean (see Fig. 1 b). The geometric mean works better for the optically thick regime, because it avoids skewing the opacity towards large values. However, in the optically thin regime the arithmetic mean is better since the intensity depends on the opacity rather linearly. While both approaches are implemented, we used the geometric mean throughout this work in order to be consistent with the approach in the original version of the DFSYNTHE code.
Sub-bins might have different widths, which are defined by the fraction of the whole bin size. Generally, the part with the greatest opacity values is subdivided into smaller sub-bins, while the part with smaller opacity values is split into a few large sub-bins. The resulting step-function in the entire bin [λ i , λ i+1 ], where each step is the averaged opacity κ s,i of the s-th sub-bin, is called opacity distribution function. For a more detailed description of ODF and the importance of different sub-bin configura- Then, the intensity in each of the sub-bins can be calculated by solving the RT equation only once using the corresponding κ s,i value. For the entire i-th bin the intensity is obtained by summing over the contributions from the sub-bins, i.e.
where I s is calculated using κ s,i , and ∆λ s is the sub-bin width. Note that by re-arranging the opacity in a wavelength interval, as is done during the sorting in the ODF approach, two important assumptions are made implicitly: i) the opacity shape in the bin does not change rapidly along the line of sight, in particular, within the region of the atmosphere where the radiation for the corresponding bin is formed (Kurucz et al. 1974) ii) the wavelength interval of the bins are small enough that changes of the Planck function as well as changes of the continuum opacity can be neglected.
As mentioned in Section 2.1 the opacity for a given elemental composition, and turbulent velocity, v turb only depends on pressure and temperature in the LTE case. Thus, the ODF table contains the information of the mean opacity in each sub-bin, for a given elemental composition and micro-turbulence and can be pre-tabulated on a T-P grid (see e.g. Castelli 2005b, and references therein). For more details on the implementation see Appendix E. Having a pre-tabulated ODF table, it is now straightforward for any stellar atmosphere to obtain the total opacity for any temperature and pressure using linear interpolation on the T-P grid and adding the continuous opacity. Note, that when line opacity from the ODF tables is added to the continuous opacity the assumption that the continuous opacity in a particular bin is independent of wavelength is made implicitly. This is currently a limiting factor for the upper limit of the bin size. One way to circumvent this limitation is to include the continuum opacity in the ODF tables.

Code structure
The MPS-ATLAS code results in one monolithic executable which depending on a control input file can perform a different set of calculations. An overview of the original three routines on which MPS-ATLAS is based, is given in Fig. 2. Note, that the original ODF generation routine 1 consists of four separate codes, resulting in four executable. These three separate routines have been designed to perform the calculations outlined in Section 2 and correspond to the three MPS-ATLAS modules (see Fig. 3 and Fig. 4), where a module is a code internal execution mode. The well-established DFSYNTHE routine (Kurucz 2005a;Castelli 2005b) corresponds to module I in the MPS-ATLAS code, and generates ODF tables (see left column in Fig. 2). Two slightly different ATLAS9 codes 2 correspond to module II and module III in MPS-ATLAS, where one calculates 1D atmosphere models and the other the emergent intensity or flux (see middle and right columns in Fig. 2). The main difference between the two ATLAS9 codes is the main routine. The code for calculating model atmospheres contains routines to perform the temperature corrections in addition to solving the RT.
The original collection of codes is intricate to run due to several codes that need consecutive execution. Furthermore, the three different codes make use of similar procedures. For example, both the DFSYNTHE and the ATLAS9 codes require calculations of the populations and continuous opacity. While these calculations can be obtained using the same functions, both codes have their own copy of these functions with slightly different implementations. This renders any modification quite difficult. Therefore, we chose to merge three codes into one code with different execution modes, that we call here modules, that take care of the different calculations. The resulting, merged code is further adjusted, parallelised, and simplified, and thus called Merged-Parallelised-Simplfied-ATLAS (MPS-ATLAS). MPS-ATLAS is mainly written in Fortran 90 (free form), but some parts are left in F77. We will rework them in future releases. We use dynamical memory allocation in several modules, for example in the ODF calculations.
The advantages of the MPS-ATLAS code are a more user friendly operation, significant computational speed-up, and wider functionality (in particular flexible ODF setup, which is discussed in Sec. 4). A schematic diagram in Fig. 3 shows the overall structure of the MPS-ATLAS code. The new code has three execution modes that correspond to the original structure of the three codes: ODF generation, model calculations, and emergent intensity or flux calculation. These three modules can be either executed consecutively or separately, which is controlled during run-time. Both model calculations and emergent intensity calculations need ODF tables as input, which is indicated by the arrow. As discussed in Sect. 2.2, the model calculation is not necessary if alternative input atmospheres are used instead, for example a ray from a 3D radiative MHD cube (such a bypass of module II is indicated by the arrow in the bottom of Fig. 3).
In the following we give an overview of the code structure in terms of input and output. There are two types of parameters (shown in Figure 4): i) overall input parameters (highlighted in green), that are specified once and do not change for a particular calculation, for example the elemental composition ii) parameters that are set in the input, but which determine a grid on which the output calculations are performed (highlighted in grey).   velocity. In addition, the temperature and pressure grid, as well as frequency range and the bin, sub-bin configuration need to be specified (more details in Appendix E). An example input file for the ODF calculation is shown in Appendix G. The output ODF is needed as input for both module II, and module III. Module II calculates a 1D atmosphere model in RE. As input parameters this module needs the elemental composition, effective temperature (T eff ), surface gravity (log g), and the mixinglength parameter, which is needed to include convection through mixing length theory (Böhm-Vitense 1958), for which overshoot can be turned on. In addition, an ODF table is required as input to account for the effect of line blanketing on the atmosphere structure. For the opacity table either the output from module I, or from a pre-calculated ODF table gird can be used. Moreover, an initial 1D model atmosphere is read, which is used as an initial starting point, and has the same format as the output model (for more details see Appendix C). The output model is obtained on an a priori specified τ Ross −grid. An atmosphere model consist of column mass, temperature, total gas pressure, electron number density (n e ), Rosseland mean opacity (κ Ross ), radiation pressure, and micro-turbulence velocity (v turb ) for each depth point. Currently, we only consider cases of height independent v turb , though for height-dependent micro-turbulence it is possible to use several ODF tables covering a range of v turb values and interpolate between them.
The third module either synthesises the emergent flux or the emergent intensity for different view angles. Here, the only actual free parameters to set are the wavelength range and viewing angles, as the input is already pre-determined by the ODF and model atmosphere. While the elemental composition should match the ODF, the wavelength grid on which the ODF is calculated sets the wavelength grid of the emergent intensities. Moreover, the ATLAS output 1D atmosphere model in RE contains all the necessary information to calculate the emergent spectrum, in particular the electron number density, and the radiative pressure at each depth point obtained with the elemental composition for which the model was calculated. On the contrary, if using an external model for example from a ray through a 3D MHD cube, several quantities might be unknown. In particular, the electron number densities, which are needed to find the populations of all other ions, have to be obtained. In order to use external models that only have the column mass, pressure, and temperature, an additional flag was introduced to recalculate the equilibrium number densities in each depth point for the given elemental composition. Example input files are discussed in more detail in Appendix G.

Functionality extensions
In the following we describe the main functionality extensions included in the MPS-ATLAS code.

Flexible wavelength grid and optimised sub-binning
The DFSYNTHE code is limited to two particular wavelengthbin grids onto which ODFs are calculated, both with the same wavelength independent sub-bin sizes. However, various applications require spectral synthesis with different spectral resolutions. Hence, it is important to have an option for changing the ODF wavelength grid. Furthermore, Cernetic et al. (2019) showed that an optimal choice of the sub-bin grid leads to significant improvements in both accuracy and speed of the calculations. In particular, they showed that for most of the wavelength range two to four sub-bins are sufficient to reach the same accuracy as the standard 12 sub-bins (which are hard coded into DFSYNTHE). Such a threefold speed up is especially important for 1.5D calculations which require an immense number of 1D calculations (about a million per cube with 1000 times 1000 horizontal grid points).
In this context, we implemented a flexible treatment of the bin grid, enabling the user to choose the overall wavelength interval on which the opacity is calculated (keeping the maximum range between 9 nm and 10 000 nm), as well as the binning and sub-binning on this interval (more details on the implementation can be found in Appendix G). The underlying high-resolution opacity is then only calculated in the chosen wavelength interval to save computational time. The flexible wavelength binning is in particularly useful if synthesised spectra need to be compared to observations of different resolving power. In order to reduce the computational cost, the number of sub-bins per bin can be reduced significantly if an optimal sub-bin border distribution is chosen. In addition, having a flexible binning and sub-binning grid allows rapid computation of ODF tables for broadband filters if only the flux through a particular filter is needed (Cernetic et al. 2019;Anusha et al. submitted).
Moreover, we have added an option to pre-tabulate the opacity on a high-resolution grid. For that, the code reads the starting and end wavelength of an interval, and uses a resolving power of R=500 000 in this interval. The opacity is then written in the same format as an ODF, but using only one sub-bin, which contains the opacity of the wavelength point. This can be used to calculate emergent intensity at high spectral resolution. Since the line opacity on a high resolution grid might have strong changes with temperature and pressure, the interpolation error is potentially greater compared to an ODF. An exemplary calculations of the Vega flux in the wavelength region of the Balmer jump shows that the error between the flux using the interpolated opacity (using a T-P grid that splits the same ranges as used in Castelli (2005b) into 100 logarithmically equidistant steps) and the calculated high-resolution opacity can reach up to 3% (see Fig. 5).
The largest deviations occur in line wings of strong lines (see Fig. 6), which indicates that the temperature steps in the considered T-P grid are insufficient to account for the line broadening sensitivity. Therefore, the T-P grid can be adjusted depending on the needed accuracy (for more details see Subsec. 4.2).

Flexible T-P grid for pre-tabulation
The DFSYNTHE code uses a predefined T-P grid of 57 temperature and 25 pressure values for pre-tabulating ODF tables. The pressure covers 12 orders of magnitude ranging from 10 −4 to 10 8 , and the temperature ranges from 2 × 10 3 − 2 × 10 6 K.
In MPS-ATLAS the number of T-P points and their values can be specified (as input parameters). This has two advantages. On the one hand, the numerical cost can be reduced, while the resolution of the T-P grid is kept, but the ranges of temperatures and pressures are adjusted to certain types of stars. On the other hand, the resolution can be increased if more accurate interpolation is needed. This is especially important if the code is used to generate pre-tabulated high resolution opacity (see Fig. 5).

Line pre-selection criterion
The most extensive line lists such as Kurucz's original line list 3 and VALD 3 (Ryabchikova et al. 2015) contain many lines that do not contribute to the opacity within the temperature and pressure range of the ODF table. To reduce the computational cost in the line opacity calculations, the lines are pre-selected based on their line core opacity and the continuous opacity for each T-P point. For that, lines which lead to an opacity in the line core of several orders of magnitude (controlled by the cutoff factor) less than the continuous opacity are discarded, see Appendix E for the selection criteria. By default the cutoff factor is set to 10 −3 . Here, we tested how the solar surface intensity is affected by choosing several smaller cutoff factors (see Fig. 7). Since a smaller cutoff factor increases the computational time significantly, the smallest value we tested is 10 −12 . The difference between a cutoff factor of 10 −9 and 10 −12 is negligible. However, overall the largest decrease in intensity can be observed around 200 nm, where it reaches up to 5%. This implies that even very weak lines have a non-negligible effect on the flux, in particular if there are many of them.

Changing line lists
Kurucz's original line list includes an immense number of theoretically computed lines which have not been measured experimentally. Over the last few decades the atomic and molecular data has been constantly updated to obtain more accurate data. Recent developments led to various line lists which are updated regularly. In particular there have been several important updates of the Vienna Atomic Line Database (VALD), which is a database focusing on atomic data relevant for modelling of stellar atmospheres (Ryabchikova et al. 2015). While the DFSYNTHE code only works with Kurucz's line list, MPS-ATLAS can be used with different line lists. For the current tests we exchanged the atomic line list, but we kept the line list for diatomic and H2O lines (Schwenke 1998) Fig. 7. Ratios of the solar emergent intensity at the disc centre calculated using ODF tables that were generated with different cutoff factors for the line pre-selection. Here cutoff factor refers to the opacity of a given line relative to the continuum opacity. If this ratio drops below the cutoff factor, this particular line is not included in the computation of the intensity spectrum. The default value of the cutoff factor is 10 −3 .
To illustrate the effect on the flux caused by the differences in the line lists, we first calculated the model atmosphere for the Sun using the standard 'little' grid ODF (Castelli 2005b) with the original line list and the VALD3 line list. Then, using these models, we calculated the intensity at disc-centre on a high-resolution wavelength grid with resolving power R=500 000 for each of the two line lists. Figure 8 shows both intensities in the range between 5340 Å-5360 Å together with the Hamburg atlases of the solar spectrum 4 (Neckel 1999;Doerr et al. 2016). For this comparison, we degraded the resolving power of the calculated flux using a Gaussian kernel. It becomes evident that the intensities obtained using the VALD3 line list agree significantly better with the solar observations, while the intensities obtained using Kurucz's original line list show quite a few lines that are not present in the solar spectrum. This is because Kurucz's line list contains significantly more lines than VALD3 and most of these excess lines have never been measured in the laboratory.
In a second step, we compare the overall energy distribution for the Sun, and a K-type star (T e ff = 4000K) obtained using different line lists. For that, we calculated the model atmospheres using the 'little' grid ODF with different line lists, and subsequently generated the overall disc integrated emergent flux using the standard 'little' grid ODF with the corresponding line list. Figure 9 shows that the fluxes calculated with the original Kurucz's line list are significantly different from the fluxes obtained using the VALD3 line lists below 450 nm. The effect in the UV is huge for the Sun, where the flux obtained using the original Kurucz's line list gives a better agreement (see a more detailed discussion in Sect. 5.2). However, for a K type star for which line opacity is even more important than for the Sun there is rather a moderate difference. This indicates that the main dif-ference between the VALD3 and the Kurucz's line lists are lines whose lower level have higher energies (i.e. mainly lines contributing to opacity at higher temperatures than in K-stars).

NH photodissociation
Comprehensive comparison of observed and modelled solar spectra showed that molecular photodissociation is an important source of continuum opacity in the UV range (Fontenla et al. 2011). Moreover, Fontenla et al. (2015) proposed that NH photodissociation is the source of the missing opacity in the UV (see, e.g., Busá et al. 2001;Short & Hauschildt 2009;Shapiro et al. 2010, for the detailed discussion of the missing opacity problem). So far only the molecules CH and OH were included in the ATLAS9 continuous opacity calculations (Castelli 2005a). We added the opacity contribution from the NH molecule using cross-sections that we obtained from Kurucz (2020). A more detailed description on how molecules are included in MPS-ATLAS is given in Appendix A, and how CH, OH, and NH photo-dissociation opacities are implemented in Appendix F. This additional opacity source changes the specific emergent intensity by at most 0.5% (see Fig. F.1). We conclude that NH is not sufficient to explain the missing opacity in the UV.

MPI parallelisation for faster computations
While the RT calculations can be sped up by using optimised ODF configurations, the calculations of the ODFs table itself remain time-consuming. To achieve faster calculations we parallelised the module for ODF generation (module I in Fig. 3) along temperature values. We note that the number of lines contributing noticeably to opacity depends on temperature. In particular, more lines have to be calculated for lower temperature values so that producing ODFs for them is the most time consuming. To account for this imbalance of computational time, we used a master-slave implementation. Here, one core distributes the T values for which the ODF must be calculated to all other cores as soon as they finish with their previous values. Consequently, while different cores calculate ODFs for different numbers of temperature values, the number of calculations is distributed between the cores more or less equally. This also has the advantage that the number of T values does not need to be divisible by the number of cores used. With this implementation, an example ODF table on the same T-P grid as used in the original DFSYN-THE code can be calculated in under 10 minutes on 36 cores on a modern high performance cluster.
For the RE calculations (module II in Fig. 3), the integrated flux over the whole spectrum is needed, where the number of frequency points can be varied. Calculations with a large number of frequency points, or if the high-resolution grid is used, are very time consuming as the RT needs to be solved in each iteration on chosen frequency grid. Thus, we also parallelised the RT calculations along wavelengths in the module for modelling atmosphere structures.
Finally, module III that calculates the emergent intensity (see module III in Fig. 3) uses a parallelisation if the intensity for more than one model has to be calculated. This is only the case for post-processing 3D calculations. In this case the atmosphere models are distributed among the available cores.

Speed-up & portability
Using the parallelised version of the ODF module we tested the MPS-ATLAS scalability. For the ODF calculation we used the standard T-P grid as in Castelli (2005b), but set the wavelength interval from 200 nm -10 000 nm, which significantly reduces the computational time. The low temperature values in the standard grid (below 5000 K) require the highest number of lines, and thus are very-time consuming compared to larger temperature values. Fig. 10 shows the speed-up for ODF generation and RE calculations. It becomes evident that there is no further speed-up of the ODF calculations for more than 10 cores. This is because the lowest temperature value (1995 K) forms the bottleneck in this example. Even if all other temperature values are calculated faster using more cores, at the end of the entire ODF calculation the temperature value that takes the longest determines the amount of time needed. The number of cores at which such a bottleneck (due to one particular temperature) appears depends on the T-P grid, and wavelength range of the ODF table. In general, the fewer P points and more T points a grid has, the higher the number of cores until the speed-up saturates.
To test the parallelisation for the RE calculations, we chose to model a 4000 K star with logg = 4.8, and M/H = 0.3. We used a standard 'little' grid ODF (Castelli 2005b), and reduced the wavelength interval to the first 1210 points on the grid, which corresponds to 9 nm -10 000 nm. The model calculation converges after 37 iterations. The speed-up of the RE calculations (module II) shows a different behaviour: the increase in speedup starts to flatten after 8 cores (see Fig. 10 b). The flattening happens because only the part of the RT calculations along the wavelength points is parallelised, while all other calculations, especially the temperature corrections are executed on one core. Here, the behaviour of the speed-up depends on the number of wavelength points that need to be calculated. The parallelisation for the emergent intensity calculations (module III) is along atmosphere models, thus it is only useful for calculating along model atmospheres from 3D cubes. Since it only needs a few communications in the beginning to distribute the input settings, we expect an almost ideal speed up. The code is fully portable. The available git version of the code includes an installation script that downloads all necessary libraries and files, such as a NetCDF library, line lists, and example input files. The code can be compiled using an openMPI and intel compiler, as well as a corresponding gnu compiler. Note that for the gnu version, a different line list format is needed, which we provide on request.

Benchmark
In order to put the MPS-ATLAS code into context, we compared emergent flux obtained by the MPS-ATLAS code to that from other codes and to observations. For the code-to-code comparison the stellar parameters were chosen to be available on both the PHOENIX-and Kurucz-model grid of calculated specific intensities. To this end, we considered three hypothetical stars of spectral classes A, F, and K (with corresponding temperatures of 8000 K, 6500 K, and 4000 K). These spectral classes have been chosen to test the performance of the code for cases when opacity is dominated by different sources. For example, in a Ktype star the opacity is mainly due to by millions of atomic and molecular lines while the main opacity source in an A-type star is the continuum and hydrogen lines. Furthermore, for comparison to observations, we used the Sun and Vega (see Table 1 for the summary of the fundamental stellar parameters of the stars we used). We chose Vega for comparison, because accurate measurements of the total flux over a large wavelength range are available.

Code-to-code comparison
We compare emergent fluxes calculated using the MPS-ATLAS code to the original fluxes from the ATLAS9 grid calculated by Kurucz 5 (Kurucz 1993), and to the modelled fluxes on the PHOENIX grid 6 for three different stellar types (K-type, Ftype and A-type; see Tabel 1 for details). In our calculations we assume the same setting as in the Kurucz grid modelled fluxes: convection is turned on, without overshoot, and the mixing length is set to 1.25. We further consider a constant microturbulence of 2 km/s. Generally, the PHOENIX code provides a more intricate setup, e.g. the models are calculated in spherical geometry, condensation is included in the equation of state, and some NLTE effects can be turned on. However, the model grid we use (Husser et al. 2013) was obtained using LTE and the RT calculations only make use of the NLTE effects by a special line profile treatment for some species (Li I, Na I, K I, Ca I, Ca II). Moreover, the mixing length and the micro-turbulence parame-5 http://kurucz.harvard.edu/grids.html 6 http://phoenix.astro.physik.uni-goettingen.de ters vary on the PHOENIX model grid, and thus can slightly differ from the one used in ATLAS9 and MPS-ATLAS (for more details see Husser et al. 2013). We preformed calculations using two different elemental compositions. For a comparable calculation to the Kurucz grid we used the same element abundances as in Kurucz's calculations. These abundances are taken from Anders & Grevesse (1989), and will be referred to as 'Anders composition'. For the comparison with PHOENIX we used the more up-to-date 'Asplund composition' taken from Asplund et al. (2009). For the model-to-model comparisons, both the model atmosphere and the emergent spectra are calculated using the standard 'little' ODF with 1221 frequency points (Castelli 2005b). We smoothed all spectral fluxes by applying the average over a 15 nm interval around each wavelength grid point. The PHOENIX highresolution spectral fluxes are first averaged over the same intervals as given by the ATLAS wavelength grid, and then smoothed with the same procedure over 15 nm intervals.
For an A-type star the emergent spectral flux returned by the MPS-ATLAS code and those from PHOENIX and ATLAS9 grids are shown in Fig. 11a. For a more detailed comparison we show the deviations of the calculated MPS-ATLAS fluxes, and the PHOENIX to the original ATLAS9 fluxes in Fig 11b. Overall there is a reasonable agreement between the three codes. The largest deviations are in the UV for wavelength shorter than 400 nm. A good example for the importance of the elemental composition can be seen in the range between 210 nm and 290 nm. differences can occur especially in the UV due to different line lists used. In the wavelength interval between 300 nm -1000 nm the PHOENIX calculation gives an overall smaller spectral flux than fluxes calculated with MPS-ATLAS and ATLAS9. The difference is approximately 5% around 400 nm and then decreases to a few percent for wavelengths greater than 1000 nm. In this wavelength interval the dominant continuum opacity contributor are H − bound-free transitions. Thus either a slight difference in the equilibrium number densities or a different implementation of the H − bound-free cross-section can potentially lead to these differences.
In addition there are three larger deviations around 820 nm, 1458nm, and 2279nm which correspond to the Paschen limit, the Brackett limit and the Pfund limit, respectively. This implies that all three codes have a different treatment of the hydrogen continuous transitions of these series. On the contrary, the Balmer series transitions, H-α (at 656nm) to H-ζ (at 389nm), show only very slight differences between the three codes. The significant deviations around 395nm are caused by the Ca II H and K doublet, where the MPS-ATLAS calculation gives more similar results to PHOENIX than to the older ATLAS9 calculations. We note, however, that the proper treatment of the Ca II H and K doublet requires NLTE modelling which is absent in all three codes.
In Fig. 12 the emergent flux for a F-type star (see Table 1 for exact fundamental parameters) is displayed. Overall the deviations in the UV are greater and reach up to 20%, while the differences in the visible and infrared are smaller. The two elemental compositions entering the MPS-ATLAS lead to the following differences in the spectra in the range 200 nm -500 nm: the flux is larger for the 'Asplund composition', and it decreases towards 450 nm whereafter the flux is smaller for the 'Asplund composition' than for the 'Anders composition'. This behaviour is most probably caused, directly or indirectly, by the effect of the composition on line blanketing. The line opacity, in particular from iron, is smaller for the 'Asplund composition' and allows more photons to escape. Then, the decrease in the visible can be explained by a compensation effect as the wavelength integrated flux has to be the same for both element compositions. The PHOENIX flux oscillates around the MPS-ATLAS flux obtained using the 'Anders composition' in the region 200 nm -400 nm, which implies that the line opacity differs significantly due to a different line list. Finally, the wave-like deviation between 500 nm -2600 nm indicates that the photospheric structures differ somewhat from each other. It becomes evident that the PHOENIX, and MPS-ATLAS model have a similar structure, while the older Kurucz's model deviates.
The deviations between the codes are larger for a K-type star (see Fig. 13). In the UV the difference between PHOENIX and other codes reaches 160%, this might be due to a smaller number of UV lines contributing to opacity at lower temperatures (i.e. below 4000 K) in the line list used in the PHOENIX code. The flux calculated using MPS-ATLAS deviates by up to 25% from Kurucz's originally computed flux. Focusing on wavelengths longward of 400 nm, the deviations oscillate, but the amplitude of the deviations decreases towards the infrared. Comparing the difference between the 'Anders composition' and the 'Asplund composition', the flux obtained using the 'Asplund composition' is closer to the PHOENIX calculation for most of the wavelength in the visible. This is also the case in the infrared, where the overall deviations become less than 10% between all four models.  Woods et al. 2009, for a detailed description). Solar irradiance (top panels) calculated using MPS-ATLAS with three different assumptions: no convection, with convection but no overshoot, and with convection and overshoot. The flux deviations between WHI observed irradiance and models in % are shown in the bottom panels.

Code-to-observation comparison
Besides comparing with the output of other codes, it is important to also check how well MPS-ATLAS reproduces observations. For the Sun a lot of accurate high-resolution data exist, as well as low resolution spectra. In contrast, for most of the stars either broadband fluxes, intermediate resolution spectra or normalised (e.g. to continuum level) high-resolution spectra for a rather small wavelength interval are available. Here we first test how well the solar flux can be matched by the MPS-ATLAS code, and in a second step we chose Vega as a comparison star.

The Sun
For the calculations presented in this subsection both the model atmosphere and the emergent spectra are calculated using the standard 'little' ODF with 1221 frequency bins (Castelli 2005b). For all considered cases we recalculated the model when changing any parameters. We always assumed a micro-turbulence of 1.5 km/s and if convection was turned on, the mixing-length is 1.25. We compare our calculations to Solar Irradiance Reference Spectra (SIRS) for the 2008 Whole Heliosphere Interval (WHI Woods et al. 2009). For a better comparison all calculated fluxes are smoothed out using a trapezoidal kernel (Harder, private comm.) to match the WHI resolution in wavelength regions where the 'little' ODF resolution is higher than that of the WHI. We kept the original resolution otherwise. Fig. 14 shows the solar irradiance at a given wavelength at one AU from the Sun as in the WHI compared to MPS-ATLAS calculations using the 'Anders composition' with and without convection, and overshoot. One can see that in the mid UV the deviations between the calculated and observed irradiance significantly exceed the uncertainty of the measurements. This is not surprising (and similar deviations are also present in the AT-LAS9 and PHOENIX calculations, see below) since it is well known that currently available line lists miss a significant number of weak lines (this excess flux constitutes the famous 'missing UV opacity' problem, see e.g. the recent review by Rutten 2019, and citations therein). We note that the solar UV spectrum is also affected by deviations from LTE, which are ignored in MPS-ATLAS calculations. However, it has been shown that ac-counting for non-LTE effects will make the deviations between the calculated and observed UV flux level even larger (see, e.g., Short & Hauschildt 2009;Shapiro et al. 2010;Tagirov et al. 2019).
In contrast to the UV spectral domain, the modelled irradiance agrees very well with the observations in the visible and infrared. Fig. 14 shows that the deviations there are mostly within the minimum measurement uncertainty. While the differences between the three calculations are very small, the agreement between the observations and the MPS-ATLAS calculation including convection with overshoot is best in the interval 450 nm -650 nm. However, in the region 650 nm -900 nm the fluxes including overshoot show the greatest deviations. The wave-like behaviour of the deviations in the range between 500 nm -1500 nm might indicate that the modelled atmosphere temperature in the region where the continuum is formed does not accurately match the actual temperature in the Sun.
The solar elemental composition has been intensively investigated for the last decades and updated following more accurate modelling (Anders & Grevesse 1989;Grevesse & Sauval 1998;Asplund et al. 2009). Changing ratios between hydrogen, helium and heavier elements, which significantly contribute to the electron number density, affect the structure of the atmosphere and the spectral synthesis in a competing way. Increasing the concentration of the electron donors without re-calculating the atmospheric structure leads to a drop of the flux, and thus a decrease of the effective temperature. However, as soon as the structure is re-calculated in RE, it compensates for the decreased flux and the effective temperature should return to its original value.
In the next step we show the difference in the irradiance for the 'Anders composition' and 'Asplund composition', both with convection but no overshoot. Fig. 15 shows flux model calculations for the different elemental compositions, together with the WHI measurements. The difference in the elemental composition leads to a redistribution of the flux in different wavelength intervals due to a change of equilibrium number densities of the species (which in turn affects the opacity). The effective temperatures obtained from the total flux are almost identical for the two compositions, being T eff = 5778.9K for the 'Asplund composition', and T eff = 5778.5 for the 'Anders composition'. Nevertheless, the irradiance calculated using the 'Asplund com-   position' has a higher flux in the UV, compared to both the WHI and the flux using the 'Anders composition', so that the former displays greater deviations from the observations. In contrast, the flux calculated using the 'Asplund composition' is lower in the visible and IR compared to that calculated using the 'Anders composition'. Note, that the same behaviour was also observed for the A-type and F-type star of solar metallicity (see discussion in Subsection 5.1 and Figs. 11-12). Both calculations match the observations within the measurement uncertainties for wavelengths greater than 450 nm.
Finally, in Fig. 16 we plot the flux obtained by MPS-ATLAS for the 'Anders composition' without overshoot, together with the original fluxes calculated by Kurucz 7 , and the PHOENIX flux calculations. Since the PHOENIX grid does not have the flux for the solar effective temperature, and surface gravity, we used linear interpolation between the closest available stellar parameters. Furthermore, before applying the same smoothing procedure using a trapezoidal kernel to match the WHI resolution, 7 http://kurucz.harvard.edu/stars.html we averaged the PHOENIX flux over the same wavelength intervals as given in the standard 'little' ODF grid. The interpolated PHOENIX flux shows a slightly worse agreement with the observations in the UV, but an overall good agreement in the visible and IR. Kurucz's original calculation agrees better with the observations in the interval 450 nm -650 nm than the PHOENIX flux. All three codes show a wave-like behaviour between 500 nm -1500 nm indicating a slight mismatch of the modelled atmosphere as discussed above.
Note that we used Kurucz's original line list for the above atmosphere model and flux calculations. The modelled flux in the UV is significantly larger than the observed values, and using the VALD3 line list leads to an even greater flux in the UV (see Fig. 9). This indicates that Kurucz's line list leads to a better agreement of the calculated UV flux to observations. This is because Kurucz's line list contains significantly more lines than VALD3 and even though most of these lines have never been measured in the laboratory they allow a better representation of the UV line haze than the VALD3 lines.

Vega
The other star for which accurate absolute spectrophotometry is available is Vega. The most up-to-date observed absolute flux consists of data from two instruments, the International Ultraviolet Explorer (IUE), and the Space Telescope Imaging Spectrograph (STIS) mounted on the Hubble Space Telescope (Bohlin & Gilliland 2004;Bohlin et al. 2014). The fundamental parameters for Vega still contain a non negligible uncertainty (Castelli & Kurucz 1994). Therefore, we tested a small range of stellar parameters to seek good agreement, where the range was chosen based on previous stellar parameter determinations (Castelli & Kurucz 1994;Catanzaro 2014). Namely, for the model comparison we calculated a small set of models with slightly different surface gravity values (logg = [3.90, 3.95, 4.00]), and different metallicity (M/H = [−0.3, −0.5]), but the same effective temperature, T eff = 9550K. The 'Anders composition' was used for the abundances.
For the search of the best stellar parameters, we used the flux calculations obtained using the standard 'little' ODF with 1221 frequency bins (Castelli 2005b). The spectral flux observed for Vega was averaged on the same wavelength grid, and subsequently we broadened the fluxes using a Gaussian kernel with a full width at half maximum (FWHM) of 3nm below 350 nm, and with a FWHM of 5nm above 350 nm. The output of MPS-ATLAS is the flux at top of the stellar atmosphere. To connect it to the observed flux one needs a scaling factor, sf = (d Vega /R Vega ) 2 , where R Vega is Vega's radius and d Vega is the distance of Vega to the observer. This factor is still uncertain and some assumptions must be made for the compar- ison (see Castelli & Kurucz 1994, where is was estimated to be sf = (1.62 ± 0.07) × 10 16 ).
Here, we take this factor into account by taking the ratio of the observed flux to the modelled flux in the wavelength region 400 -600 nm, where the ratio showed the weakest dependence on the metallicity and surface gravity values. The scaling factor we obtained is sf = 1.59 · 10 16 , which is within the uncertainty of the estimation. Overall, the best agreement of the modelled flux to the measured flux was achieved using M/H = −0.5 and logg = 3.90, but very little difference is found between the flux for logg = 3.90 and the flux for larger surface gravity of logg = 4.00.
For the comparison to PHOENIX, we downloaded the PHOENIX emergent intensities for T eff = 9400K, and T eff = 9600K, and the surface gravity values logg = 3.50 and logg = 4.00, with the metallicity M/H = −0.5. After calculating the fluxes, we had to interpolate to get to the effective temperature of T eff = 9550K, and a surface gravity of logg = 3.90. Due to the limited wavelength range for which the PHOENIX emergent intensities are provided we can not calculate the exact effective temperature from the flux.
The PHOENIX flux is averaged in exactly the same way as the IUE and STIS data, and subsequently all fluxes are broadened using a Gaussian kernel. Figure 17 shows the comparison in the UV (120 nm -350 nm) between the measured flux, the flux calculated using the MPS-ATLAS code and the PHOENIX flux. The MPS-ATLAS calculation has an overall better agreement with the measurements than the PHOENIX modelled flux, but still deviates significantly for wavelength shorter than 260 nm. Note that the elemental composition in the PHOENIX flux is the 'Asplund-composition' by default, but scaled to M/H = −0.5.
The comparison in the range 350 nm -1000 nm is displayed in Fig. 18. The MPS-ATLAS calculations show very good agreement, with only a few spectral intervals with greater deviation than the measurement uncertainty. In contrast, there is an offset in the PHOENIX flux, whose deviation from the observations shows a significant wavelength dependence that continues into the UV range (see Fig. 17). While an overall shift could be explained by the scaling factor (which was taken to be the same as for MPS-ATLAS), the wavelength dependence cannot be removed by a different normalisation. This indicates that the atmospheric structure or the continuum opacity has a slightly worse agreement. Moreover, the Paschen limit (820.4 nm) deviates significantly. The agreement of the modelled flux using MPS-ATLAS with the measurement is very good in the visible and near-infrared.
Finally, we compared calculated and observed highresolution spectra around hydrogen Balmer and Paschen lines. For that, we used the same atmospheric structure as before, but we calculated the high-resolution flux with resolving power R = 500000. Subsequently, we used Gaussian broadening to degrade the resolving power to R = 400. In Fig. 19 the observed and calculated flux together with the deviation of the calculation to the observations are shown. The H α line shows the best agreement with at most 4% deviations, while the deviations for the Paschen lines are up to 6% and for higher Balmer lines even up to about 10%. We note that in a previous comparison (Bohlin & Gilliland 2004) a better agreement (with deviations of a few percent) was achieved. However, these differences are mainly attributed to differences in the line depth of the updated STIS data (see Figure 4 in Bohlin & Gilliland (2004)), while some small deviations might be a result of modifications and updates in the MPS-ATLAS code compared to the original ATLAS9 at that time.

Summary and conclusion
We presented the structure and extended functionality of the MPS-ATLAS code. The code is based on the ATLAS9 code, but has been extended to allow flexible and faster handling of ODF calculations, as well as emergent flux calculations. This makes the code suitable for radiative transfer calculations along rays from 3D MHD cubes. Furthermore, the atmosphere model calculations were sped up and made more user-friendly. We also improved the equilibrium number density calculations, and included NH photo dissociation opacity. The code-toobservations comparison showed that MPS-ATLAS gives excellent agreement with the observed solar spectrum and Vega spectrum (better than some other widely used codes).
The source code is available on request along with a set of testing input files. A detailed explanation of the input files is given in Appendix G. Furthermore, an online tool for calculating ODFs, stellar models and fluxes will soon be released.
A fine grid of stellar models, fluxes and CLVs is being calculated (Kostogryz, et al. in prep.). This grid will cover the range of effective temperature between 3500K and 9000K in 100K steps, a range of surface gravity from logg = 3.0 to logg = 5.0, and metallcities between -5.0 to 1.5 using very fine steps around solar metallicity of 0.05 dex.
The MPS-ATLAS code has two different radiative solvers implemented. In the original ATLAS9 implementation of the RTE solver by R. L. Kurucz, the scattering part of the source function is found iteratively (hereafter, iterative solver). To minimise the computational time, pre-tabulated matrices for the evaluation of the mean intensity on a fixed optical depth grid are used (for more details see Kurucz 1969Kurucz , 1970. The alternative way to account for the scattering part of the source function is to use a Feautrier method as described in Mihalas & Mihalas (1984) (hereafter, Feautrier solver). In this approach, the second-order transfer equation, derived using Feautrier variables together with upper and lower boundary conditions, is solved. The Feautrier method allows to set the number of viewing angles, µ, which are included in the calculation, to three, four or eight (Lester & Neilson 2008).
The two solvers have a different treatment of the source function at frequencies where strong lines are formed in the upper atmosphere, i.e.the case when the top atmospheric point has τ ν > 0.2. Since for such frequencies the atmosphere is not sufficiently high in order to obtain an accurate solution of the RT and non-LTE effects are important, the Feautrier solver sets the source function and the flux to zero (see orange line in Fig. D.1). This differs from the iterative method that still solves the RT equation in such cases, but makes an approximation by considering a prolongation of the atmosphere with a source function, S , equal to theS (top) value. An example calculation of the highresolution flux for a small wavelength interval obtained using the two different solvers is shown in Fig. D.1. The two solutions are both approximations, but lead to different emergent intensities, as can be seen by looking at the Ca II K line (λ = 393.478 nm). Note, that the different treatment has a negligible effect when ODFs are used for wavelengths longward of 200 nm (see Fig. D.2).
We tested the performance of these two solvers. For that we compare the resulting effective temperature calculated from the emergent flux with the value set in the atmosphere model calculation, and we measure the computation time needed to reach a converged model. We found that, generally, the iterative solver is the most accurate but slower compared to the Feautrier solver. The optimal trade-off between computational time and accuracy is achieved when using the Feautrier solver with three µ angles. Ideally, to obtain consistent results the same solver should be used for the specific emergent intensities as is used for the radiative equilibrium calculations.  by summing the calculated intensity in each sub-bin, weighted by the sub-bin width, w s , as follows where I λ c is the intensity at λ c obtained using the capacities along the atmosphere, κ tot,s,i , which are the sum of the continuous opacity at λ c and the averaged sub-bin opacity of the sub-bin s in each atmosphere point.

Appendix F: Molecular photo-dissociation
ATLAS9 includes photo-dissociation for CH and OH included. This is achieved by using pre-tabulated cross-sections and corresponding partition functions. In addition, we implemented the NH dissociation opacity. Our implementation is, however, slightly different from that of CH and OH. The number densities for NH that are obtained from the equilibrium calculations are multiplied with the cross-sections directly. Since a negligible fraction of NH might be in an excited state, this treatment potentially overestimates the opacity.
To test the effect of NH photo-dissociation opacity on the specific intensities, we calculated the solar atmosphere model in RE with and without the NH photo-dissociation, and subsequently the corresponding emergent intensities. The ratio of the emergent intensities is shown in Fig. F.1. It indicates that the difference is below 0.5% at disc centre and below 0.3% at the limb (µ = 0.05). Thus we conclude that including NH is not sufficient to account for the missing opacity in the UV.

Appendix G: Example input
To run MPS-ATLAS several input files are required. In an overall control file, 'mpsa.control', (shown in Fig. G.1) the user has to specify which module or modules should be executed and give the names of the input files. To set a module for calculations a line with the keyword 'calc_' followed by either 'odf', 'model' or 'flux' has to be included. By adding an 'off' after the keyword the module can be switched off again without deleting the line. In the example in Fig. G.1 the code will only calculate the emergent flux module.
For every module there are three kind of input files that needs to be specified: i) the input file that contains computational flags, ii) a model atmosphere file, and iii) the file that contains the opacity tables, where for the ODF table calculation the name of the output file is specified. The order of the lines in the control file is irrelevant, but to indicate that there are no more lines to read the keyword 'endfile' has to be given.

Appendix G.1: Input files
To start a particular module, it has to be activated in the 'mpsa.control' file. Subsequently, each module requires different settings in the '.input' file, where some settings are common. The order of the keyword lines in the input file is not important, except for the last two lines which have to be 'begin' and 'end', to indicate the end of the input file and start the calculation. Example input files for each module are given in Fig. G.2-Fig. G.4. All three modules have most of the flags in common, while a few are specific for certain modules.

Appendix G.1.1: Common settings
In any of the modules, molecules can be included or excluded by the keywords 'molecules on' or 'molecules off'. If the molecules are 'on', a routine reads the file 'molecules.dat' which should be located in the INPUT folder.
For all modules the chemical composition needs to be specified (as shown in Fig. G.2 -Fig. G.4. The line starting with the keyword 'abundance scale' specifies the metallicity. The number following this keyword scales all elements heavier than helium. Consequently, a single scaling factor accounts for changes in metallicity, so that it is not necessary to change all elements individually when changing the metallicity. All element abundances can be set using the keyword 'abundance change' followed by the elemental number and the number fraction. While for hydrogen and helium the number faction is given as N element /N total , for all other elements it is given in log 10(N element /N total ). For consistency the settings in all input files should be kept the same throughout the calculation of a given star.
For the opacity calculations the continuous opacity sources as listed in Sect. 2.1 are all included by default. If any of them need to be excluded they can be switched off by adding 'opacity off' followed by keyword of the corresponding opacity source. The continuous opacity sources are grouped as follows: 1) H I bound-free transitions (bf) and free-free transitions (ff), keyword 'H1' 2) H + 2 bf and ff, keyword 'H2+' 3) H − bf and ff, keyword 'H-' 4) Rayleigh scattering on H I, keyword 'Hray' 5) He I bf and ff, keyword 'He1' 6) He II bf and ff, keyword 'He2' 7) He − ff, keyword 'He-' 8) Rayleigh scattering on He, keyword 'Heray' 9) -11) various bf and ff transitions of C, N, O, Ne, Mg, Al, Si, Ca, Fe, the molecules CH, OH and NH, and their ions, keywords 'Cool', 'Luke', and 'Hot' 12) electron scattering, keyword 'Elec' 13) Rayleigh scattering on The keywords 'user defined binning' switches between the original standard bin and sub-bin configuration used in Castelli (2005b) if set 'off', and a user-defined configuration if set 'on'. Having set a user defined binning and sub-binning in the input files, individual grids have to be provided in separate input files, where the number of grid points, their intervals, and the number of sub-bins, and their sizes are specified. Then, the code will allocate the requested number of bins and sub-bin. This switch sets the bin and sub-bin configuration for module I, but has to be consistent with the format of the ODF used in module II and III.
The two keywords 'print' and 'punch' control how much information is written out during a run. They are both followed by an integer. For 'print' the following values can be set: -0 for no print -1 for summary tables -2 prints the temperature corrections, and surfaces fluxes -3 prints everything listed as for point 2, and also τ ν , S ν and J ν for each frequency -4 prints the quantities listed under point 2 and all opacities for each frequency For 'punch' there are only four options: -0 nothing is written out -1 writes atmosphere model -2 writes punch 1 and in addition surface fluxes or intensities -5 writes punch 2 and also molecular number densities over partition functions Appendix G.1.2: ODF generation specific settings A specific switch for the ODF calculations can turn filter calculations on or off by using the keyword 'filter' (for more details on filter ODF see Anusha et al. submitted). If this switch is not set, the filter calculations are switched off by default.
Module I has additional keywords: In the default setting at least one ODF table using the micro-turbulence velocity of v turb = 0.0 km/s is calculated. To obtain additional ODF tables the keyword 'velocities values', followed by an integer indicating how many turbulent velocity values should be calculated, can be set. After the integer, the code expects the corresponding number of float numbers setting the micro-turbulence velocities. Moreover, by default the opacities in the sub-bin are averaged using the geometric average. This can be changed to an arithmetic average by setting the keyword combination 'mean arithmetic'.
Appendix G.1.3: Model calculation specific settings For the atmosphere model calculations, the desired stellar parameters, and the Rosseland grid need to be specified. This is controlled by using the keyword 'scale' which has three different options: no letter follows the 'scale' keyword; the code reads in the number of depth points of the model, the starting Rosseland mean optical-depth, τ Ross , in log10, the step-length in log10, the effective temperature, and the surface gravity, logg; if the letter 'x' follows, instead of the step-length, the maximum τ Ross is read if the letter 'b' follows, instead of the number of depth points, the starting τ Ross , the step-length, and the maximum τ Ross is read, before T eff , and log g The keyword 'iterations' has to be set for the model calculations indicating the maximal number of iterations to be executed. Note, that the calculations will be stopped if the model converges before this number is reached. For this module the effective temperature and surface gravity of the starting model, which is used for the initial re-scaling (see Sect. C), need to be specified. Thus, a line starting with the keywords 'starting model' followed by the 'teff' and 'surface gravity' values has to be included.
By default the convective flux and overshoot is turn off. To change that, the keyword 'convection' sets the convection flux calculations 'on' or 'off', and a line with 'overshoot on' turns on the additional convection flux in the overshoot region (Castelli 1996). The mixing-length value is the number following 'convection on', and is usually set between 0.0 and 2.0. Moreover, turbulence can be set on/off, by including the line 'turbulence' followed by four numbers. The turbulent velocity,