/home66/gary/public_html/cloudy/c08_branch/source/ion_solver.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ion_solver solve the bi-diagonal matrix for ionization balance */
#include "cddefines.h"
#include "yield.h"
#include "prt.h"
#include "continuum.h"
#include "iso.h"
#include "dynamics.h"
#include "grainvar.h"
#include "hmi.h"
#include "mole.h"
#include "thermal.h"
#include "thirdparty.h"
#include "conv.h"
#include "secondaries.h"
#include "phycon.h"
#include "atmdat.h"
#include "heavy.h"
#include "elementnames.h"
#include "dense.h"
#include "radius.h"
#include "ionbal.h"

void solveions(double *ion, double *rec, double *snk, double *src,
                           long int nlev, long int nmax);

void ion_solver(
        /* this is element on the c scale, H is 0 */
        long int nelem, 
        /* option to print this element when called */
        bool lgPrintIt)
{
        /* use tridiag or general matrix solver? */
        bool lgTriDiag=false;
        long int ion, 
          limit, 
          IonProduced, 
          nej, 
          ns,
          jmax=-1;
        double totsrc;
        double dennew, rateone;
        bool lgNegPop;
        static bool lgMustMalloc=true;
        static double *sink, *source , *sourcesave;
        int32 nerror;
        static int32 *ipiv;
        long ion_low, ion_range, i, j, ion_to , ion_from;
        static double sum_dense;
        /* only used for debug printout */
        static double *auger;

        double abund_total, renormnew;
        bool lgHomogeneous = true;
        static double *xmat , *xmatsave;

        DEBUG_ENTRY( "ion_solver()" );

        /* this is on the c scale, so H is 0 */
        ASSERT( nelem >= 0);
        ASSERT( dense.IonLow[nelem] >= 0 );
        ASSERT( dense.IonHigh[nelem] >= 0 );

        /* H is special because its abundance spills into three routines -
         * the ion/atom solver (this routine), the H-mole solvers (hmole), and
         * the heavy molecule solver.  xmolecules only includes the heavy mole
         * part for H only.  So the difference between gas_phase and xmolecules
         * includes the H2 part of the molecular network.  This branch does
         * this special H case, then the general case (heavy elements are
         * not part of the H2 molecular network) */

        /* >>chng 01 dec 07, define abund_total, total atom and ion abundance 
         * here, removing molecules */
        if( nelem == ipHYDROGEN )
        {
                /* Hydrogen is a special case since hmole does not include the
                 * atom/molecules - hmole collapses its network into H = H0 + H+
                 * and then forces the solution determined here for partitioning
                 * between these two */
                abund_total = dense.xIonDense[nelem][0] + dense.xIonDense[nelem][1];
        }
        else
        {
                abund_total = SDIV( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
        }

        /* >>chng 04 nov 22,
         * if gas nearly all atomic/ionic do not let source - sink terms from 
         * molecular network force system of balance equations to become 
         * inhomogeneous
         * what constitutes a source or sink IS DIFFERENT for hydrogen and the rest 
         * the H solution must couple with hmole - and its defn of source and sink.  For instance, oxygen charge
         * transfer goes into source and sink terms for hydrogen.  So we never hose source and sink for H.
         * for the heavy elements, which couple onto comole, mole.source and sink represent terms that
         * remove atoms and ions from the ionization ladder.  their presence makes the system of
         * equations inhomogeneous.  we don't want to do this when comole has a trivial effect on
         * the ionization balance since the matrix becomes unstable */
        /* >>chng 04 dec 06, limit from 0.01 to 1e-10 as per NA suggestion */
        if( nelem>ipHYDROGEN && dense.xMolecules[nelem]/SDIV(dense.gas_phase[nelem]) < 1e-10 )
        {
                for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
                {
                        mole.source[nelem][ion] = 0.;
                        mole.sink[nelem][ion] = 0.;
                }
        }

        /* protect against case where all gas phase abundances are in molecules, use the
         * atomic and first ion density from the molecule solver 
         * >>chng 04 aug 15, NA change from 10 0000 to 10 pre-coef on
         * FLT_EPSILON for stability in PDR 
         * the factor 10.*FLT_EPSILON also appears in mole_h_step in total H 
         * conservation */
        if( fabs( dense.gas_phase[nelem] -  dense.xMolecules[nelem])/SDIV(dense.gas_phase[nelem] ) <
                10.*FLT_EPSILON )
        {
                double renorm;
                /* >>chng 04 jul 31, add logic to conserve nuclei in fully molecular limit;
                 * in first calls, when searching for solution, we may be VERY far 
                 * off, and sum of first ion and atom density may be far larger 
                 * than difference between total gas and molecular densities,
                 * since they reflect the previous evaluation of the solution.  Do 
                 * renorm to cover this case */
                /* first form sum of all atoms and ions */
                realnum sum = 0.;
                for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
                        sum += dense.xIonDense[nelem][ion];
                /* now renorm to this sum - this should be unity, and is not if we have
                 * now conserved particles, due to molecular fraction changing */
                renorm = dense.gas_phase[nelem] / SDIV(sum + dense.xMolecules[nelem] );

                abund_total = renorm * sum;
        }

        /* negative residual density */
        if( abund_total < 0. )
        {
                /* >>chng 05 dec 16, do not abort when net atomic/ ionic abundance 
                 * is negative due to chem net having too much of a species - this 
                 * happens naturally when ices become well formed, but the code will 
                 * converge away from it.  Rather, we set the ionization
                 * convergence flag to "not converge" and soldier on 
                 * if negative populations do not go away, we will eventually 
                 * terminate due to convergence failures */
                /* print error if not search */
                if(!conv.lgSearch )
                        fprintf(ioQQQ,
                                "PROBLEM ion_solver - neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
                                nelem,
                                fabs(abund_total) / SDIV(dense.xMolecules[nelem]),
                                conv.nTotalIoniz );
                /* fix up is to use half the positive abundance, assuming chem is 
                 * trying to get too much of this species */
                abund_total = -abund_total/2.;

                /* say that ionization is not converged - do not abort - but if 
                 * cannot converge away from negative solution, this will become a 
                 * convergence failure abort */
                conv.lgConvIoniz = false;
                strcpy( conv.chConvIoniz, "neg ion" );
        }

        /* return if IonHigh is zero, since no ionization at all */
        if( dense.IonHigh[nelem] == 0 )
        {
                /* set the atom to the total gas phase abundance */
                dense.xIonDense[nelem][0] = (realnum)abund_total;
                return;
        }

        /* >>chng 01 may 09, add option to force ionization distribution with 
         * element name ioniz */
        if( dense.lgSetIoniz[nelem] )
        {
                for( ion=0; ion<nelem+2; ++ion )
                {
                        dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*
                                (realnum)abund_total;
                }
                return;
        }

        /* impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1
         * HIonFrac(nelem) is stripped to hydrogen */
        /* >>chng 01 oct 30, to assert */
        ASSERT( (dense.IonHigh[nelem] < nelem + 1) || 
                iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. );

        /* zero out the ionization and recombination rates that we will modify here,
         * but not the iso-electronic stages which are done elsewhere,
         * the nelem stage of ionization is he-like,
         * the nelem+1 stage of ionization is h-like */

        /* loop over stages of ionization that we solve for here, 
         * up through and including one less than nelem-NISO,
         * never actually do highest NISO stages of ionization since they
         * come from the ionization ratio from the next lower stage */
        limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);

        /* the full range of ionization - this is number of ionization stages */
        ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
        ASSERT( ion_range <= nelem+2 );

        ion_low = dense.IonLow[nelem];

        /* PARALLEL_MODE true if there can be several instances of this routine 
         * running in the same memory pool - we must have fresh instances of the 
         * arrays for each */
        if( PARALLEL_MODE || lgMustMalloc )
        {
                long int ncell=LIMELM+1;
                lgMustMalloc = false;
                if( PARALLEL_MODE )
                        ncell = ion_range;

                /* this will be "new" matrix, with non-adjacent coupling included */
                xmat=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
                xmatsave=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ));
                sink=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
                source=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
                sourcesave=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ));
                ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)ncell ));
                /* auger is used only for debug printout - it is special because with multi-electron
                 * Auger ejection, very high stages of ionization can be produced, well beyond
                 * the highest stage considered here.  so we malloc to the limit */
                auger=(double*)MALLOC( (sizeof(double)*(unsigned)(LIMELM+1) ));
        }

        for( i=0; i<ion_range;i++ )
        {
                source[i] = 0.;
        }

        /* this will be used to address the 2d arrays */
#       undef MAT
#       define MAT(M_,I_,J_)    (*((M_)+(I_)*(ion_range)+(J_)))

        /* zero-out loop comes before main loop since there are off-diagonal
         * elements in the main ionization loop, due to multi-electron processes,
         * TotIonizRate and TotRecom were already set in h-like and he-like solvers 
         * other recombination rates were already set by routines responsible for them */
        for( ion=0; ion <= limit; ion++ )
        {
                ionbal.RateIonizTot[nelem][ion] = 0.;
        }

        /* auger is used only for debug printout - it is special because with multi-electron
         * Auger ejection, very high stages of ionization can be produced, well beyond
         * the highest stage considered here.  so we malloc to the limit */
        for( i=0; i< LIMELM+1; ++i )
        {
                auger[i] = 0.;
        }

        /* zero out xmat */
        for( i=0; i< ion_range; ++i )
        {
                for( j=0; j< ion_range; ++j )
                {
                        MAT( xmat, i, j ) = 0.;
                }
        }

        {
                /* this sets up a fake ionization balance problem, with a trivial solution,
                 * for debugging the ionization solver */
                enum {DEBUG_LOC=false};
                if( DEBUG_LOC && nelem==ipCARBON )
                {
                        broken();/* within optional debug print statement */
                        dense.IonLow[nelem] = 0;
                        dense.IonHigh[nelem] = 3;
                        abund_total = 1.;
                        ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
                        /* make up ionization and recombination rates */
                        for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
                        {
                                double fac=1;
                                if(ion)
                                        fac = 1e-10;
                                ionbal.RateRecomTot[nelem][ion] = 100.;
                                for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
                                {
                                        /* direct photoionization of this shell */
                                        ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
                                }
                        }
                }
        }

        /* Now put in all recombination and ionization terms from CO_mole() that 
         * come from molecular reactions. this traces molecular process that 
         * change ionization stages with this ladder - but do not remove from 
         * the ladder */
        for( ion_to=dense.IonLow[nelem]; ion_to <= limit; ion_to++ )
        {
                for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
                {
                        /* do not do ion onto itself */
                        if( ion_to != ion_from )
                        {
                                /* this is the rate coefficient for charge transfer from ion to ion_to */
                                /*rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];*/
                                rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to];
                                ASSERT( rateone >= 0. );
                                MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
                                MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
                        }
                }
        }

        /* now get actual arrays of ionization and recombination processes,
         * but only for the ions that are done as two-level systems */
        /* in two-stage system, atom + first ion, limit is zero but must
         * include gv.GrainChTrRate[nelem][1][0] */
        /* grain charge transfer */
        if( gv.lgDustOn && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn )
        {
                long int low;
                /* do not double count this process for atoms that are in the co network - we use
                 * a net recombination coefficient derived from the co solution, this includes grain ct */
                /* >>chng 05 dec 23, add mole.lgElem_in_chemistry */
                if( mole.lgElem_in_chemistry[nelem] )
                /*if( nelem==ipHYDROGEN ||nelem==ipCARBON ||nelem== ipOXYGEN ||nelem==ipSILICON ||nelem==ipSULPHUR 
                         ||nelem==ipNITROGEN ||nelem==ipCHLORINE )*/
                {
                         low = MAX2(1, dense.IonLow[nelem] );
                }
                else
                        low = dense.IonLow[nelem];

                for( ion_to=low; ion_to <= dense.IonHigh[nelem]; ion_to++ )
                {
                        for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
                        {
                                /* do not do ion onto itself */
                                if( ion_to != ion_from )
                                {
                                        /* this is the rate coefficient for charge transfer from ion to ion_to
                                         * both grain charge transfer ionization and recombination */
                                        rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];
                                        MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
                                        MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
                                }
                        }
                }
        }

        for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
        {
                /* thermal & secondary collisional ionization */
                rateone = ionbal.CollIonRate_Ground[nelem][ion][0] +
                        secondaries.csupra[nelem][ion] +
                        /* inner shell ionization by UTA lines */
                        ionbal.UTA_ionize_rate[nelem][ion];
                ionbal.RateIonizTot[nelem][ion] += rateone;

                /* UTA ionization */
                if( ion+1-ion_low < ion_range )
                {
                        /* depopulation processes enter with negative sign */
                        MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
                        MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
                }

                /* total recombination rate */
                if( ion-1-ion_low >= 0 )
                {
                        /* loss of this ion due to recombination to next lower ion stage */
                        MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
                        MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
                }

                /* loop over all atomic sub-shells to include photoionization */
                for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
                {
                        /* direct photoionization of this shell */
                        ionbal.RateIonizTot[nelem][ion] += ionbal.PhotoRate_Shell[nelem][ion][ns][0];

                        /* this is the primary ionization rate - add to diagonal element,
                         * test on ion stage is so that we don't include ionization from the very highest
                         * ionization stage to even higher - since those even higher stages are not considered
                         * this would appear as a sink - but populations of this highest level is ensured to
                         * be nearly trivial and neglecting it production of even higher ionization OK */
                        /* >>chng 04 nov 29 RJRW, include following in this branch so only
                         * evaluated when below ions done with iso-sequence */
                        if( ion+1-ion_low < ion_range )
                        {
                                /* this will be redistributed into charge states in following loop */ 
                                MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.PhotoRate_Shell[nelem][ion][ns][0];

                                /* t_yield::Inst().nelec_eject(nelem,ion,ns) is total number of electrons that can
                                 * possibly be freed 
                                 * loop over nej, the number of electrons ejected including the primary,
                                 * nej = 1 is primary, nej > 1 includes primary plus Auger 
                                 * t_yield::Inst().elec_eject_frac is prob of nej electrons */
                                for( nej=1; nej <= t_yield::Inst().nelec_eject(nelem,ion,ns); nej++ )
                                {
                                        /* this is the ion that is produced by this ejection,
                                         * limited by highest possible stage of ionization -
                                         * do not want to ignore ionization that go beyond this */
                                        IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
                                        rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
                                                t_yield::Inst().elec_eject_frac(nelem,ion,ns,nej-1);
                                        /* >>chng 04 sep 06, above had included factor of nej to get rate, but
                                         * actually want events into particular ion */
                                        /* number of electrons ejected
                                         *(double)nej; */
                                        /* it goes into this charge state - recall upper cap due to ion stage trimming 
                                         * note that compensating loss term on diagonal was done before this
                                         * loop, since frac_elec_eject adds to unity */
                                        MAT( xmat, ion-ion_low, IonProduced-ion_low ) += rateone;

                                        /* only used for possible printout - multiple electron Auger rate  -
                                         * do not count one-electron as Auger */
                                        if( nej>1 )
                                                auger[IonProduced-1] += rateone;
                                }
                        }
                }

                /* this is charge transfer ionization of this species by hydrogen and helium */
                rateone = 
                        atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 
                        atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];
                ionbal.RateIonizTot[nelem][ion] += rateone;
                if( ion+1-ion_low < ion_range )
                {
                        MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
                        MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
                }
        }

        /* after this loop, ionbal.RateIonizTot and ionbal.RateRecomTot have been defined for the
         * stages of ionization that are done with simple */
        /* begin loop at first species that is treated with full model atom */
        j = MAX2(0,limit+1);
        /* possible that lowest stage of ionization is higher than this */
        j = MAX2( ion_low , j );
        for( ion=j; ion<=dense.IonHigh[nelem]; ion++ )
        {
                ASSERT( ion>=0 && ion<nelem+2 );
                /* use total ionization/recombination rates for species done with ISO solver */
                if( ion+1-ion_low < ion_range )
                {
                        /* depopulation processes enter with negative sign */
                        MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.RateIonizTot[nelem][ion];
                        MAT( xmat, ion-ion_low, ion+1-ion_low ) += ionbal.RateIonizTot[nelem][ion];
                }

                if( ion-1-ion_low >= 0 )
                {
                        /* loss of this ion due to recombination to next lower ion stage */
                        MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
                        MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
                }
        }

        /* this will be sum of source and sink terms, will be used to decide if
         * matrix is singular */
        totsrc = 0.;
        /*>>chng 06 nov 28 only include source from molecules if we have an estimated first 
         * solution - first test is that we have called mole at least twice,
         * second test is that we are on a later iteration */
        if( conv.nTotalIoniz > 1 || iteration > 1 )
        {

                for( i=0; i<ion_range;i++ )
                {
                        ion = i+ion_low;

                        /* these are the external source and sink terms */
                        /* source first */
                        /* need negative sign to get positive pops */
                        source[i] -= mole.source[nelem][ion];

                        totsrc += mole.source[nelem][ion];
                        /* sink next */
                        /*MAT( xmat, i, i ) += mole.sink[nelem][ion];*/
                        MAT( xmat, i, i ) -= mole.sink[nelem][ion];
                }
        }

        /* matrix is not homogeneous if source is non-zero */
        if( totsrc != 0. )
                lgHomogeneous = false;

        /* chng 03 jan 13 rjrw, add in dynamics if required here,
         * last test - only do advection if we have not overrun the radius scale */
        if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection 
                /*&& radius.depth < dynamics.oldFullDepth*/ )
        {
                for( i=0; i<ion_range;i++ )
                {                        
                        ion = i+ion_low;
                        /*MAT( xmat, i, i ) += dynamics.Rate;*/
                        MAT( xmat, i, i ) -= dynamics.Rate;
                        /*src[i] += dynamics.Source[nelem][ion];*/
                        source[i] -= dynamics.Source[nelem][ion];
                        /* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(xmat,i,i),
                          dynamics.Rate, dynamics.Source[nelem][ion]);*/
                }
                lgHomogeneous = false;
        }

        for( i=0; i< ion_range; ++i )
        {
                for( j=0; j< ion_range; ++j )
                {
                        MAT( xmatsave, i, j ) = MAT( xmat, i, j );
                }
                sourcesave[i] = source[i];
        }

        /* >>chng 06 nov 21, for very high ionization parameter sims the H molecular
        * fraction can become so small that atom = atom + molecule.  In this case we
        * will not count system as an inhomogeneous case since linear algebra package
        * will fail.  If molecules are very small, count as homogeneous.
        * Note that we have already added sink terms to the main matrix and they
        * will not be removed, a possible source of error, but they must not
        * have been significant, given that the molecular fraction is so small */
        if( !lgHomogeneous && ion_range==2 )
        {
                /* solve algebraically */
                double a = MAT( xmatsave, 0, 0 ), 
                        b = MAT( xmatsave, 1, 0 ) ,
                        c = MAT( xmatsave, 0, 1 ) ,
                        d = MAT( xmatsave, 1, 1 );
                b = SDIV(b);
                d = SDIV(d);
                double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
                if( fabs(ratio1-ratio2)/MAX2(fratio1,fratio2) <DBL_EPSILON*1e4 )
                {
                        abund_total = SDIV( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
                        lgHomogeneous = true;
                }
        }
#       if 0
        if( nelem==ipHYDROGEN && 
                fabs(dense.xMolecules[nelem]) / SDIV(dense.gas_phase[ipHYDROGEN]) <DBL_EPSILON*100. )
        {
                abund_total = SDIV( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
                lgHomogeneous = true;
        }
#       endif

        /* this is true if no source terms 
         * we will use total population and species conservation to replace one
         * set of balance equations since overdetermined */
        if( lgHomogeneous )
        {
                double rate_ioniz=1., rate_recomb /*, scale = 0.*/;
                /* Simple estimate of most abundant ion */
                jmax = 0;
                for( i=0; i<ion_range-1;i++)
                { 
                        ion = i+ion_low;
                        rate_ioniz *= ionbal.RateIonizTot[nelem][ion];
                        rate_recomb = ionbal.RateRecomTot[nelem][ion];
                        /* trips if ion rate zero, so ll the gas will be more neutral than this */
                        if( rate_ioniz == 0)
                                break;
                        /* rec rate is zero */
                        if( rate_recomb <= 0.) 
                                break;

                        rate_ioniz /= rate_recomb;
                        if( rate_ioniz > 1.) 
                        {
                                /* this is peak ionization stage */
                                jmax = i;
                                rate_ioniz = 1.;
                        }
                }

                /* replace its matrix elements with population sum */
                for( i=0; i<ion_range;i++ )
                {
                        MAT(xmat,i,jmax) = 1.;
                }
                source[jmax] = abund_total;
        }

        for( i=0; i< ion_range; ++i )
        {
                for( j=0; j< ion_range; ++j )
                {
                        MAT( xmatsave, i, j ) = MAT( xmat, i, j );
                }
                sourcesave[i] = source[i];
        }

        if( false && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 ) 
        {
                fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
                fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot[nelem][0], ionbal.RateIonizTot[nelem][1]);
                fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
                fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
        }

        {
                /* option to print matrix */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && nzone==1 && lgPrintIt )
                {
                        fprintf( ioQQQ, 
                                " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
                                nelem , ion_range,limit  , conv.nTotalIoniz );
                        if( lgHomogeneous )
                                fprintf(ioQQQ , "Homogeneous \n");
                        for( i=0; i<ion_range; ++i )
                        {
                                for( j=0;j<ion_range;j++ )
                                {
                                        fprintf(ioQQQ,"%e\t",MAT(xmatsave,i,j));
                                }
                                fprintf(ioQQQ,"\n");
                        }
                        fprintf(ioQQQ,"source follows\n");
                        for( i=0; i<ion_range;i++ )
                        {
                                fprintf(ioQQQ,"%e\t",source[i]);
                        }
                        fprintf(ioQQQ,"\n");
                }
        }

        /* first branch is tridiag, following branch is just general matrix solver */
        if( lgTriDiag ) 
        {
                bool lgLapack=false , lgTriOptimized=true;
                /* there are two choices for the tridiag solver */
                if(lgLapack) 
                {
                        /* this branch uses lapack tridiag solver, and should work 
                         * it is hardwired to not be used because not extensively tested
                         * issues - is this slower than others, and is it robust in
                         * singular cases? */
                        bool lgBad = false;
                        /* Use Lapack tridiagonal solver */
                        double *dl, *d, *du;

                        d = (double *) MALLOC((unsigned)ion_range*sizeof(double));
                        du = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
                        dl = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));

                        for(i=0;i<ion_range-1;i++) 
                        {
                                du[i] = MAT(xmat,i+1,i);
                                d[i] = MAT(xmat,i,i);
                                dl[i] = MAT(xmat,i,i+1);
                        }
                        d[i] = MAT(xmat,i,i);

#                       if 0
                        int i1, i2;
                        i1 = ion_range;
                        i2 = 1;
                        lgBad = dgtsv_wrapper(&i1, &i2, dl, d, du, source, &i2);
#                       endif
                        if( lgBad )
                                fprintf(ioQQQ," dgtsz error.\n");

                        free(dl);free(du);free(d);
                } 
                else if(lgTriOptimized)
                {
                        /* Use tridiagonal solver re-coded to avoid rounding errors
                         * on diagonal -- uses determination of jmax for the
                         * singular case, but is otherwise independent of the array
                         * filling code above */

                        for(i=0;i<ion_range;i++) 
                        {
                                source[i] = sink[i] = 0.;
                        }
                        if(nelem == ipHYDROGEN && !hmi.lgNoH2Mole)
                        {
                                for(i=0;i<ion_range;i++) 
                                {
                                        ion = i+ion_low;
                                        source[i] += mole.source[nelem][ion];
                                        sink[i] += mole.sink[nelem][ion];
                                }
                        }
                        if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth )
                        {
                                for(i=0;i<ion_range;i++) 
                                {
                                        ion = i+ion_low;
                                        source[i] += dynamics.Source[nelem][ion];
                                        sink[i] += dynamics.Rate;
                                }
                        }

                        solveions(ionbal.RateIonizTot[nelem]+ion_low,ionbal.RateRecomTot[nelem]+ion_low,
                                                                sink,source,ion_range,jmax);
                }
        } 
        else
        {
                /* this is the default solver - now get new solution */
                nerror = 0;
                /* Use general matrix solver */
                getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror);
                if( nerror != 0 )
                {
                        fprintf( ioQQQ, 
                                " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, limit=%li, nConv %li xmat follows\n",
                                nelem , 
                                elementnames.chElementSym[nelem],
                                ion_range,
                                limit , 
                                conv.nTotalIoniz );
                        for( i=0; i<ion_range; ++i )
                        {
                                for( j=0;j<ion_range;j++ )
                                {
                                        fprintf(ioQQQ,"%e\t",MAT(xmatsave,j,i));
                                }
                                fprintf(ioQQQ,"\n");
                        }
                        fprintf(ioQQQ,"source follows\n");
                        for( i=0; i<ion_range;i++ )
                        {
                                fprintf(ioQQQ,"%e\t",source[i]);
                        }
                        fprintf(ioQQQ,"\n");
                        cdEXIT(EXIT_FAILURE);
                }
                getrs_wrapper('N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror);
                if( nerror != 0 )
                {
                        fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
                                nelem , ion_range );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        {
                /*@-redef@*/
                /* this is to debug following failed assert */
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && nelem == ipHYDROGEN )
                {
                        fprintf(ioQQQ,"debuggg ion_solver1 \t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n", 
                                fnzone,
                                phycon.te,
                                dense.eden,
                                ionbal.RateIonizTot[nelem][0] , 
                                ionbal.RateRecomTot[nelem][0]);
                        fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[ipHYDROGEN][0], mole.source[ipHYDROGEN][1]);
                        fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[ipHYDROGEN][0], mole.sink[ipHYDROGEN][1]);
                        fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0],
                                ionbal.RateIonizTot[nelem][0]/ionbal.RateRecomTot[nelem][0]);
                }
        }

        for( i=0; i<ion_range; i++ )
        {
                /* is this blows, source[i] is NaN */
                ASSERT( !isnan( source[i] ) );
        }

#       define RJRW 0
        if( RJRW && 0 )
        {
                /* verify that the rates are sensible */
                double test;
                for(i=0; i<ion_range; i++) {
                        test = 0.;
                        for(j=0; j<ion_range; j++) {
                                test = test+source[j]*MAT(xmatsave,j,i);
                        }
                        fprintf(ioQQQ,"%e\t",test);
                }
                fprintf(ioQQQ,"\n");

                test = 0.;
                fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total);
                for( ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
                {
                        if( ionbal.RateRecomTot[nelem][ion] != 0 && source[ion-ion_low] != 0 )
                                fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low],
                                                                source[ion-ion_low+1]/source[ion-ion_low],
                                                                ionbal.RateIonizTot[nelem][ion]/ionbal.RateRecomTot[nelem][ion]);
                        else
                                fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]);
                        test += source[ion-ion_low];
                }
                test += source[ion-ion_low];
        }

        /* 
         * >> chng 03 jan 15 rjrw:- terms are now included for
         * molecular sources and sinks of H and H+.
         *
         * When the network is not in equilibrium, this will lead to a
         * change in the derived abundance of H and H+ when after the
         * matrix solution -- the difference between `renorm' and 1. is a
         * measure of the quality of the solution (it will be 1. if the
         * rate of transfer into Ho/H+ balances the rate of transfer
         * out, for the consistent relative abundances).
         *
         * We therefore renormalize to keep the total H abundance
         * correct -- only the molecular network is allowed to change
         * this.
         *
         * To do this, only the ion abundances are corrected, as the
         * molecular abundances may depend on several different
         * conserved species.
         *
         */
        if( lgHomogeneous )
        {
                dense.xIonDense[nelem][ion_low] = (realnum)abund_total;
                for( i=1;i < ion_range; i++ )
                {
                        dense.xIonDense[nelem][i+ion_low] = 0.;
                }
        }

        renormnew = 1.;
        /* >>chng 06 mar 17, comment out test on old full depth - keep old solution if overrun scale */
        if(iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection /*&& radius.depth < dynamics.oldFullDepth */ && 
                        nelem == ipHYDROGEN && hmi.lgNoH2Mole)
        {
                /* The normalization out of the matrix solution is correct and
                 * should be retained if: dynamics is on and the total
                 * abundance of HI & H+ isn't being controlled by the
                 * molecular network */
                renormnew = 1.;
        }
        else
        {
                dennew = 0.;
                sum_dense = 0.;

                /* find total population to renorm - also here check that negative pops did not occur */
                for( i=0;i < ion_range; i++ )
                {
                        ion = i+ion_low;
                        sum_dense += dense.xIonDense[nelem][ion];
                        dennew += source[i];
                } 

                if( dennew > 0.)
                {
                        renormnew = sum_dense / dennew;
                }
                else
                {
                        renormnew = 1.;
                }
        }
        /* check not negative, should be +ve, can be zero if species has become totally molecular.
         * this happens for hydrogen if no cosmic rays, or cr ion set very low */
        if( renormnew < 0)
        {
                fprintf(ioQQQ,"PROBLEM impossible value of renorm \n");
        }
        ASSERT( renormnew>=0 );

        /* source[i] contains new solution for ionization populations
         * save resulting abundances into main ionization density array, 
         * while checking whether any negative level populations occurred */
        lgNegPop = false;
        for( i=0; i < ion_range; i++ )
        {
                ion = i+ion_low;

                /*fprintf(ioQQQ," %li %li %.3e %.3e\n",nelem,ion,src[ion-ion_low+1],src[ion-ion_low]);
                        pop_ion_ov_neut[ion] = src[ion-ion_low+1]/src[ion-ion_low]; */
                if( source[i] < 0. )
                {
                        /* >>chng 04 dec 04, put test on neg abund here, don't print uncles value is very -ve */
                        /* >>chng 06 feb 28, from -1e-10 to -1e-9, sim func_t10 had several negative
                         * during initial search, due to extremely high ionization */
                        /* >>chng 06 mar 11, from 1e-9 to 2e-9 make many struc elements floats from double */
                        if( source[i]<-2e-9 )
                                fprintf(ioQQQ,
                                " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
                                nelem , 
                                elementnames.chElementSym[nelem],
                                ion , 
                                source[i] , 
                                TorF(conv.lgSearch) ,
                                nzone ,
                                iteration );
                        source[i] = 0.;
                        /* if this is one of the iso seq model atoms then must also zero out Lya and Pop2Ion */
                        if( ion == nelem+1-NISO )
                        {
                                long int ipISO = nelem - ion;
                                ASSERT( ipISO>=0 && ipISO<NISO );
                                StatesElem[ipISO][nelem][0].Pop = 0.;
                        }
                }

                /* use new solution */
                dense.xIonDense[nelem][ion] = (realnum)(source[i]*renormnew);
        }

        /* Zero levels with abundances < 1e-25 which which will suffer numerical noise */
        while( dense.IonHigh[nelem] > dense.IonLow[nelem] && 
                dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total )
        {
                ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
                /* zero out abundance and heating due to stage of ionization we are about to zero out */
                dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
                thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
                /* decrement counter */
                --dense.IonHigh[nelem];
        }

        /* sanity check, either offset stages of low and high ionization,
         * or no ionization at all */
        ASSERT( (dense.IonLow[nelem] < dense.IonHigh[nelem]) ||
                (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) );

        /* this should not happen */
        if( lgNegPop )
        {
                fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n", 
                  elementnames.chElementNameShort[nelem], nzone );

                fprintf( ioQQQ, " Populations were" );
                for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
                {
                        fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
                }
                fprintf( ioQQQ, "\n" );

                fprintf( ioQQQ, " destroy vector =" );
                for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot[nelem][ion-1] );
                }
                fprintf( ioQQQ, "\n" );

                /* print some extra stuff if destroy was negative */
                if( lgNegPop )
                {
                        fprintf( ioQQQ, " CTHeavy  vector =" );
                        for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                        {
                                fprintf( ioQQQ, "%9.1e", atmdat.HeCharExcIonOf[nelem][ion] );
                        }
                        fprintf( ioQQQ, "\n" );

                        fprintf( ioQQQ, " HCharExcIonOf vtr=" );
                        for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                        {
                                fprintf( ioQQQ, "%9.1e", atmdat.HCharExcIonOf[nelem][ion] );
                        }
                        fprintf( ioQQQ, "\n" );

                        fprintf( ioQQQ, " CollidRate  vtr=" );
                        for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                        {
                                fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
                        }
                        fprintf( ioQQQ, "\n" );

                        /* photo rates per subshell */
                        fprintf( ioQQQ, " photo rates per subshell, ion\n" );
                        for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                        {
                                fprintf( ioQQQ, "%3ld", ion );
                                for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
                                {
                                        fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
                                }
                                fprintf( ioQQQ, "\n" );
                        }
                }

                /* now check out creation vector */
                fprintf( ioQQQ, " create  vector =" );
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                ContNegative();
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }

        /* option to print ionization and recombination arrays
         * prt flag set with print array print arrays command */
        if( prt.lgPrtArry[nelem] || lgPrintIt )
        {
                /* say who we are, what we are doing .... */
                fprintf( ioQQQ, 
                        "\n %s ion_solver DEBUG ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n", 
                        elementnames.chElementSym[nelem],
                        elementnames.chElementName[nelem],
                        fnzone,
                        phycon.te , 
                        dense.eden,
                        dense.gas_phase[nelem],
                        abund_total ,
                        dense.xMolecules[nelem] );
                /* total ionization rate, all processes */
                fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", ionbal.RateIonizTot[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* sinks from the chemistry network */
                fprintf(ioQQQ," %s mole sink   ",
                        elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf(ioQQQ," %9.2e", mole.sink[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                if( dynamics.lgAdvection )
                {
                        /* sinks from the dynamics */
                        fprintf(ioQQQ," %s dynm sink   ",
                                elementnames.chElementSym[nelem]);
                        for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                        {
                                fprintf(ioQQQ," %9.2e", dynamics.Rate );
                        }
                        fprintf( ioQQQ, "\n" );
                }

                /* sum of all creation processes */
                fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", ionbal.RateRecomTot[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* sources from the chemistry network */
                fprintf(ioQQQ," %s mole source ",
                        elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf(ioQQQ," %9.2e", mole.source[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                if( dynamics.lgAdvection )
                {
                        /* source from the dynamcs */
                        fprintf(ioQQQ," %s dynm sour   ",
                                elementnames.chElementSym[nelem]);
                        for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                        {
                                fprintf(ioQQQ," %9.2e", dynamics.Source[nelem][ion] );
                        }
                        fprintf( ioQQQ, "\n" );
                }

                /* collisional ionization */
                fprintf( ioQQQ, " %s Coll ioniz  " ,elementnames.chElementSym[nelem] );
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", ionbal.CollIonRate_Ground[nelem][ion][0] );
                }
                fprintf( ioQQQ, "\n" );

                /* UTA ionization */
                fprintf( ioQQQ, " %s UTA ioniz   " ,elementnames.chElementSym[nelem] );
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", ionbal.UTA_ionize_rate[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* photo ionization */
                fprintf( ioQQQ, " %s Phot ioniz  " ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", 
                                ionbal.PhotoRate_Shell[nelem][ion][Heavy.nsShells[nelem][ion]-1][0] );
                }
                fprintf( ioQQQ, "\n" );

                /* auger ionization */
                fprintf( ioQQQ, " %s Auger ioniz " ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", 
                                auger[ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* secondary ionization */
                fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", 
                                secondaries.csupra[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* grain ionization - not total rate but should be dominant process */
                fprintf( ioQQQ, " %s ion on grn  "  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion][ion+1] );
                }
                fprintf( ioQQQ, "\n" );

                /* chemistry ionization - not total rate but should be dominant process 
                 * not source or sink but just increase ionization within ladder */
                fprintf( ioQQQ, " %s ion in chem "  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion][ion+1] );
                }
                fprintf( ioQQQ, "\n" );

                /* charge exchange ionization */
                fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        /* sum has units s-1 */
                        double sum = atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 
                                atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];

                        if( nelem==ipHELIUM && ion==0 )
                        {
                                sum += atmdat.HeCharExcIonTotal;
                        }
                        else if( nelem==ipHYDROGEN && ion==0 )
                        {
                                sum += atmdat.HCharExcIonTotal;
                        }
                        fprintf( ioQQQ, " %9.2e", sum );
                }
                fprintf( ioQQQ, "\n" );

                /* radiative recombination */
                fprintf( ioQQQ, " %s radiati rec "  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.RR_rate_coef_used[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* old DR recombination */
                fprintf( ioQQQ, " %s dr old  rec "  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_old_rate_coef[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* Badnell DR recombination */
                fprintf( ioQQQ, " %s drBadnel rec"  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* grain recombination - not total but from next higher ion, should
                 * be dominant */
                fprintf( ioQQQ, " %s rec on grn  "  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion+1][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* grain due to chemistry - not total but from next higher ion, should
                 * be dominant - not source or sink, just moves up or down */
                fprintf( ioQQQ, " %s rec in chem "  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", mole.xMoleChTrRate[nelem][ion+1][ion] );
                }
                fprintf( ioQQQ, "\n" );

                /* charge exchange recombination */
                fprintf( ioQQQ, " %s chr trn rec "  ,elementnames.chElementSym[nelem]);
                for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
                {
                        double sum = 
                                atmdat.HCharExcRecTo[nelem][ion]*
                                StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*
                                dense.xIonDense[ipHYDROGEN][1] +

                                atmdat.HeCharExcRecTo[nelem][ion]*
                                StatesElem[ipHE_LIKE][ipHELIUM][ipHe1s1S].Pop*
                                dense.xIonDense[ipHELIUM][1];

                        if( nelem==ipHELIUM && ion==0 )
                        {
                                sum += atmdat.HeCharExcRecTotal;
                        }
                        else if( nelem==ipHYDROGEN && ion==0 )
                        {
                                sum += atmdat.HCharExcRecTotal;
                        }
                        fprintf( ioQQQ, " %9.2e", sum );
                }
                fprintf( ioQQQ, "\n" );

                /* the "new" abundances the resulted from the previous ratio */
                fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
                for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
                {
                        fprintf( ioQQQ, " %9.2e", dense.xIonDense[nelem][ion] );
                }
                fprintf( ioQQQ, "\n" );
        }

        if( PARALLEL_MODE )
        {
                free(ipiv);
                free(source);
                free(sink);
                free(xmat);
                free(xmatsave);
                free(auger);
        }
        return;
}

/* 

         Solve an ionization level system with specified ionization and
         recombination rates between neighboring ions, and additional sink
         and source terms.  The sink array is overwritten, and the results
         appear in the source array.

         Written in matrix form, the algorithm is equivalent to the
         tridiagonal algorithm in Numerical Recipes applied to:

         / i_0+a_0     -r_0          .           .    .  \ / x_0 \   / s_0 \
         |  -i_0    i_1+a_1+r_0    -r_1          .    .  | | x_1 |   | s_1 |
         |    .        -i_1      i_2+a_2+r_1   -r_2   .  | | x_2 |   | s_2 |
         |    .          .       (etc....)               | | ... | = | ... |
         \    .          .          .                    / \     /   \     /

         where i, r are the ionization and recombination rates, s is the
         source rate and a is the sink rate.

         This matrix is diagonally dominant only when the sink terms are
         large -- the alternative method coded here prevents rounding error
         in the diagonal terms disturbing the solution when this is not the
         case.

*/

/* solveions tridiagonal solver but optimized for structure of balance matrix */
void solveions(double *ion, double *rec, double *snk, double *src,
               long int nlev, long int nmax)
{
        double kap, bet;
        long int i;

        DEBUG_ENTRY( "solveions()" );

        if(nmax != -1) 
        {
                /* Singular case */
                src[nmax] = 1.;
                for(i=nmax;i<nlev-1;i++)
                        src[i+1] = src[i]*ion[i]/rec[i];
                for(i=nmax-1;i>=0;i--)
                        src[i] = src[i+1]*rec[i]/ion[i];
        } 
        else 
        {
                i = 0;
                kap = snk[0];    
                for(i=0;i<nlev-1;i++) 
                {
                        bet = ion[i]+kap;
                        if(bet == 0.)
                        {
                                fprintf(ioQQQ,"Ionization solver error\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                        bet = 1./bet;
                        src[i] *= bet;
                        src[i+1] += ion[i]*src[i];
                        snk[i] = bet*rec[i];
                        kap = kap*snk[i]+snk[i+1];
                }
                bet = kap;
                if(bet == 0.)
                {
                        fprintf(ioQQQ,"Ionization solver error\n");
                        cdEXIT(EXIT_FAILURE);
                }
                src[i] /= bet;

                for(i=nlev-2;i>=0;i--)
                {
                        src[i] += snk[i]*src[i+1];
                }
        }
}

---