conv_base.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ConvBase main routine to drive ionization solution for all species, find total opacity
 * called by ConvIoniz */
/*lgConverg check whether ionization of element nelem has converged */
#include "cddefines.h"
#include "dynamics.h"
#include "trace.h"
#include "elementnames.h"
#include "save.h"
#include "phycon.h"
#include "secondaries.h"
#include "stopcalc.h"
#include "grainvar.h"
#include "highen.h"
#include "dense.h"
#include "hmi.h"
#include "rfield.h"
#include "pressure.h"
#include "taulines.h"
#include "rt.h"
#include "grains.h"
#include "atmdat.h"
#include "ionbal.h"
#include "opacity.h"
#include "cooling.h"
#include "thermal.h"
#include "mole.h"
#include "iso.h"
#include "conv.h"
#include "h2.h"
#include "deuterium.h"

STATIC bool lgNetEdenSrcSmall( void );

void UpdateUTAs( void );

/*lgIonizConverg check whether ionization of element nelem has converged, called by ionize,
 * returns true if element is converged, false if not */
namespace
{
        class IonizConverg
        {
                /* this is fractions [ion stage][nelem], ion stage = 0 for atom, nelem=0 for H*/
                double OldFracs[LIMELM+1][LIMELM+1];
        public:
                IonizConverg()
                        {
                                for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
                                {
                                        for (long ion = 0; ion <= nelem+1; ++ion)
                                        {
                                                OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
                                        }
                                }
                        }
                void operator()(
                        /* atomic number on the C scale, 0 for H, 25 for Fe */
                        long loop_ion ,
                        /* this is allowed error as a fractional change.  Most are 0.15 */
                        double delta ,
                        /* option to print abundances */
                        bool lgPrint )
                {
                        bool lgConverg_v = false;
                        double Abund, 
                                bigchange ,
                                change ,
                                one;
                        
                        DEBUG_ENTRY( "IonizConverg::operator()" );
                        
                        for (long nelem =ipHYDROGEN; nelem<LIMELM; ++nelem )
                        {
                                if( !dense.lgElmtOn[nelem] )
                                        continue;
                                
                                double abundold=0. , abundnew=0.;
                                long ionchg=-1;
                                
                                /* arguments are atomic number, ionization stage
                                 * OldFracs[nelem][0] is abundance of nelem (cm^-3) */
                                
                                /* this function returns true if ionization of element
                                 * with atomic number nelem has not changed by more than delta,*/
                                
                                /* check whether abundances exist yet, only do this for after first zone */
                                /*if( nzone > 0 )*/
                                /* >>chng 03 sep 02, check on changed ionization after first call through,
                                 * to insure converged constant eden / temperature models */
                                if( conv.nPres2Ioniz )
                                {
                                        /* >>chng 04 aug 31, this had been static, caused last hits on large changes
                                         * in ionization */
                                        bigchange = 0.;
                                        change = 0.;
                                        Abund = dense.gas_phase[nelem];
                                        
                                        /* loop over all ionization stages, loop over all ions, not just active ones,
                                         * since this also sets old values, and so will zero out non-existant ones */
                                        for( long ion=0; ion <= (nelem+1); ++ion )
                                        {
                                                /*lint -e727 OlsdFracs not initialized */
                                                if( OldFracs[nelem][ion]/Abund > 1e-4 && 
                                                         dense.xIonDense[nelem][ion]/Abund > 1e-4 )
                                                        /*lint +e727 OlsdFracs not initialized */
                                                {
                                                        /* check change in old to new value */
                                                        one = fabs(dense.xIonDense[nelem][ion]-OldFracs[nelem][ion])/
                                                        OldFracs[nelem][ion];
                                                        change = MAX2(change, one );
                                                        /* remember abundances for largest change */
                                                        if( change>bigchange )
                                                        {
                                                                bigchange = change;
                                                                abundold = OldFracs[nelem][ion]/Abund;
                                                                abundnew = dense.xIonDense[nelem][ion]/Abund;
                                                                ionchg = ion;
                                                        }
                                                }
                                                /* now set new value */
                                                OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
                                        }
                                        
                                        if( change >= delta )
                                        {
                                                char chConvIoniz[INPUT_LINE_LENGTH];
                                                sprintf( chConvIoniz , "%2s ion" , elementnames.chElementSym[nelem] );
                                                conv.setConvIonizFail(chConvIoniz,abundold,abundnew);
                                                ASSERT( abundold>0. && abundnew>0. );
                                        }
                                        else
                                                lgConverg_v = true;
                                }
                                else
                                {
                                        for( long ion=0; ion <= (nelem+1); ++ion )
                                        {
                                                OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
                                        }
                                        
                                }
                                
                                /* option to print abundances */
                                if( lgPrint )
                                {
                                        fprintf(ioQQQ," nz %ld loop %ld element %li converged? %c worst %ld change %g\n",
                                                          nzone, loop_ion, nelem, TorF(lgConverg_v),ionchg,bigchange);
                                        for( long ion=0; ion<(nelem+1); ++ion )
                                        {
                                                fprintf(ioQQQ,"\t%.5e", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]);
                                        }
                                        fprintf(ioQQQ,"\n");
                                }
                        }
                }
        };
}

/*ConvBase main routine to drive ionization solution for all species, find total opacity
 * called by ConvIoniz
 * return 0 if ok, 1 if abort */
int ConvBase( 
        /* this is zero the first time ConvBase is called by convIoniz,
         * counts number of call thereafter */
         long loopi )
{
        double HeatOld,
                EdenTrue_old,
                EdenFromMolecOld,
                EdenFromGrainsOld,
                HeatingOld ,
                CoolingOld;
        static double SecondOld;
        static long int nzoneOTS=-1;
#       define LOOP_ION_LIMIT   10
        long int loop_ion;
        static double SumOTS=0. , OldSumOTS[2]={0.,0.};
        double save_iso_grnd[NISO][LIMELM];
        valarray<realnum> mole_save(mole_global.num_calc);
        IonizConverg lgIonizConverg;

        DEBUG_ENTRY( "ConvBase()" );

        /* this is set to phycon.te in tfidle, is used to insure that all temp
         * vars are properly updated when conv_ionizeopacitydo is called 
         * NB must be same type as phycon.te */
        ASSERT( fp_equal( phycon.te, thermal.te_update ) );

        /* this allows zone number to be printed with slight increment as zone converged
         * conv.nPres2Ioniz is incremented at the bottom of this routine */
        fnzone = (double)nzone + (double)conv.nPres2Ioniz/100.;

        /* reevaluate pressure */
        /* this sets values of pressure.PresTotlCurr, also calls tfidle,
         * and sets the total energy content of gas, which may be important for acvection */
        PresTotCurrent();

        /* >>chng 04 sep 15, move EdenTrue_old up here, and will redo at bottom
         * to find change 
         * find and save current true electron density - if this changes by more than the
         * tolerance then ionization solution is not converged */
        /* >>chng 04 jul 27, move eden_sum here from after this loop, so that change in EdenTrue
         * can be monitored */
        /* update EdenTrue, eden itself is actually changed in ConvEdenIoniz */
        /* must not call eden_sum on very first time since for classic PDR total
         * ionization may still be zero on first call */
        if( conv.nTotalIoniz )
        {
                if( eden_sum() )
                {
                        /* non-zero return indicates abort condition */
                        ++conv.nTotalIoniz;
                        return 1;
                }
        }

        /* the following two were evaluated in eden_sum 
         * will confirm that these are converged */
        EdenTrue_old = dense.EdenTrue;
        EdenFromMolecOld = mole.elec;
        EdenFromGrainsOld = gv.TotalEden;
        HeatingOld = thermal.htot;
        CoolingOld = thermal.ctot;

        /* remember current ground state population - will check if converged */
        for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                for( long nelem=ipISO; nelem<LIMELM;++nelem )
                {
                        if( dense.lgElmtOn[nelem] )
                        {
                                /* save the ground state population */
                                save_iso_grnd[ipISO][nelem] = iso_sp[ipISO][nelem].st[0].Pop();
                        }
                }
        }

        for( long i=0; i < mole_global.num_calc; i++ )
        {
                mole_save[i] = (realnum) mole.species[i].den;
        }

        if( loopi==0 )
        {
                /* these will be used to look for oscillating ots rates */
                OldSumOTS[0] = 0.;
                OldSumOTS[1] = 0.;
                conv.lgOscilOTS = false;
        }

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, 
                        "   ConvBase called. %.2f Te:%.3e  HI:%.3e  HII:%.3e  H2:%.3e  Ne:%.3e  htot:%.3e  CSUP:%.2e Conv?%c\n", 
                  fnzone,
                  phycon.te, 
                  dense.xIonDense[ipHYDROGEN][0], 
                  dense.xIonDense[ipHYDROGEN][1], 
                  hmi.H2_total,
                  dense.eden, 
                  thermal.htot, 
                  secondaries.csupra[ipHYDROGEN][0] ,
                  TorF(conv.lgConvIoniz()) );
        }
        /* want this flag to be true when we exit, various problems will set falst */
        conv.resetConvIoniz();

        /* this routine is in heatsum.c, and zeros out the heating array */
        HeatZero();

        /* if this is very first call, say not converged, since no meaningful way to
         * check on changes in quantities.  this counter is false only on very first
         * call, when equal to zero */
        if( !conv.nTotalIoniz )
        {
                conv.setConvIonizFail( "first call", 0.0, 0.0 );
        }

        /* this will be flag to check whether any ionization stages
         * were trimmed down */
        conv.lgIonStageTrimed = false;

        /* must redo photoionization rates for all species on second and later tries */
        /* always reevaluate the rates when . . . */
        /* this flag says not to redo opac and photo rates, and following test will set
         * true if one of several tests not done*/
        opac.lgRedoStatic = false;
        if( 
                /* opac.lgOpacStatic (usually true), is set false with no static opacity command */
                !opac.lgOpacStatic || 
                /* we are in search mode */
                conv.lgSearch || 
                /* this is the first call to this zone */
                conv.nPres2Ioniz == 0 )
        {
                /* we need to redo ALL photoionization rates */
                opac.lgRedoStatic = true;
        }

        /* calculate radiative and dielectronic recombination rate coefficients */
        ion_recom_calculate();

        /* this adjusts the lowest and highest stages of ionization we will consider,
         * only safe to call when lgRedoStatic is true since this could lower the 
         * lowest stage of ionization, which needs all its photo rates */

        /* conv.nTotalIoniz is only 0 (false) on the very first call to here,
         * when the ionization distribution is not yet done */
        if( conv.nTotalIoniz )
        {
                bool lgIonizTrimCalled = false;
                static long int nZoneCalled = 0;

                fixit(); // nZoneCalled should be reinitialized for each grid point?

                /* ionization trimming only used for He and heavier, not H */
                /* only do this one time per zone since up and down cycle can occur */
                /* >>chng 05 jan 15, increasing temperature above default first conditions, also
                 * no trim during search - this fixed major logical error when sim is
                 * totally mechanically heated to coronal temperatures -
                 * problem discovered by Ronnie Hoogerwerf */
                /* do not keep changing the trim after the first call within
                 * this zone once we are deep into layer - doing so introduced very 
                 * small level of noise as some stages
                 * went up and down - O+2 - O+3 was good example, when small H+3 after He+ i-front 
                 * limit to one increase per element per zone */
                if( conv.nTotalIoniz>2 && 
                        /* only call one time per zone except during search phase,
                         * when only call after 20 times only if temperature has changed */
                        ( conv.lgSearch || nZoneCalled!=nzone) )
                {
                        lgIonizTrimCalled = true;
                        for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem )
                        {
                                if( dense.lgElmtOn[nelem] )
                                {
                                        /* ion_trim will set conv.lgIonStageTrimed true is any ion has its
                                         * lowest stage of ionization dropped or trimmed */
                                        ion_trim(nelem);
                                }
                        }
                        nZoneCalled = nzone;
                }

                /* following block only set of asserts */
#               if !defined(NDEBUG)
                /* check that proper abundances are either positive or zero */
                for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem)
                {
                        if( dense.lgElmtOn[nelem] )
                        {
                                for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
                                {
                                        ASSERT( dense.xIonDense[nelem][ion] == 0. );
                                }
                                /*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/
                                for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
                                {
                                        /* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to
                                        * trip over VERY small floats that failed on alphas, but not 386
                                        * 
                                        * in case where lower ionization stage was just lowered or
                                        * trimmed down the abundance
                                        * was set to SMALLFLOAT so test must be < SMALLFLOAT */
                                        /* >>chng 02 feb 19, add check for search phase.  During this search
                                        * models with extreme ionization (all neutral or all ionized) can
                                        * have extreme but non-zero abundances far from the ionization peak for
                                        * element with lots of electrons.  These will go away once the model
                                        * becomes stable */
                                        /* >>chng 03 dec 01, add check on whether ion trim was called 
                                         * conserve.in threw assert when iontrim not called and abund grew small */
                                        ASSERT( conv.lgSearch || !lgIonizTrimCalled ||
                                                /* this can happen if all C is in the form of CO */
                                                (ion==0 && dense.IonHigh[nelem]==0 ) ||
                                                dense.xIonDense[nelem][ion] >= SMALLFLOAT || 
                                                dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT );
                                }
                                for( long ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
                                {
                                        ASSERT( ion >= 0 );
                                        ASSERT( dense.xIonDense[nelem][ion] == 0. );
                                }
                        }
                }
#               endif
        }

        /* now check if anything trimmed down */
        if( conv.lgIonStageTrimed )
        {
                /* something was trimmed down, so say that ionization not yet stable */
                /* say that ionization has not converged, secondaries changed by too much */
                conv.setConvIonizFail( "IonTrimmed", 0.0, 0.0 );
        }

        /* reevaluate advective terms if turned on */
        if( dynamics.lgAdvection )
                DynaIonize();
        
        /* evaluate Compton heating, bound E Compton, cosmic rays */
        highen();

        // Depends on dense.eden, phycon.te and trimming level of H, He
        iso_collapsed_update();

        iso_update_rates( );

        long ionLowCache[LIMELM], ionHighCache[LIMELM];
        for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
        {
                if( dense.lgSetIoniz[nelem] )
                {
                        dense.IonLow[nelem] = 0;
                        dense.IonHigh[nelem] = nelem + 1;
                        while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
                                ++dense.IonLow[nelem];
                        while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
                                --dense.IonHigh[nelem];
                }
                ionLowCache[nelem] = dense.IonLow[nelem];
                ionHighCache[nelem] = dense.IonHigh[nelem];
        }

        /* >>chng 04 feb 15, add loop over ionization until converged.  Non-convergence
         * in this case means ionization ladder pop changed, probably because of way
         * that Auger ejection is treated - loop exits when conditions tested at end are met */
        conv.incrementCounter(CONV_BASE_CALLS);
        loop_ion = 0;
        do
        {
                conv.incrementCounter(CONV_BASE_LOOPS);
                ASSERT(lgElemsConserved());
        
                /* set the convergence flag to true, 
                 * if anything changes in ConvBase, it will be set false */
                if( loop_ion )
                        conv.resetConvIoniz();

                /* charge transfer evaluation needs to be here so that same rate 
                 * coefficient used for H ion and other recombination */
                /* fill in master charge transfer array, and zero out some arrays that track effects */

                ChargTranEval();

                /* find grain temperature, charge, and gas heating rate */
                /* >>chng 04 apr 15, moved here from after call to HeLike(), PvH */
                GrainDrive();

                /* evaluate Compton heating, bound E Compton, cosmic rays */
                highen();

                /* find corrections for three-body rec - collisional ionization */
                atmdat_3body();

                /* update UTA inner-shell ionization rates */
                UpdateUTAs();

                ASSERT(lgElemsConserved());
        
                /* evaluate molecular hydrogen level populations */
                for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
                {
                        bool lgPopsConverged = true;
                        double old_val, new_val;
                        (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
                        if( !lgPopsConverged )
                        {
                                conv.setConvIonizFail( "H2 pops", old_val, new_val);
                        }       
                }

                ASSERT(lgElemsConserved());
                
                long ion_loop = 0;
                const int nconv = 3;
                double xIonDense0[nconv][LIMELM][LIMELM+1];
                bool lgShortCircuit = false;
                for( ion_loop=0; ion_loop<nconv && !lgShortCircuit; ++ion_loop)
                {
                        conv.incrementCounter(CONV_BASE_ACCELS);
                        for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
                        {
                                for (long ion = 0; ion <= nelem+1; ++ion)
                                {
                                        xIonDense0[ion_loop][nelem][ion] = dense.xIonDense[nelem][ion];
                                }
                        }

                        double netion = 0.0;
                        for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
                        {
                                if (!dense.lgElmtOn[nelem])
                                        continue;
                                
                                ASSERT(dense.IonLow[nelem] >= ionLowCache[nelem]);
                                ASSERT(dense.IonHigh[nelem] <= ionHighCache[nelem]);
                                ion_wrapper( nelem );
                                if (nelem == ipHYDROGEN || nelem == ipOXYGEN)
                                {
                                        double hion = 
                                                atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*
                                                iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*
                                                dense.xIonDense[ipOXYGEN][1]-
                                                atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*
                                                dense.xIonDense[ipHYDROGEN][1]*
                                                dense.xIonDense[ipOXYGEN][0];
                                        if (nelem == ipHYDROGEN)
                                                netion += hion;
                                        else
                                                netion -= hion;
                                }
                                
                                if ( conv.lgUpdateCouplings )
                                        mole_update_sources();
                        }
                        if (dense.lgElmtOn[ipHYDROGEN] && dense.lgElmtOn[ipOXYGEN] && ion_loop == 1)
                        {
                                // Comparison rate is whether error is either a small fraction of the minimum ionization rate, 
                                // or a tighter check on whether the fluxes are balanced as well as is feasible numerically
                                double ion_cmp = MAX2(0.01*MIN2(ionbal.RateIonizTot(ipHYDROGEN,0), ionbal.RateIonizTot(ipOXYGEN,0)),
                                                                                         1e-4*atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*
                                                                                         iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*dense.xIonDense[ipOXYGEN][1]);

                                if (fabs(netion) > ion_cmp && 
                                         dense.xIonDense[ipOXYGEN][1]*dense.xIonDense[ipHYDROGEN][1] >
                                         1e-8*iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*dense.xIonDense[ipOXYGEN][0] )
                                {
                                        conv.setConvIonizFail( "OH CX inconsistency" , ion_cmp, netion);
                                }
                        }

                        // If not going to end anyhow, check whether changes were small
                        bool lgCanShortCircuit = (ion_loop+1 < nconv);
                        for( long nelem=ipHYDROGEN; nelem<LIMELM && lgCanShortCircuit; ++nelem )
                        {
                                if (!dense.lgElmtOn[nelem])
                                        continue;
                                for (long ion = dense.IonLow[nelem]; 
                                          ion <= dense.IonHigh[nelem]; ++ion)
                                {
                                        double x0 = xIonDense0[ion_loop][nelem][ion];
                                        double x1 = dense.xIonDense[nelem][ion];
                                        if (fabs(x0-x1) > 1e-6*(x0+x1))
                                        {
                                                lgCanShortCircuit = false;
                                                break;
                                        }
                                }
                        }
                        lgShortCircuit = lgCanShortCircuit;
                }

                if (!lgShortCircuit)
                {
                        // Apply convergence acceleration
                        for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
                        {
                                if (!dense.lgElmtOn[nelem])
                                        continue;
                                double tot0 = 0., tot1 = 0.;
                                double xIonNew[LIMELM+1];
                                for (long ion = 0; ion <= nelem+1; ++ion)
                                {
                                        double x0 = xIonDense0[nconv-2][nelem][ion];
                                        double x1 = xIonDense0[nconv-1][nelem][ion];
                                        double x2 = dense.xIonDense[nelem][ion];
                                        xIonNew[ion] = x2;
                                        tot0 += x2;
                                        // Richardson extrapolation formula to accelerate convergence
                                        // Assumes convergence to x^ is geometric, i.e. 
                                        // x_i = x^ + a d^i
                                        
                                        double extstep = 0.,predict=x2,
                                                step0 = x1-x0, step1 = x2-x1, abs1 = fabs(step1);
                                        // Protect against roundoff/noise in inner solver
                                        if ( abs1 > 1000.0*((double)DBL_EPSILON)*x2 ) 
                                        {
                                                double denom = fabs(step1-step0);
                                                double sgn = (step1*step0 > 0)? 1.0 : -1.0;
                                                // Greatest acceleration allowed is MAXACC*latest step length
                                                // Can we do better than static tuning this parameter?
                                                const double MAXACC=100.0; 
                                                double extfac = 1.0/(denom/abs1 + 1.0/MAXACC);
                                                extstep = sgn*extfac*step1;
                                                //extstep = sgn*MIN2(extfac*step1,0.01*x2);
                                                predict = x2+extstep;
                                                if (predict > 0.0)
                                                        xIonNew[ion] = predict;
                                        }
                                        if ( 0 )
                                                //if ( nelem == ipHYDROGEN || (nelem == ipIRON && ion <=2 ) )
                                                if ( (nelem == ipNICKEL && ion <=2 ) )
                                                        fprintf(ioQQQ,"Extrap %3ld %3ld %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
                                                                          nelem,ion,
                                                                          x0,x0-xIonDense0[nconv-3][nelem][ion],x1-x0,x2-x1,extstep,predict);
                                        tot1 += xIonNew[ion];
                                }
                                if ( tot1 > SMALLFLOAT )
                                {
                                        double scal = tot0/tot1;
                                        for (long ion = 0; ion <= nelem+1; ++ion)
                                        {
                                                dense.xIonDense[nelem][ion] = scal*xIonNew[ion];
                                        }
                                        for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
                                        {
                                                double renorm;
                                                iso_renorm(nelem, ipISO, renorm);
                                        }
                                }
                        }

                        bool lgPostExtrapSolve = true;
                        if (lgPostExtrapSolve)
                        {
                                for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
                                {
                                        if (!dense.lgElmtOn[nelem])
                                                continue;
                                        ion_wrapper( nelem );
                                }
                        }
                }
                SetDeuteriumIonization( dense.xIonDense[ipHYDROGEN][0], dense.xIonDense[ipHYDROGEN][1] );

                /*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures
                 * which cascaded downstream if we did not abort */
                /* return if one of above solvers has declared abort condition */
                if( lgAbort )
                {
                        ++conv.nTotalIoniz;
                        return 1;
                }
                
                ASSERT(lgElemsConserved());
        
                /* drive chemistry network*/
                mole_drive();

                ASSERT(lgElemsConserved());
        
                /* all elements have now had ionization reevaluated - in some cases we may have upset
                 * the ionization that was forced with an "element ionization" command - here we will 
                 * re-impose that set ionization */
                /* >>chng 04 oct 13, add this logic */
                for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
                {
                        if( dense.lgSetIoniz[nelem] )
                        {
                                dense.IonLow[nelem] = 0;
                                dense.IonHigh[nelem] = nelem + 1;
                                while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
                                        ++dense.IonLow[nelem];
                                while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
                                        --dense.IonHigh[nelem];
                        }
                }

                /* redo ct ion rate for reasons now unclear */
                ChargTranEval();

                /* lgIonizConverg is a function to check whether ionization has converged
                 * check whether ionization changed by more than relative error
                 * given by second number */
                /* >>chng 04 feb 14, loop over all elements rather than just a few */
                lgIonizConverg(loop_ion, conv.IonizErrorAllowed  , false );

                if( deut.lgElmtOn )
                {
                        static double OldDeut[2] = {0., 0.};
                        for( long ion=0; ion<2; ++ion )
                        {
                                if( fabs(deut.xIonDense[ion] - OldDeut[ion] ) > 0.2*conv.IonizErrorAllowed*fabs(OldDeut[ion]) )
                                {
                                        conv.setConvIonizFail( "D ion" , OldDeut[ion], deut.xIonDense[ion]);
                                }
                                OldDeut[ion] = deut.xIonDense[ion];
                        }
                }

                if( fabs(EdenTrue_old - dense.EdenTrue) > conv.EdenErrorAllowed/2.*fabs(dense.EdenTrue) )
                {
                        conv.setConvIonizFail( "EdTr cng" , EdenTrue_old, dense.EdenTrue);
                }

                for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
                {
                        for( long nelem=ipISO; nelem<LIMELM; ++nelem )
                        {
                                lgStatesConserved( nelem, nelem-ipISO, iso_sp[ipISO][nelem].st, iso_sp[ipISO][nelem].numLevels_local, 1e-3f, loop_ion );
                        }
                }
        
                if( trace.nTrConvg>=4 )
                {
                        /* trace ionization  */
                        fprintf( ioQQQ, 
                                "    ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" ,
                                loop_ion , 
                                TorF( conv.lgConvIoniz()) ,
                                conv.chConvIoniz() );
                }

                ++loop_ion;
        }
        /* loop is not converged, less than loop limit, and we are reevaluating */
        while( !conv.lgConvIoniz() && loop_ion < LOOP_ION_LIMIT && rfield.lgIonizReevaluate);

        if( conv.lgConvIoniz() )
                lgNetEdenSrcSmall();

        /* >>chng 05 oct 29, move CT heating here from heat_sum since this sometimes contributes
         * cooling rather than heat and so needs to be sorted out before either heating or cooling 
         * are derived first find net heating - */
        thermal.char_tran_heat = ChargTranSumHeat();
        /* net is cooling if negative */
        thermal.char_tran_cool = MAX2(0. , -thermal.char_tran_heat );
        thermal.char_tran_heat = MAX2(0. ,  thermal.char_tran_heat );

        /* get total cooling, thermal.ctot = does not occur since passes as pointer.  This can add heat.
         * it calls coolSum at end to sum up the total cooling */
        CoolEvaluate(&thermal.ctot );

        /* get total heating rate - first save old quantities to check how much it changes */
        HeatOld = thermal.htot;

        /* HeatSum will update ElecFrac, 
         * secondary electron ionization and excitation efficiencies,
         * and sum up current secondary rates - remember old values before we enter */
        SecondOld = secondaries.csupra[ipHYDROGEN][0];

        /* update the total heating - it was all set to zero in HeatZero at top of this routine,
         * occurs before secondaries bit below, since this updates electron fracs */
        HeatSum();

        /* test whether we can just set the rate to the new one, or whether we should
         * take average and do it again.  secondaries.sec2total was set in hydrogenic, and
         * is ratio of secondary to total hydrogen destruction rates */
        /* >>chng 02 nov 20, add test on size of csupra - primal had very close to underflow */
        if( (secondaries.csupra[ipHYDROGEN][0]>SMALLFLOAT) && (secondaries.sec2total>0.001) && 
                fabs( 1. - SecondOld/SDIV(secondaries.csupra[ipHYDROGEN][0]) ) > 0.1 && 
                SecondOld > 0. &&  secondaries.csupra[ipHYDROGEN][0] > 0.)
        {
                /* say that ionization has not converged, secondaries changed by too much */
                conv.setConvIonizFail( "SecIonRate", SecondOld,
                                                                          secondaries.csupra[ipHYDROGEN][0] );
        }

#       if 0
        static realnum hminus_old=0.;
        /* >>chng 04 apr 15, add this convergence test */
        if( conv.nTotalIoniz )
        {
                realnum hminus_den = findspecieslocal("H-")->den; 
                if( fabs( hminus_old-hminus_den ) > fabs( hminus_den * conv.EdenErrorAllowed ) )
                {
                        conv.setConvIonizFail( "Big H- chn", hminus_old, hminus_den );
                }
                hminus_old = hminus_den;
        }
#       endif

        if( HeatOld > 0. && !thermal.lgTemperatureConstant )
        {
                /* check if heating has converged - tolerance is final match */
                if( fabs(1.-thermal.htot/HeatOld) > conv.HeatCoolRelErrorAllowed*.5 )
                {
                        conv.setConvIonizFail( "Big d Heat", HeatOld, thermal.htot);
                }
        }

        /* abort flag may have already been set - if so bail */
        if( lgAbort )
        {

                return 1;
        }

        /* evaluate current opacities 
         * rfield.lgOpacityReevaluate normally true, 
         * set false with no opacity reevaluate command, an option to only 
         * evaluate opacity one time per zone */
        if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate  || conv.lgIonStageTrimed )
                OpacityAddTotal();

        /* >>chng 02 jun 11, call even first time that this routine is called -
         * this seemed to help convergence */

        /* do OTS rates for all lines and all continua since
        * we now have ionization balance of all species.  Note that this is not
        * entirely self-consistent, since destruction probabilities here are not the same as
        * the ones used in the model atoms.  Problems??  if near convergence
        * then should be nearly identical */
        if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || rfield.lgIonizReevaluate ||
                conv.lgIonStageTrimed || conv.lgSearch || nzone!=nzoneOTS )
        {
                RT_OTS();
                nzoneOTS = nzone;

                /* remember old ots rates */
                OldSumOTS[0] = OldSumOTS[1];
                OldSumOTS[1] = SumOTS;
                /*fprintf(ioQQQ," calling RT_OTS zone %.2f SumOTS is %.2e\n",fnzone,SumOTS);*/

                /* now update several components of the continuum, this only happens after
                 * we have gone through entire solution for this zone at least one time. 
                 * there can be wild ots oscillation on first call */
                /* the rel change of 0.2 was chosen by running hizqso - smaller increased
                 * itrzn but larger did not further decrease it. */
                RT_OTS_Update(&SumOTS);
                /*fprintf(ioQQQ,"RT_OTS_Updateee\t%.3f\t%.2e\t%.2e\n", fnzone,SumOTS , OldSumOTS[1] );*/
        }
        else
                SumOTS = 0.;

        /* now check whether the ots rates changed */
        if( SumOTS> 0. )
        {
                /* the ots rate must be converged to the error in the electron density */
                /* how fine should this be converged?? originally had above, 10%, but take
                 * smaller ratio?? */
                if( fabs(1.-OldSumOTS[1]/SumOTS) > conv.EdenErrorAllowed )
                {
                        /* this branch, ionization not converged due to large change in ots rates.
                         * check whether ots rates are oscillating, if loopi > 1 so we have enough info*/
                        if( loopi > 1 )
                        {
                                /* here we have three values, are they changing sign? */
                                if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
                                {
                                        /* ots rates are oscillating */
                                        conv.lgOscilOTS = true;
                                }
                        }

                        conv.setConvIonizFail( "OTSRatChng" , OldSumOTS[1], SumOTS);
                }

                /* produce info on the ots fields if either "trace ots" or 
                 * "trace convergence xxx ots " was entered */
                if( ( trace.lgTrace && trace.lgOTSBug ) ||
                        ( trace.nTrConvg && trace.lgOTSBug ) )
                {
                        RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
                }
                /*fprintf(ioQQQ,"DEBUG opac\t%.2f\t%.3e\t%.3e\n",fnzone,
                        dense.xIonDense[ipNICKEL][0] ,
                        dense.xIonDense[ipZINC][0] );*/
                {
                        /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
                        enum {DEBUG_LOC=false};
                        if( DEBUG_LOC && (nzone>110) )
                        {
#                               if 0
#                               include "lines_service.h"
                                DumpLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(2,0) );
#                               endif
                                /* last par 'l' for line, 'c' for continua, 'b' for both,
                                 * the numbers printed are:
                                 * cell i, [i], so 1 less than ipoint
                                 * anu[i], 
                                 * otslin[i], 
                                 * opacity_abs[i],
                                 * otslin[i]*opacity_abs[i],
                                 * rfield.chLineLabel[i] ,
                                 * rfield.line_count[i] */
                        }
                }
        }
        {
                /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
                enum {DEBUG_LOC=false};
                if( DEBUG_LOC && (nzone>200) )
                {
                        fprintf(ioQQQ,"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n",
                                nzone,
                                iso_sp[0][1].trans(15,3).Emis().ots(),
                                TauLines[26].Emis().ots(),
                                opac.opacity_abs[2069]);
                }
        }

        /* option to print OTS continuum with TRACE OTS */
        if( trace.lgTrace && trace.lgOTSBug )
        {
                /* find ots rates here, so we only print fraction of it,
                 * SumOTS is both line and continuum contributing to ots, and is multiplied by opacity */
                /* number is weakest rate to print */
                RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
        }

        // all RT routines called 
        RT_line_all(  );

        /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
        /* reevaluate pressure */
        /* this sets values of pressure.PresTotlCurr, also calls tfidle */
        PresTotCurrent();

        ASSERT(lgElemsConserved());

        /* update some counters that keep track of how many times this routine
         * has been called */

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, 
                        "   ConvBase return. fnzone %.2f nPres2Ioniz %li Te:%.3e  HI:%.3e  HII:%.3e  H2:%.3e  Ne:%.3e  htot:%.3e  CSUP:%.2e Conv?%c reason:%s\n", 
                        fnzone,
                        conv.nPres2Ioniz,
                        phycon.te, 
                        dense.xIonDense[ipHYDROGEN][0], 
                        dense.xIonDense[ipHYDROGEN][1], 
                        hmi.H2_total,
                        dense.eden, 
                        thermal.htot, 
                        secondaries.csupra[ipHYDROGEN][0] ,
                        TorF(conv.lgConvIoniz()) ,
                        conv.chConvIoniz());
        }

        /* this counts number of times we are called by ConvPresTempEdenIoniz, 
         * number of calls in this zone so first call is zero
         * reset to zero each time ConvPresTempEdenIoniz is called */
        ++conv.nPres2Ioniz;

        /* this is abort option set with SET PRESIONIZ command,
         * test on nzone since init can take many iterations
         * this is seldom used except in special cases */
        if( conv.nPres2Ioniz > conv.limPres2Ioniz && nzone > 0)
        {
                fprintf(ioQQQ,"PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz. ");
                fprintf(ioQQQ,"Their values are %li and %li.\n",conv.nPres2Ioniz , conv.limPres2Ioniz);
                lgAbort = true;
                return 1;
        }

        /* various checks on the convergence of the current solution */
        if( eden_sum() )
        {
                /* non-zero return indicates abort condition */
                return 1;
        }

        /* is electron density converged? */
        if( fabs(EdenTrue_old - dense.EdenTrue) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
        {
                conv.setConvIonizFail( "eden chng", EdenTrue_old, dense.EdenTrue);
        }

        /* check on molecular electron den */
        if( fabs(EdenFromMolecOld - mole.elec) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
        {
                conv.setConvIonizFail( "edn chnCO", EdenFromMolecOld, dense.EdenTrue);
        }

        if( gv.lgGrainElectrons )
        {
                /* check on grain electron den */
                if( fabs(EdenFromGrainsOld - gv.TotalEden) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
                {
                        conv.setConvIonizFail( "edn grn e", EdenFromGrainsOld, gv.TotalEden);
                }

                /* check on sum of grain and molecular electron den - often two large numbers that cancel */
                if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (mole.elec-gv.TotalEden) ) > 
                        fabs(dense.EdenTrue * conv.EdenErrorAllowed/4.) )
                {
                        conv.setConvIonizFail( "edn mole-grn", 
                                                                                  (EdenFromMolecOld-EdenFromGrainsOld),
                                                                                  (mole.elec-gv.TotalEden));
                }
        }

        /* check on heating and cooling if vary temp model 
         * >>chng 08 jul 01, over the code's entire history it had tested whether
         * this is a constant temperature simulation and did not do this test if
         * the thermal solution was not done.  There are some cases where we do
         * want to specify the temperature and then find the heating or cooling -
         * this is done in calculations of cooling curves for instance.  With this
         * change the heating/cooling are converged even in a constant temperature
         * sim.  this does make CT sims run more slowly but with greater accuracy
         * if heating or cooling is reported */
        if( !thermal.lgTemperatureConstant )
        {
                if( fabs(HeatingOld - thermal.htot)/thermal.htot > conv.HeatCoolRelErrorAllowed/2. )
                {
                        conv.setConvIonizFail( "heat chg", HeatingOld, thermal.htot);
                }

                if( fabs(CoolingOld - thermal.ctot)/thermal.ctot > conv.HeatCoolRelErrorAllowed/2. )
                {
                        conv.setConvIonizFail( "cool chg", CoolingOld, thermal.ctot);
                }
        }

        /* check whether molecular abundances are stable */
        for( long i=0; i < mole_global.num_calc; ++i )
        {
                // done relative to total nuclei density so that we can do pure metal plasmas.
                if( fabs(mole.species[i].den-mole_save[i])/dense.xNucleiTotal-1. >
                        conv.EdenErrorAllowed/2.)
                {
                        char chConvIoniz[INPUT_LINE_LENGTH];
                        sprintf( chConvIoniz, "ch %-4.4s",mole_global.list[i]->label.c_str() );
                        conv.setConvIonizFail( chConvIoniz,
                                                                                  mole_save[i]/dense.xNucleiTotal,
                                                                                  mole.species[i].den/dense.xNucleiTotal);
                }
        }

        /* >>chng 05 mar 26, add this convergence test - important for molecular or advective
         * sims since iso ion solver must sync up with chemistry */
        /* remember current ground state population - will check if converged */
        for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                for( long nelem=ipISO; nelem<LIMELM;++nelem )
                {
                        if( dense.lgElmtOn[nelem] )
                        {
                                /* only do this check for "significant" levels of ionization */
                                /*lint -e644 var possibly not init */
                                if( dense.xIonDense[nelem][nelem-ipISO]/dense.gas_phase[nelem] > 1e-5 )
                                {
                                        if( fabs(iso_sp[ipISO][nelem].st[0].Pop()-save_iso_grnd[ipISO][nelem])/SDIV(iso_sp[ipISO][nelem].st[0].Pop())-1. >
                                                conv.EdenErrorAllowed)
                                        {
                                                char chConvIoniz[INPUT_LINE_LENGTH];
                                                sprintf( chConvIoniz,"iso %2li %2li",ipISO, nelem );
                                                conv.setConvIonizFail(chConvIoniz,
                                                                                                         save_iso_grnd[ipISO][nelem],
                                                                                                         iso_sp[ipISO][nelem].st[0].Pop());
                                        }
                                }
                                /*lint +e644 var possibly not init */
                        }
                }
        }

        /* this counts how many times ConvBase has been called in this iteration,
         * located here at bottom of routine so that number is false on first
         * call, set to 0 in when iteration starts - used to create itr/zn 
         * number in printout often used to tell whether this is the very 
         * first attempt at solution in this iteration */
        ++conv.nTotalIoniz;

        /* first sweep is over if this is not search phase */
        if( !conv.lgSearch )
                conv.lgFirstSweepThisZone = false;

        /* this was set with STOP NTOTALIONIZ command - only a debugging aid 
         * by default is zero so false */
        if(StopCalc.nTotalIonizStop && conv.nTotalIoniz> StopCalc.nTotalIonizStop )
        {
                /* non-zero return indicates abort condition */
                fprintf(ioQQQ , "ABORT flag set since STOP nTotalIoniz was set and reached.\n");
                return 1;
        }

        if( save.lgTraceConvergeBase )
        {
                if( iteration > 1 && save.lgTraceConvergeBaseHash )
                {
                        static int iter_punch=-1;
                        if( iteration !=iter_punch )
                                fprintf( save.ipTraceConvergeBase, "%s\n",save.chHashString );
                        iter_punch = iteration;
                }

                fprintf( save.ipTraceConvergeBase, 
                "%li\t%.4e\t%.4e\t%.4e\n",
                nzone , thermal.htot , thermal.ctot , dense.eden  );
        }

        return 0;
}

void UpdateUTAs( void )
{
        DEBUG_ENTRY( "UpdateUTAs()" );

        /* only reevaluate this on first pass through on this zone */
        if( conv.lgFirstSweepThisZone )
        {
                /* inner shell ionization */
                for( long nelem=0; nelem< LIMELM; ++nelem )
                {
                        for( long ion=0; ion<nelem+1; ++ion )
                        {
                                ionbal.UTA_ionize_rate[nelem][ion] = 0.;
                                ionbal.UTA_heat_rate[nelem][ion] = 0.;
                        }
                }
                /* inner shell ionization by UTA lines */
                /* this flag turned off with no UTA command */
                if( ionbal.lgInnerShellLine_on )
                {
                        for( long i=0; i < nUTA; ++i )
                        {
                                /* rateone is inverse lifetime of level against autoionization */
                                double rateone = UTALines[i].Emis().pump() * UTALines[i].Emis().AutoIonizFrac();
                                ionbal.UTA_ionize_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone;
                                /* heating rate in erg atom-1 s-1 */
                                ionbal.UTA_heat_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone*UTALines[i].Coll().heat();
                                {
                                        /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
                                        /*@-redef@*/
                                        enum {DEBUG_LOC=false};
                                        /*@+redef@*/
                                        if( DEBUG_LOC /*&& UTALines[i].nelem==ipIRON+1 && (UTALines[i].IonStg==15||UTALines[i].IonStg==14)*/ )
                                        {
                                                fprintf(ioQQQ,"DEBUG UTA %3i %3i %.3f %.2e %.2e %.2e\n",
                                                        (*UTALines[i].Hi()).nelem() , (*UTALines[i].Hi()).IonStg() , UTALines[i].WLAng() ,
                                                        rateone, UTALines[i].Coll().heat(), 
                                                        UTALines[i].Coll().heat()*dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] );
                                        }
                                }
                        }
                }
        }

        return;
}

STATIC bool lgNetEdenSrcSmall( void )
{
        fixit(); // this routine needs to be enabled
        return true;
        DEBUG_ENTRY( "lgNetEdenSrcSmall()" );

        if( conv.lgSearch )
                return true;
        fixit(); // grain rates not well tested below
        if( gv.lgDustOn() ) 
                return true;

        // Check for consistency of explicit source and sink rates for
        // electrons with population derived from neutrality
        double ionsrc = 0., ionsnk = 0.;
        for( long nelem=0; nelem < LIMELM; ++nelem )
        {
                if( !dense.lgElmtOn[nelem] )
                        continue;
                ionsrc += ionbal.elecsrc[nelem];
                ionsnk += ionbal.elecsnk[nelem];
                for( long ion_from = 0; ion_from <= nelem + 1; ++ion_from )
                {
                        for( long ion_to = 0; ion_to <= nelem + 1; ++ion_to )
                        {
                                if( ion_to-ion_from > 0 )
                                {
                                        ionsrc += gv.GrainChTrRate[nelem][ion_from][ion_to] *
                                                dense.xIonDense[nelem][ion_from] * (ion_to-ion_from);
                                }
                                else if( ion_to-ion_from < 0 )
                                {
                                        ionsnk += gv.GrainChTrRate[nelem][ion_from][ion_to] *
                                                dense.xIonDense[nelem][ion_from] * (ion_from-ion_to);
                                }
                        }
                }
        }
        long ipMElec = findspecies("e-")->index;
        const double totsrc = ionsrc + mole.species[ipMElec].src;
        const double totsnk = ionsnk + mole.species[ipMElec].snk*dense.EdenTrue;
        const double diff = (totsrc - totsnk);
        const double ave = ( fabs(totsrc) + fabs(totsnk) )/2.;
        fixit(); // Need to tighten up e- population convergence criterion
        const double error_allowed = 0.05 * ave; //conv.EdenErrorAllowed * ave;
        if( fabs(diff) > error_allowed )
        {
                enum {DEBUG_LOC=false};
                if( DEBUG_LOC )
                {
                        fprintf(ioQQQ,"PROBLEM large NetEdenSrc nzone %li\t%e\t%e\t%e\t%e\n",
                                nzone, 
                                totsrc/SDIV(totsnk),
                                dense.EdenTrue,
                                ionsrc + mole.species[ipMElec].src,
                                ionsnk + mole.species[ipMElec].snk*dense.EdenTrue);
                }
                conv.setConvIonizFail( "NetEdenSrc", diff, error_allowed);
                return false;
        }
        else
                return true;
}