/home66/gary/public_html/cloudy/c08_branch/source/conv_init_solution.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 */
/*ConvInitSolution drive search for initial temperature, for illuminated face */
/*lgTempTooHigh returns true if temperature is too high */
/*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */
#include "cddefines.h"
#include "physconst.h"
#include "trace.h"
#include "struc.h"
#include "rfield.h"
#include "mole.h"
#include "dense.h"
#include "stopcalc.h"
#include "heavy.h"
#include "wind.h"
#include "geometry.h"
#include "thermal.h"
#include "radius.h"
#include "phycon.h"
#include "pressure.h"
#include "conv.h"

/* derivative of net cooling wrt temperature to check on sign oscillations */
static double dCoolNetDTOld = 0;

static double OxyInGrains , FracMoleMax;
/* determine whether chemistry is important - what is the largest
 * fraction of an atom in molecules?  Also is ice formation 
 * important
 * sets above two file static variables */

/*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */
STATIC bool lgCoolHeatCheckConverge( 
                double *CoolNet )
{
        static double HeatOld=-1. , CoolOld=-1.;
        DEBUG_ENTRY( "lgCoolHeatCheckConverge()" );

        double HeatChange = thermal.htot - HeatOld;
        double CoolChange = thermal.ctot - CoolOld;
        /* will use larger of heating or cooling in tests for relative convergence
         * because in constant temperature models one may have trivially small value */
        double HeatCoolMax = MAX2( thermal.htot , thermal.ctot );

        /* is the heating / cooling converged?
         * get max relative change, determines whether converged heating/cooling */
        double HeatRelChange = fabs(HeatChange)/SDIV( HeatCoolMax );
        double CoolRelChange = fabs(CoolChange)/SDIV( HeatCoolMax );
        bool lgConverged = true;
        if( MAX2( HeatRelChange , CoolRelChange ) > conv.HeatCoolRelErrorAllowed )
                lgConverged = false;

        *CoolNet = thermal.ctot - thermal.htot;

        HeatOld = thermal.htot;
        CoolOld = thermal.ctot;
        return lgConverged;
}

/* call lgCoolHeatCheckConverge until cooling and heating are converged */
STATIC bool lgCoolNetConverge( 
                double *CoolNet , 
                double *dCoolNetDT )
{
        const long int LOOP_MAX=20;
        long int loop = 0;
        bool lgConvergeCoolHeat = false;

        DEBUG_ENTRY( "lgCoolNetConverge()" );

        while( loop < LOOP_MAX && !lgConvergeCoolHeat && !lgAbort )
        {
                if( ConvEdenIoniz() )
                {
                        /* this is an error return, calculation will immediately stop */
                        lgAbort = true;
                }
                lgConvergeCoolHeat = lgCoolHeatCheckConverge( CoolNet );
                if( trace.lgTrace || trace.nTrConvg>=3 )
                        fprintf(ioQQQ,"    lgCoolNetConverge %li calls to ConvEdenIoniz, converged? %c\n",
                        loop , TorF( lgConvergeCoolHeat ) );
                ++loop;
        }

        if( thermal.ConstTemp > 0 )
        {
                /* constant temperature model - this is trick so that temperature
                 * updater uses derivative to go to the set constant temperature */
                *CoolNet = phycon.te - thermal.ConstTemp;
                *dCoolNetDT = 1.f;
        }
        else if( phycon.te < 4000.f )
        {
                /* at low temperatures the analytical cooling-heating derivative
                 * is often of no value - the species densities change when the
                 * temperature changes.  Use simple dCnet/dT=(C-H)/T - this is 
                 * usually too large and results is smaller than optical dT, but
                 * we do eventually converge */
                *dCoolNetDT = thermal.ctot / phycon.te;
        }
        else if( thermal.htot / thermal.ctot > 2.f )
        {
                /* we are far from the solution and the temperature is much lower than
                 * equilibrium - analytical derivative from cooling evaluation is probably 
                 * wrong since coolants are not correct */
                *dCoolNetDT = thermal.ctot / phycon.te;
        }
        else 
        {
                *dCoolNetDT = thermal.dCooldT - thermal.dHeatdT;
                if( dCoolNetDTOld * *dCoolNetDT < 0. )
                {
                        /* derivative is oscillating */
                        *dCoolNetDT = *CoolNet / phycon.te;
                        fprintf(ioQQQ,"PROBLEM CoolNet derivative oscillating in initial solution \n");
                }
        }
        return lgConvergeCoolHeat;
}

/*ChemImportance find fraction of heavy elements in molecules, O in ices */
STATIC void ChemImportance( void )
{
        long int i;

        DEBUG_ENTRY( "ChemImportance()" );

        FracMoleMax = 0.;
        long int ipMax = -1;
        for( i=0; i<mole.num_comole_calc; ++i )
        {
                /* is element included in the chemistry, and is it a molecule? */
                if( dense.lgElmtOn[COmole[i]->nelem_hevmol] && COmole[i]->n_nuclei>1 )
                {
                        realnum frac = COmole[i]->hevmol / SDIV(dense.gas_phase[COmole[i]->nelem_hevmol]);
                        frac *= (realnum) COmole[i]->nElem[COmole[i]->nelem_hevmol];
                        if( frac > FracMoleMax )
                        {
                                /* remember largest fraction */
                                FracMoleMax = frac;
                                ipMax = i;
                        }
                }
        }

        /* fraction of all oxygen in ices on grains */
        OxyInGrains = 0;
        for(i=0;i<mole.num_comole_calc;++i)
        {
                OxyInGrains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN]; 
        }
        /* this is now fraction of O in ices */
        OxyInGrains /= SDIV(dense.gas_phase[ipOXYGEN]);

        {
                /* option to print out entire matrix */
                enum {DEBUG_LOC=false};
                if( DEBUG_LOC )
                {
                        fprintf(ioQQQ,"DEBUG fraction molecular %.2e species %li O ices %.2e\n" , 
                                FracMoleMax , ipMax ,OxyInGrains );
                }
        }

        return;
}

double FindTempChangeFactor(void)
{
        double TeChangeFactor;

        DEBUG_ENTRY( "FindTempChangeFactor()" );

        /* find fraction of atoms in molecules - need small changes
         * in temperature if fully molecular since chemistry
         * network is complex and large changes in solution would
         * cause linearization to fail */
        /* this evaluates the global variables FracMoleMax and
         * OxyInGrains */
        ChemImportance();

        /* Following are series of chemical 
         * tests that determine the factor by which
         * which can change the temperature.  must be VERY small when 
         * ice formation on grains is important due to exponential sublimation
         * rate.  Also small when molecules are important, to keep chemistry
         * within limits of linearized solver 
         * the complete linearization that is implicit in all these solvers
         * will not work when large Te jumps occur when molecules/ices
         * are important - solvers can't return to solution if too far away
         * keep temp steps small when mole/ice is important */
        if( OxyInGrains > 0.99  )
        {
                TeChangeFactor = 0.999;
        }
        else if( OxyInGrains > 0.9  )
        {
                TeChangeFactor = 0.99;
        }
        /* >>chng 97 mar 3, make TeChangeFactor small when close to lower bound */
        else if( phycon.te < 5. ||
                /*>>chng 06 jul 30, or if molecules/ices are important */
                OxyInGrains > 0.1  )
        {
                TeChangeFactor = 0.97;
        }
        /*>>chng 07 feb 23, add test of chemistry being very important */
        else if( phycon.te < 20. || FracMoleMax > 0.1 )
        {
                /* >>chng 07 feb 24, changed from 0,9 to 0.95 to converge 
                 * pdr_coolbb.in test */
                TeChangeFactor = 0.95;
        }
        else
        {
                /* this is the default largest step */
                TeChangeFactor = 0.8;
        }
        return TeChangeFactor;
}

/*ConvInitSolution drive search for initial temperature, for illuminated face */
int ConvInitSolution(void)
{
        long int i, 
          ionlup, 
          nelem ,
          nflux_old,
          nelem_reflux ,
          ion_reflux;
        /* current value of Cooling - Heating */
        bool lgConvergeCoolHeat;

        double 
          TeChangeFactor, 
          TeBoundLow, 
          TeBoundHigh;

        DEBUG_ENTRY( "ConvInitSolution()" );

        /* set NaN for safety */
        set_NaN( TeBoundLow );
        set_NaN( TeBoundHigh );

        /* this counts number of times ConvBase is called by PressureChange, in current zone */
        conv.nPres2Ioniz = 0;
        /* this counts how many times ConvBase is called in this iteration
         * and is flag used by various routines to understand they are 
         * being called the first time*/
        conv.nTotalIoniz = 0;
        /* ots rates not oscillating */
        conv.lgOscilOTS = false;

        lgAbort = false;
        dense.lgEdenBad = false;
        dense.nzEdenBad = 0;
        /* these are variables to remember the biggest error in the
         * electron density, and the zone in which it occurred */
        conv.BigEdenError = 0.;
        conv.AverEdenError = 0.;
        conv.BigHeatCoolError = 0.;
        conv.AverHeatCoolError = 0.;
        conv.BigPressError = 0.;
        conv.AverPressError = 0.;

        /* these are equal if set dr was used, and we know the first dr */
        if( fp_equal( radius.sdrmin, radius.sdrmax ) )
        {
                radius.drad = MIN2( radius.sdrmax  ,  radius.drad );
                radius.drad_x_fillfac = radius.drad * geometry.FillFac;
        }

        /* the H+ density is zero in sims with no H-ionizing radiation and no cosmic rays,
         * the code does work in this limit.  But it is unphysical for the H0 density
         * to be zero.  This could only happen in the fully molecular limit, and the
         * code does not go to this limit at this stage, due to simple approximations. */
        ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);

        if( trace.nTrConvg )
        {
                fprintf( ioQQQ, "\nConvInitSolution entered \n" );
        }

        /********************************************************************
         *
         * this is second or higher iteration, reestablish original temperature
         *
         *********************************************************************/
        if( iteration != 1 )
        {
                /* this is second or higher iteration on multi-iteration model */
                if( trace.lgTrace || trace.nTrConvg )
                {
                        fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n", 
                                iteration, struc.testr[0] );
                }

                if( trace.lgTrace || trace.nTrConvg )
                {
                        fprintf( ioQQQ, " search set true\n" );
                }

                /* search phase must be turned on so that variables such as the ots rates,
                 * secondary ionization, and auger yields, can converge more quickly to 
                 * proper values */
                conv.lgSearch = true;

                /* this is the temperature and pressure from the previous iteration */
                /*phycon.te = TempPrevIteration;*/
                /*>>chng 06 jun 09, reset to this rather than init value set in this routine */
                TempChange( struc.testr[0] , false);

                /* the initial pressure should now be valid */
                /* this sets values of pressure.PresTotlCurr */
                PresTotCurrent();

                /* now get final temperature */
                if( ConvPresTempEdenIoniz() )
                {
                        /* this is an error return, calculation will immediately stop */
                        lgAbort = true;
                        return(1);
                }

                if( trace.lgTrace || trace.nTrConvg )
                {
                        fprintf( ioQQQ, " ========================================================================================================================\n");
                        fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te);
                        fprintf( ioQQQ, " ========================================================================================================================\n");
                }
                conv.lgSearch = false;

                /* now get final temperature */
                if( ConvTempEdenIoniz() )
                {
                        /* this is an error return, calculation will immediately stop */
                        lgAbort = true;
                        return(1);
                }

                /* the initial pressure should now be valid 
                 * this sets values of pressure.PresTotlCurr */
                PresTotCurrent();
        }

        else
        {
                /********************************************************************
                 *
                 * do first te from scratch 
                 *
                 *********************************************************************/

                /* say that we are in search phase */
                conv.lgSearch = true;

                if( trace.lgTrace )
                {
                        fprintf( ioQQQ, " ConvInitSolution called, new temperature.\n" );
                }

                /* coming up to here Te is either 4,000 K (usually) or 10^6 */
                TeBoundLow = phycon.TEMP_LIMIT_LOW/1.001;

                /* set temperature floor option  - StopCalc.TeFloor is usually
                * zero but reset this this command - will go over to constant 
                * temperature calculation if temperature falls below floor */
                TeBoundLow = MAX2( TeBoundLow , StopCalc.TeFloor );

                /* highest possible temperature - must not get this high since
                 * parts of code will abort if value is reached. 
                 * divide by 1.2 to prevent checking on temp > 1e10 */
                TeBoundHigh = phycon.TEMP_LIMIT_HIGH/1.2;

                /* set initial temperature, options are constant temperature,
                 * or approach equilibrium from either high or low TE */
                double TeFirst;
                if( thermal.ConstTemp > 0 )
                {
                        /* constant temperature , set 4000 K then relax to desired value
                         * for cold temperatures, to slowly approach fully molecular solution.
                         * if constant temperature is higher than 4000., go right to it */
                        TeFirst = MAX2( 4000. , thermal.ConstTemp );
                        /* true allow ionization stage trimming, false blocks it */
                }
                else if( thermal.lgTeHigh )
                {
                        /* approach from high TE */
                        TeFirst = 1e6;
                }
                else
                {
                        /* approach from low TE */
                        TeFirst = 4000.;
                }
                TempChange(TeFirst , false);
                if( trace.lgTrace || trace.nTrConvg>=2 )
                        fprintf(ioQQQ,"ConvInitSolution doing initial solution with temp=%.2e\n", 
                        phycon.te);

                /* this sets values of pressure.PresTotlCurr */
                PresTotCurrent();

                thermal.htot = 1.;
                thermal.ctot = 1.;

                /* call cooling, heating, opacity, loop to convergence
                 * this is very first call to it, by default is at 4000K */

                double CoolNet=0, dCoolNetDT=0;
                /* do ionization twice to get stable solution, evaluating initial dR each time */
                lgConvergeCoolHeat = false;
                for( ionlup=0; ionlup<2; ++ionlup )
                {
                        if( trace.lgTrace || trace.nTrConvg>=2 )
                                fprintf( ioQQQ, "ConvInitSolution calling lgCoolNetConverge "
                                "initial setup %li with Te=%.3e C=%.2e H=%.2e d(C-H)/dT %.3e\n", 
                                ionlup,phycon.te  , thermal.ctot , thermal.htot,dCoolNetDT );
                        /* do not flag oscillating d(C-H)/dT until stable */
                        dCoolNetDTOld = 0.;
                        lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , &dCoolNetDT );
                        if( lgAbort )
                                break;
                        /* set thickness of first zone, this affects the pumping rates due
                         * to correction for attenuation across zone, so must be known 
                         * for ConvEdenIoniz to get valid solution - this is why it
                         * is in this loop */
                        radius_first();
                }
                if( !lgConvergeCoolHeat )
                        fprintf(ioQQQ," PROBLEM ConvInitSolution did not converge the heating-cooling during initial search.\n");

                if( lgAbort )
                {
                        /* we hit an abort */
                        fprintf( ioQQQ, " DISASTER PROBLEM ConvInitSolution: Search for "
                                "initial conditions aborted - lgAbort set true.\n" );
                        ShowMe();
                        /* this is an error return, calculation will immediately stop */
                        return(1);
                }

                /* we now have error in C-H and its derivative - following is better
                 * derivative for global case where we may be quite far from C==H */
                if( thermal.ConstTemp<=0 )
                        dCoolNetDT = thermal.ctot / phycon.te;
                bool lgConvergedLoop = false;
                long int LoopThermal = 0;
                /* initial temperature is usually 4000K, if the dTe scale factor is 0.95
                 * it will take 140 integrations to lower temperature to 3K,
                 * and many more if ices are important.  */
                const long int LIMIT_THERMAL_LOOP=300;
                double CoolMHeatSave=-1. , TempSave=-1. , TeNew=-1.,CoolSave=-1;
                while( !lgConvergedLoop && LoopThermal < LIMIT_THERMAL_LOOP )
                {
                        /* change temperature until we sign of C-H changes */
                        CoolMHeatSave = CoolNet;
                        dCoolNetDTOld = dCoolNetDT;
                        CoolSave = thermal.ctot;
                        TempSave = phycon.te;

                        /* find temperature change factor, a number less than 1*/
                        TeChangeFactor = FindTempChangeFactor();
                        ASSERT( TeChangeFactor>0. && TeChangeFactor< 1. );

                        TeNew = phycon.te - CoolNet / dCoolNetDT;

                        TeNew = MAX2( phycon.te*TeChangeFactor , TeNew );
                        TeNew = MIN2( phycon.te/TeChangeFactor , TeNew );
                        /* change temperature */
                        TempChange(TeNew , true);
                        lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , & dCoolNetDT );

                        if( trace.lgTrace || trace.nTrConvg>=2 )
                                fprintf(ioQQQ,"ConvInitSolution %4li TeNnew=%.3e "
                                "TeChangeFactor=%.3e Cnet/H=%.2e dCnet/dT=%10.2e Ne=%.3e"
                                " Ograins %.2e FracMoleMax %.2e\n",
                                LoopThermal , TeNew , TeChangeFactor ,
                                CoolNet/SDIV(thermal.htot) ,  dCoolNetDT,
                                dense.eden , OxyInGrains , FracMoleMax );
                        if( lgAbort )
                        {
                                /* we hit an abort */
                                fprintf( ioQQQ, " Search for initial conditions aborted - "
                                        "lgAbort set true, lgConvergeCoolHeat=%c.\n",
                                        TorF(lgConvergeCoolHeat) );
                                /* this is an error return, calculation will immediately stop */
                                return(1);
                        }

                        /* keep changing temperature until sign of CoolNet changes
                         * in constant temperature case CoolNet is delta Temperature
                         * so is zero for last pass in this loop 
                         * this is for constant temperature case */
                        if( fabs(CoolNet)< SMALLDOUBLE )
                                /* CoolNet is zero */
                                lgConvergedLoop = true;
                        else if( (CoolNet/fabs(CoolNet))*(CoolMHeatSave/fabs(CoolMHeatSave)) <= 0)
                                /* change in sign of net cooling */
                                lgConvergedLoop = true;
                        else if( fabs(CoolNet)/MAX2( thermal.ctot , thermal.htot ) <conv.HeatCoolRelErrorAllowed*10. )
                                lgConvergedLoop = true;
                        else if( phycon.te <= TeBoundLow || phycon.te>=TeBoundHigh)
                                lgConvergedLoop = true;
                        ++LoopThermal;
                }

                if( !lgConvergedLoop )
                        fprintf(ioQQQ,"PROBLEM ConvInitSolution: temperature solution not "
                        "found in initial search, final Te=%.2e\n",
                        phycon.te );

                /* interpolate temperature where C=H if not constant temperature sim 
                 * will have set constant temperature mode above if we hit temperature
                 * floor */
                if( thermal.ConstTemp<=0 && 
                        ! fp_equal( TempSave , phycon.te ) )
                {
                        if( trace.lgTrace || trace.nTrConvg>=2 )
                                fprintf(ioQQQ," Change in sign of C-H found, Te brackets %.3e "
                                "%.3e Cool brackets %.3e %.3e ",
                                TempSave , phycon.te , CoolMHeatSave , CoolNet);
                        /* interpolate new temperature assuming Cool = a T^alpha,
                         * first find alpha */
                        double alpha = log(thermal.ctot/CoolSave) / log( phycon.te/TempSave);
                        if( fabs(alpha) < SMALLFLOAT )
                                /* alpha close to 0 means constant temperature */
                                TeNew = phycon.te;
                        else
                        {
                                /* next find log a - work in logs to avoid infinities */
                                if( thermal.ctot<0. || thermal.htot<0. )
                                        TotalInsanity();
                                double Alog = log( thermal.ctot ) - alpha * log( phycon.te );
                                /* the interpolated temperature where heating equals cooling */
                                double TeLn = (log( thermal.htot ) - Alog) / alpha ;

                                /* a sanity check */
                                if( TeLn < log( MIN2(phycon.te , TempSave )) )
                                        TeNew = MIN2(phycon.te , TempSave );
                                else if( TeLn > log( MAX2(phycon.te , TempSave )) )
                                        TeNew = MAX2(phycon.te , TempSave );
                                else
                                        TeNew = pow( EE , TeLn );
                        }

                        ASSERT( TeNew>=MIN2 ( TempSave , phycon.te ) ||
                                TeNew<=MAX2 ( TempSave , phycon.te ));

                        if( trace.lgTrace || trace.nTrConvg>=2 )
                                fprintf(ioQQQ," interpolating to Te %.3e \n\n",
                                TeNew);

                        /* change temperature */
                        TempChange(TeNew , true);
                }

                if( ConvTempEdenIoniz() )
                {
                        /* this is an error return, calculation will immediately stop */
                        lgAbort = true;
                        return(1);
                }

                /* this sets values of pressure.PresTotlCurr */
                PresTotCurrent();

                if( trace.lgTrace || trace.nTrConvg )
                {
                        fprintf( ioQQQ, "\nConvInitSolution: end 1st iteration search phase "
                                "finding Te=%.3e\n\n" , phycon.te);
                }
                conv.lgSearch = false;

                if( trace.lgTrace )
                {
                        fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n", 
                          phycon.te );
                }
        }

        /* we now have a fairly good temperature and ionization 
         * iteration is 1 on first iteration - remember current pressure
         * on first iteration so we can hold this constant in constant 
         * pressure simulation.
         *
         * flag dense.lgDenseInitConstant false if
         * *constant pressure reset* is used - 
         * default is true, after first iteration initial density is used for
         * first zone no matter what pressure results.  Small changes occur due
         * to radiative transfer convergence, 
         * when set false with *reset* option the density on later iterations 
         * can change to keep the pressure constant */
        if( iteration==1 || dense.lgDenseInitConstant )
        {
                double PresNew = pressure.PresTotlCurr;

                /* option to specify the pressure as option on constant pressure command */
                if( pressure.lgPressureInitialSpecified )
                        /* this is log of nT product - if not present then set zero */
                        PresNew = pressure.PressureInitialSpecified;
                if( trace.lgTrace )
                {
                        fprintf( ioQQQ, 
                                "     PresTotCurrent 1st zn old values of PresTotlInit:%.3e"
                                " to %.3e \n",
                                pressure.PresTotlInit,
                                PresNew);
                }

                pressure.PresTotlInit = PresNew;
        }

        /* now bring current pressure to the correct pressure - must do this
         * before initial values are saved in iter start/end */
        ConvPresTempEdenIoniz();

        /* this counts number of times ConvBase is called by PressureChange, in
         * current zone these are reset here, so that we count from first zone 
         * not search */
        conv.nPres2Ioniz = 0;

        dense.lgEdenBad = false;
        dense.nzEdenBad = 0;

        /* save counter upon exit so we can see how many iterations were
         * needed to do true zones */
        conv.nTotalIoniz_start = conv.nTotalIoniz;

        /* these are variables to remember the biggest error in the
         * electron density, and the zone in which it occurred */
        conv.BigEdenError = 0.;
        conv.AverEdenError = 0.;
        conv.BigHeatCoolError = 0.;
        conv.AverHeatCoolError = 0.;
        conv.BigPressError = 0.;
        conv.AverPressError = 0.;

        /*>>chng 06 jun 09, only reset on first try - otherwise have stable solution? */
        if( iteration == 1 )
        {
                /* now remember some things we may need even in first zone, these
                * are normally set towards end of zone calculation in RT_tau_inc */
                struc.testr[0] = (realnum)phycon.te;
                /* >>chng 02 May 2001 rjrw: add hden for dilution */
                struc.hden[0] = dense.gas_phase[ipHYDROGEN];
                /* pden is the total number of particles per unit vol */
                struc.DenParticles[0] = dense.pden;
                struc.heatstr[0] = thermal.htot;
                struc.coolstr[0] = thermal.ctot;
                struc.volstr[0] = (realnum)radius.dVeff;
                struc.drad_x_fillfac[0] = (realnum)radius.drad_x_fillfac;
                struc.histr[0] = dense.xIonDense[ipHYDROGEN][0];
                struc.hiistr[0] = dense.xIonDense[ipHYDROGEN][1];
                struc.ednstr[0] = (realnum)dense.eden;
                struc.o3str[0] = dense.xIonDense[ipOXYGEN][2];
                struc.DenMass[0] = dense.xMassDensity;
                struc.drad[0] = (realnum)radius.drad;
        }

        /* check that nflux extends above IP of highest ionization species present.
         * for collisional case it is possible for species to exist that are higher
         * IP than the limit to the continuum.  Need continuum to encompass all
         * possible emission - to account for diffuse emission
         * NB 
         * on second iteration of multi-iteration model this may result in rfield.nflux increasing
         * which can change the solution */
        nflux_old = rfield.nflux;
        nelem_reflux = -1;
        ion_reflux = -1;
        for( nelem=2; nelem < LIMELM; nelem++ )
        {
                /* do not include hydrogenic species in following */
                for( i=0; i < nelem; i++ )
                {
                        if( dense.xIonDense[nelem][i+1] > 0. )
                        {
                                if( Heavy.ipHeavy[nelem][i] > rfield.nflux )
                                {
                                        rfield.nflux = Heavy.ipHeavy[nelem][i];
                                        nelem_reflux = nelem;
                                        ion_reflux = i;
                                }
                        }
                }
        }

        /* was the upper limit to the continuum updated? if so need to define
         * continuum variables to this new range */
        if( nflux_old != rfield.nflux )
        {
                /* zero out parts of rfield arrays that were previously undefined */
                rfield_opac_zero( nflux_old-1 , rfield.nflux );

                /* if continuum reset up, we need to define gaunt factors through high end */
                /*tfidle(false);*/
                /* this calls tfidle, among other things */
                /* this sets values of pressure.PresTotlCurr */
                PresTotCurrent();

                /* redo ionization and update opacities */
                if( ConvBase(1) )
                {
                        /* this is catastrophic failure */
                        lgAbort = true;
                        return( 1 );
                }

                /* need to update continuum opacities */
                if( trace.lgTrace )
                {
                        fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n",
                                nflux_old , rfield.nflux ,
                                rfield.anu[nflux_old] , rfield.anu[rfield.nflux] );
                        fprintf(ioQQQ," caused by element %li ion %li \n",
                          nelem_reflux ,ion_reflux );
                }
        }

        /* wind with negative velocity is flow from PDR - change density slightly
         * and call pressure solver to see if it returns to original density */
        if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. )
        {
                /* >>chng 04 apr 23, change pressure and let solver come back to correct
                 * pressure.  This trick makes sure we are correctly either super or sub sonic 
                 * since solver will jump from one branch to the other if necessary */
#               define PCHNG    0.98
                /* this ignores return condition, assume must be ok if got this far */
                if( ConvPresTempEdenIoniz() )
                {
                        /* this is an error return, calculation will immediately stop */
                        lgAbort = true;
                        return(1);
                }

                pressure.PresTotlInit *= PCHNG;
                if( ConvPresTempEdenIoniz() )
                {
                        /* this is an error return, calculation will immediately stop */
                        lgAbort = true;
                        return(1);
                }

                pressure.PresTotlInit /= PCHNG;
                if( ConvPresTempEdenIoniz() )
                {
                        /* this is an error return, calculation will immediately stop */
                        lgAbort = true;
                        return(1);
                }

#               undef PCHNG
        }
        /* this is the only valid return and lgAbort should be false*/
        return( lgAbort );
}

---