Interactive comment on “ ASIS v 1 . 0 : an adaptative solver for the simulation of atmospheric chemistry ” by Daniel Cariolle et al

The manuscript of Cariolle et al. describes the development and testing of a chemical solver suitable for implementation in models of atmospheric chemistry (ASIS). It appears that this solver is not intended to be available to the community due to the restrictive nature of the "code availability" (Section 7). Rather, it appears that the authors have chosen for their own reasons to develop an in-house solution instead of adopting an open source solver package such as KPP, and wish to use this publication to describe its basic features.


Introduction
In chemical transport models (CTMs) or general circulation models (GCMs), the description of atmospheric chemistry has rapidly increased in complexity.Early model developments were devoted to the study of the stratospheric and upper tropospheric compositions focusing on the gasphase reactions that control the ozone distribution.Emphasis has since been put on tropospheric chemistry due to its oxidant properties and its possible impact on climate via the lifetime of several greenhouse gases and the distribution of secondary-formed aerosols.
Large-scale models now include chemical schemes that deal with about a hundred species and with several hundred reactions in gas phase and in heterogeneous phases (solid and liquid).Most of those species undergo transport processes, like advection, diffusion, and convection.As a result the models include the solution of complex coupled systems which cannot be handled in a single operator.In practice the various processes are decomposed in a series of operators that are solved numerically in sequence.For example, the time evolution of the species are first calculated taking into account advection, then diffusion, convection, and so on.Among these processes the evolution of the species due to chemical transformations is a key component of the models.
The models have to solve coupled ordinary differential equation (ODE) systems that describe the adopted chemical mechanism.These ODE systems are of the non-linear form: where C represents the vector of the local species concentrations, P(t, C) and L(t, C) the production and loss term matrices.The stiffness of the systems comes from the wide range of values that can take the production and loss terms.Small values of the loss term correspond to stable species having long lifetimes (e.g., CH 4 , N 2 O, . . . ) whereas large values correspond to radical species (e.g., O( 1 D), OH, Cl, . . . ) with short lifetimes.Typical atmospheric situations lead to species lifetimes ranging from milliseconds to years.Since the other physical processes can change the conditions and compositions of the air masses (i.e., surface emissions, transport at all scales, day-night transitions, etc. . .), the chemical system is often out of chemical equilibrium and the ODE system to be solved can be very stiff.
Adequate algorithms must then be used to deal with the stiffness of the ODE systems and to achieve good accuracy.Existing algorithms vary in formulation and complexity.They can be classified as explicit or implicit schemes (see, e.g., Sandu et al., 1997a), and use single or multistep or multistage methods (Sandu et al., 1997b).Multistage implicit algorithms based on Gear (Hindmarsh, 1980) and Rosenbrock (Hairer and Wanner, 1991) formulations are the most accurate but require a significant amount of computation, which limits their use in comprehensive atmospheric chemistry models.To reduce the computational cost the solver described in the present article is based on a one-step implicit algorithm.Its characteristics are detailed in the next sections along with comparisons with other implicit schemes.
For atmospheric applications some numerical properties of those algorithms should be particularly sought for the following: i. Mass conservation: atmospheric models are often integrated for long-term simulations (up to several decades for global climate simulations), and small trends and anomalies are investigated.Any bias or trend in the atmospheric composition due to numerical algorithms must therefore be avoided.It is therefore essential that the algorithms chosen to solve the chemical systems preserve mass.All the atoms or elementary groups of atoms (e.g., nitrogen oxides) must be conserved.
ii. Accuracy: it is of course always desirable to obtain a numerical solution that is as accurate as possible, although the uncertainties associated with the other operators and the fact that they are integrated successively in time introduce a significant degree of inaccuracy.This leads also to transient evolutions in the chemical system, especially for the short-lived radicals, that have no real physical basis.It is not always necessary to obtain a very accurate numerical solution during those transient evolutions if they do not last long and have little impact on the solutions for the other longer lived species.The key point is to design an algorithm where the accuracy can be chosen a priori by the user and controlled during the course of the numerical integration.
iii.Positivity: it is highly desirable to maintain positivity of the concentrations.Otherwise instability might arise when coupled to other operators dealing with advection or convection of the minor species.Some algorithms maintain positivity of the solution by construction, others introduce clipping of the negative values at the expense of local mass conservation.Negative values can be tolerated if they are small and transient, and if they have little impact on the algorithms used to account for the other physical processes.
iv. Adaptability and flexibility: the adopted solver should cope with a variety of chemical mechanisms with the possibility to easily add or remove species and reactions.It is also desirable to give to the user a minimum of free parameters to tune.The solver should also run efficiently on a large variety of computers without having to rewrite large parts of the code.This can be obtained with extensive use of mathematical libraries that are often optimized for the computer being used.
This article describes a solver for the simulation of gas-phase atmospheric chemistry, the adaptive semi-implicit scheme (ASIS), that has most of the desirable properties discussed above.Section 2 gives the basic formulation of the scheme and discusses its characteristics in comparison with other solvers currently used within atmospheric models.Section 3 gives results from box-model simulations and comparisons with other state-of-the-art algorithms, and Sects.4 and 5 detail the implementation of ASIS within the MOCAGE CTM (Michou and Peuch, 2002;Josse at al., 2004;Teyssèdre et al., 2007) and the GCM of planet Mars (Lefèvre et al., 2004) of the Laboratoire de Météorologie Dynamique (LMD).Possible future extensions of the solver are discussed in the final section.
2 ASIS: description of the chemical solver

Implicit discretization and numerical methods
The integration of Eq. ( 1) cannot be done using a simple one-step explicit scheme with the left-hand-side terms evaluated at time t, as numerical stability would require the use of time steps lower than the shortest species lifetime.Since some radicals have lifetimes lower than a few milliseconds in the atmosphere, too many iterations would be required to obtain simulations for 100 days or more with affordable computer time.Several explicit methods have been developed to address this issue which are based on classification of the species according to their lifetimes.For instance, with the QSSA method (Hesstvedt et al., 1978) the fast species with lifetimes much lower than the time step are often assumed to be at equilibrium (C = P/L) the intermediate species are obtained using an exponential solution of Eq. ( 1), and the long-lived species are computed using the simple explicit solution.Other explicit schemes gain in accuracy with the use of multistep algorithms with predictor-corrector evaluations of the concentration at time t +δt, for instance the CHEMEQ solver of Young and Boris (1997) with subsequent developments by Mott et al. (2000).Limitations of these explicit schemes are that they often do not conserve mass and that the choice of species classifications is somewhat arbitrary.Mass conservation can be improved using the technique of "species lumping" where additional equations are introduced for linear combinations of species concentrations to reduce the stiffness or enforce conservation for a chemical family.The drawback of those approaches is that the algorithm becomes problem-dependent and requires a very good knowledge of the chemical system, especially when updating the constant rates or the list of reacting species.One possibility to increase the time step is to treat part of the right-hand side of Eq. ( 1) implicitly, for instance keeping the evaluations of P and L at time t but C at time t + δt: where C t+1 is the concentration vector of the species at time t + δt.The second term of the left hand side of this equation forms a diagonal matrix so the numerical solution of Eq. ( 2) is straightforward.With this discretization the numerical solution is positive and unconditionally stable.The mass conservation is however not maintained due to the fact that for a given reaction between two species the value of the associated tendency is different for each species.
One way to alleviate this problem is to discretize Eq. (2) fully implicitly in time using the simple Euler-backward (EB) method: Solution of Eq. ( 3) requires an evaluation of the terms L(t + 1, C) and P(t + 1, C) that can be obtained using C t+1 from the solution of Eq. ( 2).In practice Eq. ( 3) is solved iteratively with successive evaluations of C t+1 , for instance using the iterative Newton method.A correcting term to the iterative solution can also be added to increase the accuracy (Stott and Harwood, 1993;Carver and Stott, 2000).
Still, the mass conservation can only be obtained if a good convergence of the solution is reached and additional constraints, such as species lumping or equilibrium assumptions for the shorter-lived species, are often used to increase accuracy and to speed up convergence.Schemes of this type are, for example, used within the MOZART model (Emmons et al., 2010), the ECHAM-HAMMOZ model (Pozzoli et al., 2008), the TM5 model (Huijnen et al., 2010), the UKCA climate-composition model (O'Connor et al., 2014), and the MOCAGE CTM.
The implicit methods described above to solve the ODE chemical system are all one time step: only concentrations at time t are used to evaluate the concentrations at time t + δt.Although the numerical stiff ODE field is largely developed and accurate ODE solvers are available and have been used for atmospheric chemistry problems, many of them involve multiple steps or stages.Consequently, several evaluations of C at various past or intermediate time steps are used to obtain the concentration at time t + δt.A direct extension of the simple EB method is to use a higher-order backward differentiation formula (BDF) to solve Eq. ( 1).Based on that approach, Verwer (1994) has developed the TWOSTEP atmospheric chemical solver, which uses a second-order BDF formula combined with a Gauss-Seidel iteration technique to solve the resulting implicit system.This solver can be very efficient but it is not naturally mass conserving.It is, for example, implemented in the CHIMERE model (Menut et al., 2013).
Mass-conserving, multistep or multistage, and high-order accurate implicit methods exist to solve the ODE stiff system.Among the methods based on BDF, Gear's predictorcorrector method has been adapted to atmospheric chemical systems, for example, the SMVGEAR code (Jacobson and Turco, 1994) implemented in the GEOSCHEM CTM (Bey et al., 2001).More recently, the Rosenbrock's method (Rosenbrock, 1963) is becoming widely used in atmospheric chemistry modeling (Sandu et al., 1997b) despite the fact that its computational cost is still rather high compared to approaches based on low-order BDF methods.The implementation in chemical models of Rosenbrock's and other highorder methods has been eased by the development of the Kinetic PreProcessor (KPP) by Sandu and Sander (2006), which allows the choice of an integration method and generates the adequate codes accordingly.
When the chemical scheme involves more than 100 species and over 200 reactions, the implicit multistage methods are still computationally expensive, especially if they are to be used within global 3-D models with horizontal resolutions on the order of 1 • × 1 • , with several tens of vertical levels and for simulations lasting for several years.The increase of the computational cost comes from the need to solve at each stage a linear system on the order of the number of species, and this cost varies non-linearly (often quadratically; see, for example, Golub and Van Loan, 2013) with the number of species.

Formulation of the ASIS solver
The approach adopted for ASIS is to restrict the algorithm to a single implicit step combined with a specific evaluation of the Jacobian matrix of the chemical fluxes, J = f (C) = ∂f/∂C.

D. Cariolle et al.: ASIS: an adaptive solver
The starting point comes from the decomposition of chemical tendencies in three terms: The first term of the right-hand side corresponds to the chemical productions or destructions due to first-order reaction rates with constants K.The second term arises from thermal decompositions and/or photodissociations of species with a rate D, and the last term accounts for external tendencies that come from other physical processes than chemistry.For example, the surface emissions affecting the lowest levels of the model will result in species tendencies F.
The time discretization of Eq. ( 4) is then performed with a semi-implicit scheme for the first term adapted to each reaction and time step and an implicit discretization for the second one, and the external tendencies are assumed to be constant over the time step δt: with 0 ≤ t l,m ≤ 1. Equation ( 5) can be recast with terms containing species concentrations at time t + 1 on the left-hand side and the others on the right-hand side.
that can be reformulated in a matrix form: with the matrix M, an approximation of the Jacobian J , containing species concentrations at time t and values of l,m also evaluated at time t.
Compared to other one-step semi-implicit schemes like SIS (i.e., Ramarosson et al., 1994), one specificity of our scheme lies in the evaluation of t l,m .Let us consider the system of a single reaction between species C l and C m with a reaction rate constant K l,m .If the initial values of the concentrations are equal (C 0 l = C 0 m = C 0 ), the exact solution of the system gives a hyperbolic decay for the concentrations: This solution is obtained exactly using the discretization given by Eq. ( 6), with l,m = 1/2.If C 0 l C 0 m , the evolution of the lowest concentration C m shows a quasiexponential decay with an e-folding time τ = 1/(K l,m C 0 l ) while the concentration C l reaches its steady-state value C 0 l − C 0 m that does not depart strongly from its initial value.In that case, in order to maximize the time step and to increase the stability of the scheme, there is an advantage in treating the evolution of the shorter-lived species C m as implicitly as possible by giving more weight to the term C t l C t+1 m in Eq. ( 6).This is obtained if l,m tends towards 1.Those simple considerations lead us to introduce the following function for l,m that depends on the concentrations at time t: with β ≥ 1.With this formulation the value of l,m has the required properties: as a function of the concentrations.Large values of β favor the implicit treatment for the lowest concentrations.Positive values of β lower than 1 could also be used but they may not discriminate enough the treatment of species according to their concentrations.For the situations studied in this paper the numerical simulations did not show a large sensitivity to this parameter, which was fixed to 1 hereafter.Furthermore, the use of Eq. ( 9) to calculate t l,m and evaluate M, the approximate Jacobian matrix, gives interesting properties to our scheme: -The oscillations from odd to even time steps that can appear in the numerical solution of Eq. ( 6) when the semi-implicit scheme is centered and symmetrical (i.e., Suhre and Rosset, 1994), as would be if the fixed value t l,m = 1/2 was adopted, is damped with the evaluation of t l,m by Eq. ( 9).
-Since the largest terms contributing to the evolution of the shortest-lived species are treated implicitly, the system increases in stability.Larger time steps can be used and positive values for the concentrations are more easily preserved.
-All the species are treated in the same manner without any a priori considerations on lifetimes or abundances.For instance in the case of the Earth composition, O 2 is treated like the other species even if its chemical sources and sinks are negligible.Since the concentration of O 2 is much larger than the other species concentrations, any species reacting with O 2 will be treated implicitly.This is the case, for example, with atomic oxygen O reacting with O 2 to form O 3 .The corresponding term in Eq. ( 6) for the O tendency will be O The option to treat all the species in the same manner simplifies the programming of the scheme and allows the solver to be easily adapted to various chemical systems.An example is given in Sect. 5 with the simulation of the atmospheric composition of the planet Mars.
Once the matrix M is evaluated and the time step is determined (see the next section), the solver computes the solution to the system of linearized Eq. ( 7), which becomes a possible computational bottleneck.Our approach is to use standard methods and well-optimized software libraries.
Our baseline option is to use the direct solver DGESV of the Lapack library that solves system (6) by lower-upper (LU) decomposition.Therefore no extra specific routine associated with the chemical mechanism is needed and the optimization on the computer used is left to the implementation of the Lapack library.As reported below, this option works well and gives accurate results even for comprehensive mechanisms involving hundreds of species or more.
To reduce the computational cost, other options for the solution of the linear system have been investigated.Two iterative solvers have been tested.The first one is an implementation of the Gauss-Seidel algorithm.This algorithm has been used with success to solve stiff systems from chemical kinetics (Verwer, 1994;Menut et al., 2013).For the cases studied in the following sections, the Gauss-Seidel algorithm was found to be efficient with a good rate of convergence in most cases.Although in specific situations where the system is largely driven out of equilibrium, for instance during daynight transitions and for large surface emissions, the number of iterations could increase by 1 order of magnitude to obtain the required accuracy.
A second iterative algorithm has been implemented, the generalized minimal residual method (GMRES).The method approximates the solution by a vector in a Krylov subspace with minimal residual norm.The Arnoldi iteration algorithm is used to find this vector.The GMRES method was developed by Saad and Schultz (1986) and further described by Saad (2003).In order to accelerate the convergence, preconditioning techniques are used.An efficient technique was obtained by introducing the matrix B using the lower triangular part of the A = I − Mδt matrix to compute an approximation of A −1 and apply GMRES to the solution of the right-preconditioned linear system ABC * = C t + Fδt where C t+1 = BC * .For the implementation discussed hereafter, the GMRES method needs fewer iterations than the Gauss-Seidel one, especially in situations where the Gauss-Seidel algorithm shows slower convergence, and was found to speed up the computation by at least a factor of 2 compared to the DGESV implementation.

Time stepping
Since the time discretization adopted to solve system (7) is first-order accurate, the choice of the time step δt is important to obtain a solution with a desired accuracy.In our applications the evolution of the species over rather large time intervals t is required.The value t is determined by other physical processes than chemistry, for instance advection, convection, or vertical diffusion, and is often too large to be used directly to solve Eq. ( 7) without encountering nu-merical instabilities and loss of accuracy.For example, in the 3-D model results discussed in Sect.4, the time interval t = 15 min is determined by horizontal advection whereas the chemical time step has to be decreased to a few seconds in situations where the chemical state is driven far from a quasi-steady-state.
Therefore a variable stepsize strategy has to be implemented with the time interval t divided in n successive integrations of the chemical system with time steps δt n .
The choice of δt n is made iteratively using a strategy similar to the one described by Verwer (1994).First a local error indicator E is computed: with where γ = δt n /δt k+1 , δt k+1 is a first-guess time step, C k+1 m is the concentration of the species m at the iteration k + 1, and ATOL and RTOL are absolute and relative error tolerance.C k+1 m is evaluated using Eq. ( 2) with C n m as initial concentration and the time step δt k+1 .E k+1 depends on the curvature of the solution, a measure of the departure of the solution from linearity.If E k+1 ≤ 1, the time step δt k+1 is adopted (δt n+1 = δt k+1 ), otherwise a new time step is estimated by the following: Then a new value C k+2 m is evaluated followed by the computation of E k+2 , and so on until convergence.In practice the convergence is obtained within a few iterations, less than 5 in the cases reported thereafter.Those iterations have a low computational cost because the solution of Eq. ( 2) at each iteration involves only diagonal matrices.Once the value of δt n+1 is determined, the concentration C n+1 m is obtained by the solution of Eq. ( 7).
For the first iteration species, concentrations at two consecutive times and a first-guess time step are needed.To avoid storing concentrations at consecutive times we assume that at the beginning of the iterative process the system is in a steady state, C n m = C n−1 m in Eq. ( 10), and the first-guess time step is set to its largest possible value t.To secure the iterative process a minimum time step, δt min , is also prescribed in order to limit the number of iterations.The value of this minimum time step is left to the user who has to choose a value consistent with the error tolerance parameters.

Tests and validation
To validate and evaluate the performances of ASIS and the associated numerical codes, several case studies have been www.geosci-model-dev.net/10/1467/2017/Geosci.Model Dev., 10, 1467-1485, 2017 used.All the cases reported in this section are based on the RACMOBUS chemical scheme used within the MOCAGE CTM.RACMOBUS is a combination of the REPROBUS scheme adapted to the stratosphere and the free troposphere (Lefèvre et al., 1994) and the RACM scheme (Stockwell et al., 1997) that treats the urban polluted earth atmosphere with the addition of volatile organic compounds (VOCs) and their degradation products.Table 1 lists the chemical species taken into account; the overall scheme includes about 120 species linked by 200 gas-phase reactions and photodissociations.The photodissociation rates are calculated every 15 min using the tropospheric ultraviolet and visible (TUV) radiation model version 5.2 (Madronich and Flocke, 1998) for conditions corresponding to the equinox at 30 • latitude.Two test cases are used to evaluate the accuracy and performance of the ASIS scheme.The first one is based on the FLUX test case described by Crassier et al. (2000).It corresponds to a ground-level situation in a polluted urban area.The list of species and fluxes emitted at the surface is given in Table 2.The emissions are injected in a boundary layer with a 2000 m constant thickness weighted by an emission factor of 0.6.This leads to a constant tendency F in Eq. ( 4) for the emitted species.The initial concentrations are given in Table 3, the atmospheric temperature is set to 298 K, and the ground pressure is 1000 hPa.
The second case, STRATO, is representative of situations encountered in the middle stratosphere.The initial concentrations for this case are given in Table 3.The atmospheric temperature is 215 K and the pressure is 50 hPa.For both cases the integration starts at midnight and stops 24 h after, and the photodissociation rates are updated every 15 min.

The FLUX case
To assess the performances of ASIS, two reference simulations have been obtained for the FLUX case using Rosenbrock's and Gear's BDF solvers (referred to hereafter as R1 and G1).Those solvers use the "ode23s" and "ode15s" codes, respectively, from the Matlab ODE suite (Shampine and Reichelt, 1997;Ashino et al., 2000).For the Rosenbrock's scheme a three-stage algorithm is used and the simulations are third-order accurate.For the Gear's scheme the third-order accurate option was also chosen.In the most accurate simulations, R1 and G1 (see Table 4), the relative tolerance RTOL is set to 0.001 and the absolute tolerance ATOL equals 10 4 molecules cm −3 for all species.With the Rosenbock's solver additional simulations with higher tolerance values have been performed, R2 with RTOL = 0.01 and R3 with RTOL = 0.025.
The same FLUX case is integrated using the ASIS solver.In a first simulation noted A1, ASIS uses a RTOL value of 0.001 and a minimum time step of 1 s.For the solution of the linear system associated with ASIS, the DGESV code of the Lapack library is used.To compare with the Rosenbrock's solver a second simulation A2 has been obtained with the same settings as A1 but with a higher relative tolerance value of 0.01, and a third one A3 with a tolerance value of 0.025.For all experiments the FLUX case is integrated over 24 h.The settings used in the overall simulations are given in Table 4.
Figures 1 and 2 show the evolution of some key species for each experiment and the relative differences from the R1 experiment.Those results are representative of all of the species.As expected, the R1, G1, and A1 simulations give very close results.R1 and G1 show relative differences below 0.1 %, consistent with the value chosen for RTOL.A1 results are comparable with differences in the 0.1-0.2% range, except at the beginning of the simulation when the chemical state is out of equilibrium and during day-night transitions.In those situations the differences between A1 and R1 or G1 can reach 0.5 %.As expected from the choice of a higher value for RTOL, the A2 experiment shows less accuracy but is still in the range of 0.5 % compared to the other experiments.The A3 experiment has differences below 2 % with the other experiments.For most of the atmospheric simulations an accuracy below 1 % is sufficient for the longest-lived species, and even larger values are acceptable for short-lived species if they are transient, given the uncertainties in the representation of the other processes and the inaccuracies introduced by their solution by a series of successive operators.
The efficiency of ASIS can be first evaluated by comparison of the mean time steps (Table 4).For simulations R1 and G1 the mean time steps are between 25 and 40 s.Since ASIS uses a first-order scheme to maintain good accuracy, the mean time step is lowered, in the 5 s range for the A1 experiment.However, ASIS is a one-stage scheme (only one linear system is solved by time step) compared to R1 and G1 that need three or more stages.The amount of computation is therefore comparable.When the relative tolerance is increased, the mean time step of ASIS increases.For the A2 experiment it is 25 s, identical to R1, and up to 49 s for A3.Since for most atmospheric simulations a relative tolerance of 0.01 to 0.025 seems to be sufficient, the ASIS solver gives acceptable solutions with less computation than the higherorder schemes.
The efficiency of ASIS in terms of CPU time has been evaluated within the Matlab environment.Table 4 gives the ratio of CPU time used for each simulation relative to the R1 case.The ode23s code used to run the R cases needs the implementation of two subroutines, one that computes the species tendencies and one that gives the Jacobian of the system.If the latter is not provided, the ode23s code computes an approximation of the Jacobian by differentiation and the CPU cost increases by a factor of 2 to 10.As can be seen from Table 4, the CPU cost of ASIS is comparable to or lower than the ode23s cost for relative tolerance values larger than 0.01.
An important point to mention is that within the Matlab environment the CPU cost does not come from the linear algebra parts of the algorithms but from the evaluation of tendencies and Jacobian matrices.Therefore it is very dependent upon the chemical system and the details of the programming of the associated subroutines.The situation is quite different within the Fortran environment.With the Fortran version of ASIS the CPU cost for the calculation of the approximated Jacobian (the matrix M of Eq. 7) is negligible compared to the linear algebra computations.This is because the compiler efficiently handles the associated subroutine (fill_matrix, see Sect.7) that contains frequent indirect addressing.It is not possible to evaluate if this is also the case with all the codes based on Rosenbrock's algorithm, but if it is so, ASIS should perform well when the mean time steps are comparable since it needs fewer linear algebra computations.
For the A experiments, ASIS uses the DGESV code for the solution of the linear systems.To save computational time two iterative solvers have been tested, one using the Gauss-Seidel algorithm, the other the GMRES method.Both solvers used the same criterion for convergence (tolerance for convergence set to 10 −14 ).For the GMRES method the preconditioning technique described in Sect.2.1 is implemented.With those settings the experiment A2 has been repeated.
The results are practically identical to the solution obtained using the DGESV code; the differences between the solutions are below 0.02 % for all the species concentrations.The simulation with the Gauss-Seidel algorithm shows good efficiency in terms of mean number of iterations, but requires 6 to 10 times more iterations when the system is driven out of equilibrium during day-night transitions.Using GMRES was found to be more stable and efficient, with less than 10 iterations needed to solve the linear systems and half as much computational time (using the Fortran version of the code) compared to the simulation using DGESV.
From the simulations of this FLUX case, which is rather representative of situations encountered in polluted earth boundary layers, it can be concluded that the ASIS solver performs well compared to higher-order schemes when moderate accuracy is required.Apart from tolerance parameters and the choice of a minimum time step, no specific tuning is required.The one-step implicit scheme gains in efficiency when coupled to the GMRES iterative solver used for the solution of the linear systems.

The STRATO case
The STRATO case differs from the FLUX case in the dominant chemical regimes involved.In the FLUX case the VOC decomposition during day and night dominates the system.With the STRATO case the chemistry is dominated by NO x , HO x , Cl x catalytic cycles, and the ozone content.The stiffness of the system is less stringent, and rapid variations in the www.geosci-model-dev.net/10/1467/2017/Geosci.Model Dev., 10, 1467-1485, 2017  The left column shows the volume mixing ratios and the right one the differences in percentage (%) relative to the R1 experiment.The color code is the following: blue for G1, orange for A1, red for A2, and purple for A3.concentrations of the species are tightly linked to the variations of the insulation at sunrise and sunset.For this case, two simulations have been performed.The first one, RS1, uses the Rosenbrock's algorithm with settings similar to experiment R1.For the second one, AS2, the ASIS solver is used with settings similar to experiment A2 and with the iterative linear solver GMRES.The two simulations show results consistent with the findings for the FLUX case.The mean time steps are very similar for both experiments, 49 s for RS2 and 41.4 s for AS2.As expected, the time step decreases in ASIS at sunrise and sunset when the stiffness of the system is at maximum. Figure 3 shows the number of time steps for every 15 min interval for the 24 h simulation of experiment AS2.Apart from the very beginning of the simu- lation that starts in a situation out of equilibrium, the largest values are found at sunrise and sunset in the 150-200 range.It corresponds to time steps of about 4.5 s.
In terms of accuracy, the AS2 experiment gives results that depart less than 1 % compared to RS1.This is consistent with the chosen value of 0.01 for RTOL. Figure 4 shows that the lowest accuracy is found at sunrise and sunset when the short-lived radical species have the largest variations.Those transition situations are the most difficult because not only must the accuracy be maintained but spurious numerical oscillations must be avoided.ASIS performs well here and adapts its time step automatically to reach the required accuracy.The numerical treatment adopted to calculate an approximation of the Jacobian (Eq. 9 with β = 1) contributes greatly to the reduction of numerical oscillations without significant degradation of the accuracy of the solution.
Equally, the approximations in the Jacobian are efficient to prevent the development of negative mixing ratios.In the two cases FLUX and STRATO, we did not encounter any significant (larger than ATOL) negative values during the course of the simulation, and all the concentrations remain positive at the end of the 15 min intervals before the photodissociation rates are updated.
In summary, the results of the two test cases confirm the properties targeted in the design of ASIS.At the moderate accuracy required for atmospheric simulations, the ASIS solver compares well with higher-order schemes, and limits the computational cost while assuring mass conservation.The next sections illustrate how it performs in more realistic situations with implementations in state-of-the-art global CTMs for Earth and Mars atmospheres.

Implementation within the MOCAGE model
For this study we have used the global version of the MOCAGE CTM with a horizontal resolution of 2 • × 2 • and 47 levels in the vertical from ground to 5 hPa.The chemical scheme is RACMOBUS, identical to the one used for the test cases of Sect. 3. In addition to chemistry and transport by the large-scale winds and by convection, the model includes the main processes that contribute to the sources and sinks of the species: surface emissions, scavenging by rain, and dry and wet depositions.The time steps for these processes is 15 min; photodissociation and chemical rate constants are updated at the same frequency.We report here simulations over 3 months from the beginning of August to the end of October 2011.Wind and temperature fields come from the operational weather analyses of the ECMWF.They are updated every 3 h and linearly interpolated in between these time intervals.
The reference simulation (referred to hereafter as MR) uses the original solver for chemistry, an iterative semiimplicit scheme with assumptions of equilibrium for shortlived species and species lumping for NO x and Cl x families.The chemical time step varies with altitude but is kept constant during the model integration.It increases from 20 s in the planetary boundary layer to 15 min in the stratosphere.
The simulation with the ASIS solver, MA, uses the same configuration for MOCAGE as MR except that the original chemical solver is replaced by ASIS with settings similar to experiment A3: RTOL = 0.025, ATOL = 10 4 molecules cm −3 , a minimum time step of 5 s, and the GMRES solver for the solution of the linear systems.
The characteristics of the ASIS functioning implemented within MOCAGE can be first examined by the diagnostic of the number of sub-time-steps for chemistry.Figure 5 shows this number for three different levels for a date corresponding to 15 September at midday.In the midstratosphere, at 50 hPa, the number of sub-time-steps varies in accordance with what was found for the STRATO test case.At midday or midnight the chemical system is in quasi-steady-state and this number is small, below 3. Close to the terminators, this number increases up to 40-60, highlighting the change of regime of the chemical system when the photodissociation is activated or deactivated.During these transition phases the stiffness of the system increases and the sub-time-steps decrease to maintain the required accuracy.Also barely noticeable is an increase of the number of sub-time-steps over the Antarctic coast at the edge of the polar vortex.In these regions the heterogeneous reactions acting at the surface of polar clouds are activated, driving the concentrations of the chlorine species out of equilibrium.It leads to a reduction of the sub-timesteps to cope with the rapid variations of the chemical composition of the air masses.
In the middle troposphere the same behavior is encountered near the terminators, with a tendency to maintain reduced sub-time-steps during longer periods after sunrise or before sunset (Fig. 5).An increase of the number of subtime-steps is also encountered over the African continent at low latitudes.Those regions are prone to convective activity, and injection of species by convection is activated, leaving air masses far from chemical steady state.Since the chemical evolution of the species is calculated after the transport processes, ASIS starts with a situation far from a chemical equilibrium and the number of sub-time-steps increases.
At the surface, Fig. 5 shows the same characteristics as in the midtroposphere with an increase of the number of sub-time-steps at the terminators and over the continents.Over the continents the surface emissions play a larger role than convection in destabilizing the chemical system.Within MOCAGE the emissions are calculated according to inven- tories and deposited in the boundary layer.This is treated as an isolated process that changes the concentrations.As a result ASIS starts with situations out of chemical equilibrium and adopts small sub-time-steps, about 20 s compared to 60 s over the oceans.
Except for noticeable cases that are discussed hereafter, the species distributions of the simulation with the ASIS solver (referred to hereafter as MA) are close to those obtained in MR.As an illustration, Fig. 6 shows the zonally averaged distributions of O 3 , CO, OH, and HNO 3 for the month of September.In most altitudes the differences are below 10 %, with the largest differences in Southern Hemisphere high latitudes in the lower stratosphere.Similar differences are found for the other species, except for the NO x species in the lower troposphere and the chlorine species in the high latitudes of the Southern Hemisphere during the formation of the stratospheric ozone hole.
In the lower troposphere, examination of the code of the MR simulation reveals that approximations and steady-state assumptions are made for the computation of the nighttime NO 2 /NO 3 /N 2 O 5 system.These approximations are valid in the stratosphere but fail in the lower troposphere where the pressure and temperature are larger.As a result the MR simulation underevaluates the concentrations of NO 3 and N 2 O 5 .
Figure 7 shows, for example, the distributions of NO 2 , HNO 3 , and N 2 O 5 at the lower level near the surface averaged over the month of October.In the region of surface emissions, the MR simulation strongly underestimates the N 2 O 5 concentrations.Since the chemical scheme adopted for the present simulations does not include the formation of HNO 3 by the hydrolysis reaction of N 2 O 5 on aerosol surface (Dentener and Crutzen, 1993), it does not have a major influence on the other species.Nevertheless, the maximum values for HNO 3 and NO 2 are larger in the MA simulation than in the MR simulation.
Another significant difference between MR and MA is found in the simulation of the ClO x system in the lower stratosphere at high Southern Hemisphere latitudes.In late August and early September, the solar radiation comes back at high latitudes and the lower stratospheric O 3 is destroyed by catalytic cycles involving chlorine radicals (Solomon, 1999).The chlorine radical concentrations are enhanced by the heterogeneous reactions on polar stratospheric clouds' (PSCs) surface that convert HCl and ClONO 2 into Cl 2 that is photodissociated to form the chlorine radicals.In addition, the catalytic destruction of O 3 also involves the bromine species.
In the air masses prone to heterogeneous reactions on PSC, the composition changes rapidly at sunrise and non-linear processes, like the formation of Cl 2 O 2 , a key species for the O 3 destruction, play a major role.As a result the chemical system is very stiff and the ASIS solver diminishes the chemical time step to a few seconds to maintain good accuracy.In these transient situations the original code in MR does not change its settings and a fixed time step of 15 min is used.
As a result, the MR simulation produces a much more pronounced ozone depletion over Antarctica than the MA simulation.MR calculates ozone column contents as low as 100 Dobson Units (DU), whereas the MA simulation maintains values in the range of 150 DU.This is well illustrated in Fig. 8, which shows the evolution of the total ozone columns over two Antarctic stations, Dumont d'Urville and Dome C. For these two stations the measurements done by SAOZ instruments (Pommereau and Goutail, 1988) at sunrise and sunset are also presented (data available at http: //saoz.obs.uvsq.fr/).
Starting around Julian day 220, the MR and MA simulations start to diverge.Over Dumont d'Urville, the station that sees first the return of the sunlight, the ozone decrease is about 50 % larger in the MR simulation than for MA.By day 260 the ozone column is just above 150 DU whereas it is in the 200 DU range in the MA simulation.Clearly the MA simulation is in better agreement with the SAOZ measurements.
The same behavior is seen for the Dôme C station.The ozone depletion starts slightly later, around day 240.In the MR simulation the depletion is very pronounced and the ozone column diminishes rapidly in a few days from 240 to 150 DU, and further decreases at a slower rate to reach a minimum of 100 DU at day 260.The MA simulation shows a more continuous decrease from day 240 to 260, with an ozone column reaching a minimum of 150 DU.The MA simulation is here again in very good agreement with the SAOZ observations.
Implementation of the ASIS solver within MOCAGE has thus revealed two weaknesses of the original model.One problem is in a limitation on the validity of assumptions made to compute the night-time distribution of the NO x species.It can be solved by adequate coding.The other one is a lack of accuracy in the solution of the chemical system in specific situations in the lower stratosphere.This can certainly be avoided by a drastic reduction of the time step, but it would need the implementation of a time-varying time step strategy somewhat similar to the one adopted for ASIS.
Clearly the implementation of ASIS within MOCAGE is very beneficial to the model simulations and increases the confidence on the model results.In addition, further evolution of the model with adoption of different chemical schemes or addition of new reactions is very easy with ASIS.
There is, however, a price to pay in terms of computer time.Overall the MA simulation takes 4.7 times more computational time than the MR simulation.This number could certainly be decreased by further tuning of the parameters of the solver, RTOL, ATOL, and δt min , and maybe also by the use of the Gauss-Seidel algorithm instead of GMRES in situations where the solution of the linear system converges easily.
Our experience with ASIS shows that, since various processes are computed by a series of operators, the solver starts new time steps with situations often out of chemical equilibrium and must use small sub-time-steps.To alleviate this, one possibility is that tendencies from these operators would be computed and stored rather than used to update the species concentrations.The tendencies can then be used to solve the system though their introduction in the term F of Eq. ( 7).We have tested this option for the species emissions at the surface and found that the number of sub-time-steps is decreased by a factor of 2 in the lower troposphere.It remains to be seen if other processes can be treated that way.Emissions are the most straightforward because the resulting tendencies are positive and cannot lead to the calculation of negative concentrations.
Another issue lies in the parallelization of the computations.In the reference simulation the computational cost is equal for each grid point at a given level and good parallelization is obtained with an equally spaced latitudinal band decomposition (and use of openMP directives).When ASIS is used the computational cost in each grid point depends on the state of the chemical system.As illustrated in Fig. 5, in the stratosphere and upper troposphere more computer time is needed near the terminators and in case of PSC-induced chemistry.In the lower troposphere more computer time is spent in grid points influenced by surface emissions, and convective and boundary-layer transport processes.A speedup of 15 was, however, obtained for the MA simulation on our cluster computer (using 1 node and 16 cores of our BULL computer) with a decomposition that groups more longitudes in the Southern Hemisphere than in the Northern Hemisphere, near the poles.But further tuning would be required if more nodes are to be used.This tuning could vary with season and additional parallelization could be introduced with domain decomposition on the vertical.

Implementation within the LMD Mars model
To illustrate the versatility of the ASIS solver, we present results of the implementation of ASIS in the LMD Mars model with photochemistry (Lefèvre et al., 2004).This Mars GCM describes the evolution of 19 species (Table 5) by means of 54 chemical or photolytic reactions.The bulk atmosphere of Mars is composed of 95 % CO 2 with trace amounts of H 2 O.As a result, the only processes that initiate Martian photochemistry are the photolysis of CO 2 and H 2 O by ultraviolet solar light.Therefore, the photochemistry of the lower at- mosphere of Mars can be summarized by the interactions between the oxygenated species O( 1 D), O, and O 3 produced by CO 2 photolysis and the hydrogen radicals H, OH, and HO 2 produced by H 2 O photolysis.These processes are similar to those occurring in the Earth's mesosphere, with comparable conditions of pressure and temperature.
In the standard version described in Lefèvre at al. (2004), the LMD GCM with photochemistry uses the EB method expressed in Eq. (3) to solve its chemical system.As mentioned earlier, this method is positive, stable, and can be computationally effective but does not maintain mass conservation.Iterative evaluations of C t+1 are performed in the lower atmosphere of Mars to reduce this problem.In the Mars thermosphere, another option is used in the LMD model, which consists of shortening the time step δt according to the species with the shortest lifetime (González-Galindo et al., 2009).In both cases, species lumping and assumptions of photo-  chemical equilibrium are used to increase accuracy and avoid very small time steps.However, conditions of photochemical equilibrium change at night and are also very dependent on altitude.For instance, on Mars, the lifetimes of O( 3 P) and H vary between less than 1 s near the surface and several years at 100 km.Such stark variation prevents the assumption of photochemical equilibrium or use of Eq. ( 3) throughout the atmosphere.Thus, despite its apparent simplicity, the EB method may complicate the problem by requiring dif-ferent treatments for specific species or specific parts of the atmosphere.
Figure 9 compares the results obtained with the EB and ASIS solvers applied to a box-model version of the LMD Mars model.The atmospheric pressure (temperature) is 5.4 hPa (212 K) at the surface and 0.2 hPa (140 K) at 30 km.In both cases the integration starts at noon, and stops after one Martian solar day of 24 h 40 mn.The photodissociation rates are calculated every 15 min using the TUV radiation model adapted to Mars.The time step of the EB solver is fixed to δt = 7.5 min as done in the Mars GCM.ASIS uses the variable stepsize strategy described in Sect.2.3, bounded by a maximum value of 15 min and the minimum time step of 10 s.RTOL is fixed to 0.05 and the ATOL density corresponds to a volume mixing ratio of 10 ppt.ATOL is therefore variable with altitude.The solution of the linear systems associated with ASIS is done using the DGESV direct solver.We found that these settings were adequate to reach a satisfying compromise between accuracy and computing time.
At the surface, Fig. 9 shows that the ASIS solver calculates an O 3 mixing ratio that is lower by 3 to 6 % compared to the EB solver.This difference is related to the lack of accuracy in the treatment of the HO x species in the EB solver, which assumes that OH and HO 2 are at photochemical equilibrium at all times within the HO x family.This assumption is close to reality during the day, but becomes problematic at sunrise and sunset and is wrong at night, when the HO 2 lifetime can reach several hours at the surface.As a result, the OH mixing ratio calculated by the EB solver is overestimated by a factor of 10 compared to ASIS, which does not require any a priori assumption on chemical lifetimes and provides an accurate solution throughout sunset and nighttime.At sunrise and sunset ASIS reduces the chemical time step down to 10 s to solve the sharp transitions in the concentrations of short-lived species H, OH, O, and NO.Outside these critical (but short) periods, the Martian settings of RTOL and ATOL allow time steps that increase rapidly and may reach δt = 15 mn without sacrificing the accuracy.Thus, in the example of Fig. 9, at the surface level, the number of chemical time steps performed by ASIS over 1 Martian day is only 12 % larger than in the EB simulation.
The box-model simulations at 30 km are performed at the hygropause level where the production rate of HO x radicals by H 2 O photolysis is the largest.This results in a maximum stiffness of the system at sunrise and sunset, when the H 2 O photolysis rate varies rapidly.Those critical day-night transitions show large differences between the ASIS and the EB simulations.In the EB run, ozone is integrated implicitly by Eq. ( 3) at night and is assumed to be at photochemical equilibrium within the O x family during the day.This abrupt change in treatment contrasts with the smooth transition carried out with the time-step-adaptive scheme ASIS.At the price of a strong reduction of the time step to maintain the required accuracy, ASIS calculates an O 3 mixing ratio that is respectively 35 % larger and 20 % smaller than in the EB run at sunrise and sunset.Both solvers give the same results during the day.However, the more accurate description of the O 3 increase at sunset by ASIS induces a 5 % difference with the EB solver that persists into the night.Regarding OH, the simulation at 30 km confirms the weakness of the steadystate approximation for HO x at night in the EB scheme.In ASIS, the stiffness of the system diagnosed by the solver remains high in the first hours following sunset (due to strong curvature of the solution for H, not shown here) and leads to a reduction of the time step to about 30 s.The nighttime OH mixing ratio is larger by a factor of 2 to 4 than in the EB simulation.For this extreme case of stiffness in the Mars atmospheric chemistry, the total number of chemical time steps executed by ASIS over one Martian day is 65 % larger than in the EB simulation.
In its 3-D implementation, ASIS is called by the LMD GCM at each physical time step t = 15 mn.The ASIS settings in the GCM are identical to those of the box model presented earlier, i.e., the solver may select any sub-time-step value between t and the minimum value δt min = 10 s.To compare the GCM performances with ASIS and with the EB method, two simulations of 150 Martian solar days have been performed with each method starting with an identical initial situation.Figure 10 shows the number of sub-time-steps per physical time step of 15 mn in a GCM simulation of northern spring (Ls = 70 • ) using ASIS.For the three levels presented here (surface, 30 km, 80 km), the number of sub-time-steps is equal to 1 or 2 for a large fraction of time.This is the case when the chemical system is in equilibrium, far from the terminators at night or during the day.As in the MOCAGE model, at the terminators the number of sub-time-steps increases dramatically to cope with the change of chemical regime at the day-night transitions.The maximum number (40-50) is found at sunrise at 30 km and is essentially driven by the abrupt changes in OH and O 3 already seen in Fig. 10.At the surface, an increase in the number of sub-time-steps is also visible near the North Pole.This is related to fast heterogeneous reactions of HO x species on water-ice clouds (Lefèvre et al., 2008), a process similar to that occurring with chlorine on Earth stratospheric clouds.In those cases ASIS adopts a smaller time step to resolve with good accuracy a system that is locally away from chemical equilibrium.
Figure 11 compares at 30 km the results of GCM simulations using either the EB or the ASIS solver.Both schemes give distributions of O 3 and OH that are in general very close during daytime and away from the terminators.At the terminators, ASIS calculates O 3 amounts that are about 50 % larger than EB at sunrise and 25 % smaller at sunset.These large differences are similar to those found with the boxmodel runs (Fig. 9) but are limited in time and space.However, the better description of O 3 by ASIS across the terminators may be crucial when comparing the GCM to Martian ozone measurements performed at the terminators by the solar occultation technique.Regarding OH, the GCM results confirm the poor description of the HO x chemistry by the EB scheme at the terminators.At night, OH values calculated by EB are more than 30 % smaller than with ASIS.The amount of nighttime OH is small relative to daytime values.Thus, the bias in the EB scheme does not significantly affect the oxidizing capacity of Mars simulated by the GCM.Nevertheless, similarly to ozone, the more accurate description of OH and the Martian nighttime chemistry in general is an important advantage brought by ASIS for the interpretation of the numerous observations of nightglow or measurements by stellar occultation carried out on that planet.

Conclusions
The ASIS solver has been designed to cope with the various situations encountered within the numerical simulation of the atmospheric chemistry.The main properties of the solver are mass conservation, an approximation of the Jacobian matrix of the chemical fluxes that stabilizes the associated system of differential equations, a time-stepping varying module to control accuracy, and a code implementation that allows an easy adaptation to various chemical schemes.In box-model test cases, the numerical solutions obtained with the ASIS solver were found to be in good agreement with those of multistep algorithms like Rosenbrock's and Gear's methods.
The ASIS solver has been implemented in 3-D models of the Earth (MOCAGE) and Mars (LMD model) planets.The results with MOCAGE using ASIS reveals two weaknesses of the original semi-implicit solver.One is related to the calculation of the partitioning of the NO x species at the surface and the other to an overestimation of the ozone depletion in the Antarctic stratospheric vortex in Spring.In the simulation of the Mars atmosphere ASIS gives more accurate simulations during day-night transitions and at night for the HO x species.These results stress the importance of having accurate enough numerical solutions, otherwise differences between model simulations and observations could be wrongly attributed to missing chemistry or misrepresentation of some physical processes.
The model simulations show the benefit of using a chemical solver with good properties such as mass conservation and controlled accuracy.This objective can be achieved using multistep high-order algorithms, but the computational cost of those schemes increases rapidly with the number of species considered.Since ASIS is implicit and one-step, a single linear system has to be solved for each iteration.For this, direct or iterative algorithms can be used.The direct methods based on LU decomposition see their computational cost increasing at least quadratically with the number of species, whereas the cost of iterative solvers increases rather linearly.Within ASIS we found that the GMRES iterative algorithm is stable and efficient, and is competitive in terms of CPU cost compared to the direct DGESV algorithm.
In atmospheric models the computational cost is a key issue, and parallelization of the computations must be efficient to reduce the elapsed time spent for the simulations.As pointed out earlier, the amount of computation spent by ASIS to solve the chemical system can vary significantly from one grid point to another.This renders the work balancing of tasks more difficult if a domain decomposition strategy is adopted to implement the parallelization.As already discussed with the surface emissions, one possibility to diminish the number of iterations and the heterogeneity in the CPU used at each grid point is to account for nonchemical tendencies in the species continuity equations (term F of Eq. 4).Rather than updating the concentrations after each process, the resulting tendencies could be added and integrated within ASIS.This strategy has been adopted by, for example, Menut et al. (2013) for the CHIMERE model; it remains to be seen if the stability and the positivity of the solution can be maintained.
The present version of the ASIS solver addresses the evolution of the concentrations in gas phase only.For some applications the aqueous phase associated with the presence of clouds must be also considered (e.g., Leriche et al., 2013).The chemistry module has to solve both gaseous-and aqueous-phase chemistry as well as mass transfer reactions between gas and liquid phases.There is a priori no difficulty in adding the prognostic concentrations in the water phase to the system of equations and making a linearization similar to what is done in Eq. ( 6).However, the addition of aqueous reactions tends to increase the stiffness of the numerical ODE (Audiffren et al., 1998), so the performances of ASIS could diminish and may result in reduced time steps and increased computer time.
In conclusion, the ASIS solver can deal with many situations encountered in modeling atmospheric chemistry for a computational cost affordable by CTMs and GCMs that include comprehensive chemical schemes.Evolution of the ASIS solver to treat aqueous-phase chemistry is planned in the near future.
Code availability.The Fortran code to run the ASIS solver on the FLUX case described is Sect. 3 is available as a Supplement to the present article and can be downloaded from the CERFACS server (www.cerfacs.fr).
The code associated with the chemistry model includes subroutines that define the mechanism and those more specific to the ASIS solver.At this stage we have not developed an external driver or a pre-processor that would generate specific codes based on the adopted mechanism.This choice was done because our experience is that the maintenance of the driver outputs can be somewhat cumbersome when many developers work in parallel on a CTM.In addition, the code generated by the driver must be optimized often for the computer used and adapted to the CTM.It is therefore not used directly, which introduces further constraints on the maintenance of the overall code.
Once the definition of species and reactions is completed, the calculation of the matrices (Eq.7) is done by the fill_matrix routine, the time steps are monitored by the define_dt routine, and the solution of the linear systems by the Solvesys routine.
The Supplement related to this article is available online at doi:10.5194/gmd-10-1467-2017-supplement.

D
. Cariolle et al.: ASIS: an adaptive solver ∂C/∂t = f (t, C) = P(t, C) − L(t, C).C, with σ = −1 if m = k and l = k, σ = −2 if l = m = k,and σ = 1 if l and m = k , where C k is the concentration of the k species.

Figure 1 .
Figure 1.FLUX case.Time evolution of selected species (O 3 , NO 2 , NO 3 , OH) for the R1, G1, A1, A2, and A3 experiments.The left column shows the volume mixing ratios and the right one the differences in percentage (%) relative to the R1 experiment.The color code is the following: blue for G1, orange for A1, red for A2, and purple for A3.

Figure 3 .
Figure 3. STRATO case.Number of time steps of the ASIS solver for each interval of 15 min in experiment AS2.

Figure 4 .
Figure 4. STRATO case.Time evolution of selected species (CH 2 O, NO 2 , NO 3 , OH) for the RS1, and AS2 experiments.(a) shows the mixing ratios or concentrations and (b) the differences in percentage (%) for AS2 relative to the RS1 experiment.

Figure 5 .
Figure 5. Number of sub-time-steps per time step of 15 min in the MA simulation for the 15 September at midday.Three levels are presented representative of the stratosphere (50 hPa), the midtroposphere (540 hPa), and the surface.

Figure 6 .
Figure 6.Zonal mean distributions of O 3 , CO, OH, and HNO 3 for the month of September.The panels in (a) show results for the reference simulation, MR; those in (b) show results of the MA simulation with the use of the ASIS solver.

Figure 7 .
Figure 7. Monthly mean distributions of NO 2 , HNO 3 , and N 2 O 5 at the surface for the month of October after a 2 month integration.The panels in (a) show results of the reference simulation, MR; those in (b) show results of the MA simulation with the use of the ASIS solver.The MR simulation underevaluates the N 2 O 5 in the lower troposphere.

Figure 8 .
Figure 8. Evolution of the total ozone column over the Dumont d'Urville and Dome C Antarctic stations.The dots are the observations of the SAOZ instrument, the orange line is the evolution calculated in the reference simulation, MR, and the red line the same output from the simulation MA using the ASIS solver.

Figure 9 .
Figure 9.Comparison of the Euler-backward (EB) and ASIS solvers applied to the Mars box-model version.The left column shows the mixing ratios of O 3 and OH and the right one the ratio between the ASIS and EB experiments.Results are presented at 30 km (a) and at the surface (b), for equatorial conditions in northern spring (solar longitude Ls = 70 • ).Local noon is at day 0.5.

Figure 10 .
Figure 10.Number of sub-time-steps per time interval of 15 min in the LMD Mars GCM in northern spring (instantaneous result at solar longitude Ls = 70 • and day 150 of the simulation).Three altitude levels are represented at 80 km (a), 30 km (b), and the surface (c).Local noon is located at longitude zero.The white contour represents topography, with a 4 km interval.

Figure 11 .
Figure 11.Distribution of O 3 (left) and OH (right) at 30 km calculated by the LMD Mars GCM in northern spring (instantaneous result at solar longitude Ls = 70 • and day 150 of the simulation).Top: Euler-backward (EB) solver.Middle: ASIS solver.Bottom: relative difference (%) between ASIS and EB, using thresholds of 10 ppbv for O 3 and 10 pptv for OH.Local noon is located at longitude zero.The white contour represents topography, with a 4 km interval.Off-scale differences in O 3 and OH are limited to ±50 % and ±30 % respectively.
tions.The reactions are classified in three groups: 1/ A → b B + c C 2/ A + A → b B + c C 3/ a A + b B → c C + d D

Table 2 .
VOC emissions in the FLUX test case.

Table 3 .
Initial conditions for the FLUX and STRATO test cases.

Table 4 .
List of the different settings used by the 0-D model for the FLUX test case.The mean time step is denoted by δt m and CPU-Ratio is the CPU time used in each simulation relative to R1.The simulations are performed in the MATLAB environment.