/home66/gary/public_html/cloudy/c08_branch/source/radius_next.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 */
/*radius_next use adaptive logic to find next zone thickness */
/*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
/*GrainRateDr called by radius_next to find grain heating rate dr */
#include "cddefines.h"
#include "lines_service.h"
#include "iso.h"
#include "geometry.h"
#include "h2.h"
#include "mole.h"
#include "hyperfine.h"
#include "opacity.h"
#include "dense.h"
#include "heavy.h"
#include "grainvar.h"
#include "elementnames.h"
#include "conv.h"
#include "rfield.h"
#include "dynamics.h"
#include "thermal.h"
#include "hmi.h"
#include "coolheavy.h"
#include "timesc.h"
#include "doppvel.h"
#include "stopcalc.h"
#include "colden.h"
#include "phycon.h"
#include "rt.h"
#include "trace.h"
#include "wind.h"
#include "punch.h"
#include "taulines.h"
#include "pressure.h"
#include "iterations.h"
#include "struc.h"
#include "radius.h"

#if 0
/*ChkRate called by radius_next to check how rates of destruction of various species changes */
STATIC void ChkRate(
          /* element number on physical scale */
          long int nelem, 
          /* change in destruction rate */
          double *dDestRate, 
          /* old and new destruction rates */
          double *DestRateOld,
          double *DestRateNew,
          /* stage of ionization on the physical scale */
          long int *istage);
#endif

/*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
STATIC void ContRate(double *freqm, 
  double *opacm);

/*GrainRateDr called by radius_next to find grain heating rate dr */
STATIC void GrainRateDr(double *grfreqm, 
  double *gropacm);

/*radius_next use adaptive logic to find next zone thickness 
 * return 0 if ok, 1 for abort */
int radius_next(void)
{
        char chLbl[11];
        bool lgDoPun;

        long int level , ipStrong;

        double thickness_total , drThickness , DepthToGo , AV_to_go;
        int mole_dr_change;

        double DrGrainHeat, 
          GlobDr, 
          SaveOHIIoHI, 
          SpecDr, 
          Strong, 
          TauDTau, 
          TauInwd, 
          drSolomon_BigH2 ,
          dEfrac, 
          dHdStep, 
          dRTaue, 
          dTdStep, 
          dnew, 
          drConPres, 
          drH2_heat_cool = 0. ,
          dH2_heat_cool = 0.,
          drH2_abund = 0. ,
          dr_mole_abund = 0.,
          dH2_abund=0.,
          dCO_abund=0.,
          drDepth, 
          drDest, 
          drEfrac, 
          drFail, 
          drFluc, 
          drHMase, 
          drHe1ovHe2,
          drHion, 
          drInter, 
          drLeiden_hack ,
          drLineHeat, 
          drSphere, 
          drTab, 
          drdHdStep, 
          drdTdStep, 
          drThermalFront ,
          drmax, 
          dVelRelative,
          fac2, 
          factor, 
          freqm, 
          grfreqm=0., 
          gropacm=0., 
          hdnew, 
          opacm, 
          recom_dr_last_iter ,
          OldDR ,
          winddr, 
          WindAccelDR,
          x;

        double change_heavy_frac=-1. , change_heavy_frac_big , dr_change_heavy ,
                frac_change_heavy_frac_big, Efrac_old , Efrac_new;
        long int ichange_heavy_nelem=-1 , nelem , ion , ichange_heavy_ion=-1;

        static double OHIIoHI, 
          OldHeat = 0., 
          OldTe = 0.,
          OlddTdStep = 0.,
          OldH2Abund=0.,
          OldWindVelocity=0.,
          Old_He_atom_ov_ion = 0,
          Old_H2_heat_cool;
        static long int iteration_last=-1;

        static double BigRadius = 1e30;
        bool lgFirstCall;

        DEBUG_ENTRY( "radius_next()" );


        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         *
         * changes in logic
         * 95 oct 19, drSphere now 3% of radius, was 2%, fewer zone
         * 95 oct 19, subtracted grain opacity from total opacity used
         * to get thickness in routine ContRate
         *
         *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         *
         * free statement labels >= 13
         *
         *-----------------------------------------------------------------------
         *
         * this sub determines the thickness of the next zone
         * if is called one time for each zone
         * flag lgNxtDROff is true if this is initialization of radius_next,
         * is false if we are to use logic to find dr
         *
         *----------------------------------------------------------------------- */

        /* >>chng 03 sep 21 - decide whether this is the first call */
        if( (iteration != iteration_last) && (nzone==0) )
        {
                /* this is the first call in this iteration */
                iteration_last = iteration;
                lgFirstCall = true;
        }
        else
                lgFirstCall = false;

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

        /* save current dr */
        OldDR = radius.drad;

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /*  1    '' radius_next keys from change in H ionization'',e11.3)')
         * check on change in hydrogen ionizaiton */
        if( lgFirstCall )
        {
                if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
                {
                        OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
                }
                else
                {
                        OHIIoHI = 0.;
                }
                SaveOHIIoHI = OHIIoHI;
                drHion = BigRadius;
                /* else if(hii.gt.0. .and. hi.gt.0. .and. OHIIoHI.gt.0. ) then
                 * >>chng 97 jul 9, for deep in PDR vastly now ionz H slowed down works */
        }

        else if( (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) && OHIIoHI > 1e-15 )
        {
                double atomic_frac = (dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]);
                /* ratio of current HII/HI to old value - < 1 when becoming more neutral */
                /* this is now change in atomic fraction */
                x = 1. - atomic_frac /OHIIoHI;
                if( atomic_frac > 0.05 && atomic_frac < 0.9 )
                {
                        /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
                         * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
                        /* >>chng 02 jun 05 from 0.2 to 0.05 poorly resolved i-front, also added two-branch logic*/
                        drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) );
                }
                else if( x > 0. )
                {
                        /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
                         * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
                        drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) );
                }
                else
                {
                        drHion = BigRadius;
                }
                SaveOHIIoHI = OHIIoHI;
                OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
        }

        else
        {
                SaveOHIIoHI = OHIIoHI;
                if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
                {
                        OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
                }
                else
                {
                        OHIIoHI = 0.;
                }
                drHion = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* "radius_next keys from H maser, dt, ij=" possible hydrogen maser action */
        if( rt.dTauMase < -0.01 )
        {
                /* maser so powerful, must limit inc in tay
                 * >>chng 96 aug 08, from 0.3 to 0.1 due to large maser crash */
                drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase);
        }
        else
        {
                drHMase = BigRadius;
        }

        /* >>chng 03 nov 09, try doing he in following, not above */
        Old_He_atom_ov_ion = 0.;
        drHe1ovHe2 = BigRadius;

        /* check on N0 - 1 ionization changes,
         * >>chng 03 jun 06, add this test due to smashing into H ifront in blr89.in
         */
        dr_change_heavy = BigRadius;
        change_heavy_frac_big = -1.;
        frac_change_heavy_frac_big = -1.;
        /* >chng 03 jun 09, from 0.05 to 0.1, initial tests with zoning */
#       define CHANGE_ION_HEAV  0.2f
#       define CHANGE_ION_HHE   0.15f
        if( nzone > 4 )
        {
                for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
                {
                        if( dense.lgElmtOn[nelem] )
                        {
                                realnum change;
                                /* this is the limit to the ionization we will check -
                                 * also used in prt_comment to check on whether oscillations occurred */
                                realnum frac_limit;
                                if( nelem<=ipHELIUM || nelem==ipCARBON )
                                {
                                        frac_limit = 1e-4f;
                                        change = CHANGE_ION_HHE;
                                }
                                else
                                {
                                        /* this var is used to print warnings,
                                         * make sure we converge below it */
                                        frac_limit = struc.dr_ionfrac_limit/2.f;
                                        change = CHANGE_ION_HEAV;
                                }
                                /* >>chng 03 dec 09, go up to full range of ion, not just =2 */
                                for( ion=0; ion<=nelem+1; ++ion )
                                {
                                        realnum abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
                                        realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
                                        if( abundnew > frac_limit && abundold > frac_limit )
                                        {
                                                /* NB must make sure this test is not done when nzone-x <0 */
                                                /* -2 because previous zone, and nzone is off by one (on physical, not C, scale) */
                                                /*realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);*/
                                                realnum abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]);
                                                realnum abundoldest = struc.xIonDense[nelem][ion][nzone-5]/SDIV( struc.gas_phase[nelem][nzone-5]);
                                                /* this is fractional change */
                                                /* >>chng 04 feb 28, take min of old and new abund, to try to pick up
                                                 * rapid changing Ca+ sooner */
                                                change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew);
                                                /* want fractional change to be less than this factor */
                                                if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
                                                        /* >>chng 03 dec 07, add test that abund is not oscillating */
                                                        /* also test that abundance is increasing - we are headed into a front */
                                                        ((abundnew-abundold)>0.)   && 
                                                        ((abundold-abundolder)>0.) && 
                                                        ((abundolder-abundoldest)>0.) )
                                                {
                                                        ichange_heavy_nelem = nelem;
                                                        ichange_heavy_ion = ion;
                                                        change_heavy_frac_big = change_heavy_frac;
                                                        frac_change_heavy_frac_big = abundnew;
                                                        /* >>chng 03 dec 06, from min of 0.5 to min of 0.25, crash into He i-front 
                                                         * in hizqso.in */
                                                        /* >>chng 04 mar 03, min had become 0.1, forced oscillations in nova.in
                                                         * in silicon, zoning changed greatly, causing change in diffuse lin
                                                         * pumping.  put back to 0.25 */
                                                        dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac );
                                                }
                                        }
                                }
                        }
                }
        }

        /* if Leiden hacks are on then use increase in dust extinction as control 
         * >>chng 05 aug 13, add this */
        if(!co.lgUMISTrates)
        {
                /* Draine field is only defined over narrow range in FUV - must not let change
                 * in extinction become too large - 
                 * prefactor is change in optical depth */
                drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point) / SDIV( geometry.FillFac * rfield.opac_mag_V_point );
        }
        else
        {
                drLeiden_hack = BigRadius;
        }
        /* >>chng 04 feb 15, kill this block - not used */

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* check how heating is changing
         * '' radius_next keys from change in heating; current, delta='', */
        if( nzone <= 1 || thermal.lgTemperatureConstant )
        {
                drdHdStep = BigRadius;
                dHdStep = FLT_MAX;
        }
        else
        {
                dHdStep = fabs(thermal.htot-OldHeat)/thermal.htot;
                if( dHdStep > 0. )
                {
                        if( dense.gas_phase[ipHYDROGEN] >= 1e13 )
                        {
                                /* drdHdStep = drad * MAX( 0.8 , 0.05/dHdStep ) */
                                drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep);
                        }
                        else if( dense.gas_phase[ipHYDROGEN] >= 1e11 )
                        {
                                /* drdHdStep = drad * MAX( 0.8 , 0.075/dHdStep ) */
                                drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep);
                        }
                        else
                        {
                                /* changed from .15 to .12 for outer edge of coolhii too steep dT
                                 * changed to .10 for symmetry, big change in some rates, 95nov14
                                 * changed from .10 to .125 since parispn seemed to waste zones
                                 * >>chng 96 may 21, from .125 to .15 since pn's still waste zones */
                                drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep);
                        }
                }
                else
                {
                        drdHdStep = BigRadius;
                }
        }
        OldHeat = thermal.htot;

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* pressure due to incident continuum if in eos */
        if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn )
        {
                if( nzone > 2 && pressure.pinzon > 0. )
                {
                        /* pinzon is pressrue from acceleration onto previos zone
                         * want this less than some fraction of total pressure */
                        /* >>chng 06 feb 01, change from init pres to current total pressure
                         * in const press high U ulirgs current pressure may be quite larger
                         * than init pressure due to continuum absorption */
                        drConPres = 0.05*pressure.PresTotlCurr/(wind.AccelTot*
                          dense.xMassDensity*geometry.FillFac);
                }
                else
                {
                        drConPres = BigRadius;
                }
        }
        else
        {
                drConPres = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* check how temperature is changing
         *  1    '' radius_next keys from change in temperature; current, delta='', */
        if( nzone <= 1 )
        {
                drdTdStep = BigRadius;
                dTdStep = FLT_MAX;
                OlddTdStep = dTdStep;
        }
        else
        {
                /* change in temperature; current=      */
                dTdStep = (phycon.te-OldTe)/phycon.te;
                /* >>chng 02 dec 08, desired change in temperature per zone must not
                 * be smaller than allower error in temperature.  For now use relative error
                 * in heating - - cooling balance.  Better would be to also use c-h deriv wrt t
                 * to get slope */
                x = conv.HeatCoolRelErrorAllowed*2.;
                x = MAX2( 0.01 , x ); 
                x = MIN2( 0.05 , x );
                /* >>chng 02 dec 11 rjrw change back to 0.03 -- improve dynamics.dRad criterion instead */
                x = 0.03;
                /* >>chng 02 dec 07, do not do this if there is mild te jitter, 
                 * so that dT is changing sign - this happens
                 * in ism.ini, very stable temperature with slight noise up and down */
                if( dTdStep*OlddTdStep > 0. )
                {
                        /* >>chng 96 may 30, variable depending on temp
                         * >>chng 96 may 31, allow 0.7 smaller, was 0.8
                         * >>chng 97 may 05, from 0.7 to 0.5 stop from punching through thermal front
                         * >>chng 04 mar 30, from 0.7 to 0.5 stop from punching through thermal front,
                         * for some reason factor was 0.7, not 0.5 as per previous change */
                        double absdt = fabs(dTdStep);
                        drdTdStep = radius.drad*MAX2(0.5,x/absdt);
                }
                else
                {
                        drdTdStep = BigRadius;
                }
        }
        OlddTdStep = dTdStep;
        OldTe = phycon.te;

        /* >>chng 02 oct 06, only check on opacity - interaction if not
         * constant temperature - there were constant temperature models that
         * extended to infinity but hung with last few photons and this test.
         * better to ignore this test which is really for thermal models */
        /* >>chng 06 mar 20, do not call if recombination logic in place */
        if( !thermal.lgTemperatureConstant && !dynamics.lgRecom )
        {
                /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
                /* find freq where opacity largest and interaction rate is fastest
                * "cont inter nu=%10.2e opac=%10.2e\n" */
                ContRate(&freqm,&opacm);

                /* use optical depth at max interaction energy
                * >>chng 96 jun 06 was drChange=0.15 changed to 0.3 for high Z models
                * taking too many zones
                * drInter = drChange / MAX(1e-30,opacm*FillFac) */

                drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin);
        }
        else
        {
                drInter = BigRadius;
                freqm = 0.;
                opacm = 1.;
        }

#       if 0
        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* check on changes in destruction rates for various atoms */
        ChkRate(ipCARBON,&dDRCarbon,&DestOldCarb, &DestNewCarb,&icarstag);
        ChkRate(ipNITROGEN,&dDRNitrogen,&DestOldNit, &DestNewNit,&initstag);
        ChkRate(ipOXYGEN,&dDROxygen,&DestOldOxy, &DestNewOxy,&ioxystag);
        ChkRate(ipIRON,&dDRIron,&DestOldIron, &DestNewIron,&iironstag);

        dDestRate = MAX4(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen);

        if( fp_equal( dDRCarbon, dDestRate ) )
        {
                dDestRate = dDRCarbon;
                DestOld = DestOldCarb;
                DestNew = DestNewCarb;
                istage = icarstag;
                strcpy( chDestAtom, "Carbon  " );
        }

        else if( fp_equal( dDRNitrogen, dDestRate ) )
        {
                dDestRate = dDRNitrogen;
                DestOld = DestOldNit;
                DestNew = DestNewNit;
                istage = initstag;
                strcpy( chDestAtom, "Nitrogen" );
        }

        else if( fp_equal( dDROxygen, dDestRate ) )
        {
                dDestRate = dDROxygen;
                DestOld = DestOldOxy;
                DestNew = DestNewOxy;
                strcpy( chDestAtom, "Oxygen  " );
                istage = ioxystag;
        }

        else if( fp_equal( dDRIron, dDestRate ) )
        {
                dDestRate = dDRIron;
                DestOld = DestOldIron;
                DestNew = DestNewIron;
                istage = iironstag;
                strcpy( chDestAtom, "Iron    " );
        }

        else
        {
                fprintf( ioQQQ, " insanity following ChkRate\n" );
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }

        /*  radius_next keys from change in dest rates, atom= */
        if( dDestRate > 0. )
        {
                /* if( te.gt.40 000. ) then
                 * added different accuracy for hot gas since tend to jump over
                 * big te range for small change in heat (intrinsically unstable)
                 * drDest = drad * MAX( 0.5 , 0.10/dDestRate )
                 * else
                 * was drChange, changed to .15 for parishii go through HeII=HeI I front
                 * drDest = drad * MAX( 0.5 , 0.15/dDestRate )
                 * >>chng 95 dec 18 from above to below to stop oscillations
                 * >>chng 95 dec 27 from min of .5 to .75 to stop zone size from changing rapidly
                 * drDest = drad * MAX( 0.8 , 0.20 /dDestRate )
                 * >>chng 96 jan 14 from .2 to .25 to use less zones
                 * >>chng 96 may 30 from min of 0.8 to 0.5 to prevent crashing into He i-front */
                drDest = radius.drad*MAX2(0.5,0.20/dDestRate);
                /* endif */
        }
        else
        {
                drDest = BigRadius;
        }
#       endif
        drDest = BigRadius;/*03 dec 12 is this needed? */

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* check whether change in wind velocity constrains DRAD 
         * WJH 22 May 2004: disable when we are near the sonic point since
         * the velocity may be jumping all over the place but we just want
         * to push through it as quickly as we can */
        if( wind.windv!=0. && !pressure.lgSonicPoint && !pressure.lgStrongDLimbo )
        {
                double v = fabs(wind.windv);
                /* this is relative change in velocity */
                dVelRelative = fabs(wind.windv-OldWindVelocity)/
                        MAX2(v,0.1*timesc.sound_speed_isothermal);

#               define WIND_CHNG_VELOCITY_RELATIVE      0.01
                /*fprintf(ioQQQ,"DEBUG rad %.3f vel %.2e dVelRelative %.2e", 
                        log10(radius.Radius) , wind.windv ,  dVelRelative );*/

                /* do not use this on first zone since do not have old velocity */
                if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE  && nzone>1 )
                {
                        /* factor less than one if change larger than WIND_CHNG_VELOCITY_RELATIVE */
                        double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative;
                        /*fprintf(ioQQQ," DEBUG factor %.2e", factor );*/
                        factor = MAX2(0.8 , factor );
                        winddr = radius.drad * factor;
                }
                else
                {
                        winddr = BigRadius;
                }
                /*fprintf(ioQQQ," DEBUG \n" );*/

                WindAccelDR = BigRadius;

                /* set dr from advective term,
                 * the 1/20 came from looking at one set of structure plots */
                if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax+1)
                {
                        /* >>chng 02 dec 11, from 5 to 20 */
                        winddr = MIN2( winddr , dynamics.dRad / 20. );
                        /*>>chng 04 oct 06, set dVelRelative to dynamics.dRad since dVelRelative is printed as part
                         * of reason for choosing this criteria, want to reflect valid reason. */
                        dVelRelative = dynamics.dRad;
                }
        }
        else
        {
                winddr = BigRadius;
                WindAccelDR = BigRadius;
                dVelRelative = 0.;
        }
        /* remember old velocity */
        OldWindVelocity = wind.windv;

        /* inside out globule */
        if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
        {
#               define  DNGLOB  0.10
                if( radius.glbdst < 0. )
                {
                        fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured,  sorry.\n" );
                        fprintf( ioQQQ, " This is routine radius_next, GLBDST is%10.2e\n", 
                          radius.glbdst );
                        cdEXIT(EXIT_FAILURE);
                }
                factor = radius.glbden*pow(radius.glbrad/radius.glbdst,radius.glbpow);
                fac2 = radius.glbden*pow(radius.glbrad/(radius.glbdst - (realnum)radius.drad),radius.glbpow);
                if( fac2/factor > 1. + DNGLOB )
                {
                        /* DNGLOB is relative change in density allowed this zone, 0.10 */
                        GlobDr = radius.drad*DNGLOB/(fac2/factor - 1.);
                }
                else
                {
                        GlobDr = BigRadius;
                }
                /* GlobDr = GLBDST * DNGLOB * (GLBRAD/GLBDST)**(-GLBPOW) / GLBPOW */
                GlobDr = MIN2(GlobDr,radius.glbdst/20.);
        }
        else
        {
                GlobDr = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        hdnew = 0.;
        if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 )
        {
                /* one of the special density laws, first get density at possible next dr */
                if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
                {
                        hdnew = dense_fabden(radius.Radius+radius.drad,radius.depth+
                          radius.drad);
                }
                else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
                {
                        hdnew = dense_tabden(radius.Radius+radius.drad,radius.depth+
                          radius.drad);
                }
                else
                {
                        fprintf( ioQQQ, " dlw insanity in radius_next\n" );
                        cdEXIT(EXIT_FAILURE);
                }
                drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]);
                drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab));
        }
        else
        {
                drTab = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* special density law */
        if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
        {
                dnew = fabs(dense_fabden(radius.Radius+radius.drad,radius.depth+
                  radius.drad)/dense.gas_phase[ipHYDROGEN]-1.);
                /* DNGLOB is relative change in density allowed this zone, 0.10 */
                if( dnew == 0. )
                {
                        SpecDr = radius.drad*3.;
                }
                else if( dnew/DNGLOB > 1.0 )
                {
                        SpecDr = radius.drad/(dnew/DNGLOB);
                }
                else
                {
                        SpecDr = BigRadius;
                }
        }
        else
        {
                SpecDr = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* check grain line heating dominates
         * this is important in PDR and HII region calculations
         * >>chng 97 jul 03, added following check */
        if( thermal.heating[0][13]/thermal.htot > 0.2 )
        {
                /* >>chng 01 jan 03, following returns 0 when NO light at all,
                 * had failed with botched assert */
                GrainRateDr(&grfreqm,&gropacm);
                DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin);
        }
        else
        {
                gropacm = 0.;
                grfreqm = 0.;
                DrGrainHeat = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* check if line heating dominates
         * this is important in high metallicity models */
        if( thermal.heating[0][22]/thermal.htot > 0.2 )
        {
                FndLineHt(&level,&ipStrong,&Strong);
                if( Strong/thermal.htot > 0.1 )
                {
                        if( level == 1 )
                        {
                                /* a level1 line was the heat source (usual case)
                                 * drLineHeat = MAX(1.0,TauLines(ipLnTauIn,ipStrong)*0.2) /
                                 *  1      TauLines(ipLnDTau,ipStrong) */
                                TauInwd = TauLines[ipStrong].Emis->TauIn;
                                TauDTau = TauLines[ipStrong].Emis->PopOpc * TauLines[ipStrong].Emis->opacity / 
                                        DoppVel.doppler[TauLines[ipStrong].Hi->nelem-1];
                        }
                        else if( level == 2 )
                        {
                                /* a atom_level2 line was the heat source
                                 * (bad case since not as well treated)
                                 * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
                                 *  1      WindLine(ipLnDTau,ipStrong) */
                                TauInwd = TauLine2[ipStrong].Emis->TauIn;
                                TauDTau = TauLine2[ipStrong].Emis->PopOpc * TauLine2[ipStrong].Emis->opacity / 
                                        DoppVel.doppler[TauLine2[ipStrong].Hi->nelem-1];
                        }
                        else if( level == 3 )
                        {
                                /* a 12CO line was the heat source
                                 * (bad case since not as well treated)
                                 * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
                                 *  1      WindLine(ipLnDTau,ipStrong) */
                                TauInwd = C12O16Rotate[ipStrong].Emis->TauIn;
                                TauDTau = C12O16Rotate[ipStrong].Emis->PopOpc * C12O16Rotate[ipStrong].Emis->opacity / 
                                        DoppVel.doppler[C12O16Rotate[ipStrong].Hi->nelem-1];
                        }
                        else if( level == 4 )
                        {
                                /* a 13CO line was the heat source
                                 * (bad case since not as well treated)
                                 * drLineHeat = MAX(1.0,WindLine(ipLnTauIn,ipStrong)*0.2) /
                                 *  1      WindLine(ipLnDTau,ipStrong) */
                                TauInwd = C13O16Rotate[ipStrong].Emis->TauIn;
                                TauDTau = C13O16Rotate[ipStrong].Emis->PopOpc * C13O16Rotate[ipStrong].Emis->opacity / 
                                        DoppVel.doppler[C13O16Rotate[ipStrong].Hi->nelem-1];
                        }
                        else if( level == 5 )
                        {
                                /* >>chng 03 dec 07, this branch had been left off, caught by Hiroaki Oyaizu */
                                /* a hyperfine transition */
                                TauInwd = HFLines[ipStrong].Emis->TauIn;
                                TauDTau = HFLines[ipStrong].Emis->PopOpc * HFLines[ipStrong].Emis->opacity / 
                                        DoppVel.doppler[HFLines[ipStrong].Hi->nelem-1];
                        }
                        else
                        {
                                /* this is insane, since Strong was set, but not level */
                                fprintf( ioQQQ, " PROBLEM radius_next Strong line heating set, but not level.\n" );
                                TotalInsanity();
                        }
                        /* in following logic cannot use usual inverse opacity,
                         * since line heating competes with escape probability,
                         * so is effective at surprising optical depths */
                        if( TauDTau > 0. )
                        {
                                drLineHeat = MAX2(1.,TauInwd)*0.4/TauDTau;
                        }
                        else
                        {
                                drLineHeat = BigRadius;
                        }
                }
                else
                {
                        TauInwd = 0.;
                        drLineHeat = BigRadius;
                        ipStrong = 0;
                        Strong = 0.;
                }
        }
        else
        {
                TauInwd = 0.;
                drLineHeat = BigRadius;
                ipStrong = 0;
                level = 0;
                Strong = 0.;
        }

        /* >>chng 03 mar 03, add this logic */
        /* do not let change in cooling/heating due to H2 become too large */
        drH2_heat_cool = BigRadius;
        if( lgFirstCall )
        {
                Old_H2_heat_cool = 0.;
        }
        else if( !thermal.lgTemperatureConstant )
        {
                /* this is case where temperature has not been set */
                /* compare total heating - cooling due to h2 with total due to everything */
                double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
                if( H2_heat_cool > 0.1 )
                {
                        dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
                        dH2_heat_cool = SDIV(dH2_heat_cool);
                        /* >>chng 05 oct 27, had been taking 20% of original radius - this caused zoning
                         * to become very fine and may not have been needed - change from 0.2 to 0.5 */
                        /*drH2_heat_cool = radius.drad*MAX2(0.2,0.05/dH2_heat_cool);*/
                        drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
                }
                else
                {
                        drH2_heat_cool = BigRadius;
                }
        }
        Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;

        /* >>chng 03 mar 04, add this logic */
        /* do not let change in H2 and heavy element molecular abundances become too large 
         * "change in heav ele mole abundance, d(mole)/elem" */
        drH2_abund = BigRadius;
        dr_mole_abund = BigRadius;
        mole_dr_change = -1;
        if( nzone>=4 )
        {
                /* first do H2 abundance */
                double Old_H2_abund;
                /* >>chng 04 dec 15, do not use special logic when large h2 is turned on */
                int i;
                Old_H2_abund = struc.H2_molec[ipMH2g][nzone-3] + struc.H2_molec[ipMH2s][nzone-3];
                /* >>chng 03 jun 16, limit from 0.01 to 0.001, some models fell over H2 front due to
                 * large zoning, when large H2 was just this caused oscillations in solomon process */
                /* >>chng 03 nov 22, from > 0.001 to > 3e-4, models that start almost in H2 have
                 * rapid increase in H2 at shallow depths, start sensing this sooner */
                /* >>chng 03 dec 10, from 3e-4 to 1e-4, to get smaller chagned in leiden1.in test */
                /* radius_next keys from change in H2 abundance, d(H2)/H */
                /* >>chng 04 mar 11, start sensing H2 earlier since otherwise step size
                 * needs to become way too small tooo quickly.  limit changed from 1e-4 to 1e-6 */
                /* >>chng 04 jun 29, fromo > 1e-6 to >1e-8, some pdr's had too large a change in H2 */

                if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
                {
                        double fac = 0.2;
                        /* this is percentage change in H2 density - "change in H2 abundance" */
                        dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
                        /* in testing th85ism the dH2_abund did come out zero */
                        /* >>chng 03 jun 16, change d(H2) from 0.05 to 0.1, fine resolution of H2/H exposed
                         * small numerical oscillations in Solomon process */
                        /* >>chng 03 nov 22, smallest possible ratio of dr(next)/dr changed from
                         * 0.2 to 0.05, models that started almost in H2 front need to be able to sense it */
                        /*drH2_abund = radius.drad*MAX2(0.05,fac/SDIV(dH2_abund) );*/
                        /* >>chng 04 mar 15, with such small possible changes in dr, only 0.05, a thermal front
                         * can easily cause large changes in H2 and T that are not due to zoning, but to the
                         * discontinuity.  make smallest change larger to prevent hald due to dr */
                        /* >>chng 05 aug 23, thermal front allowed dr to become much too small
                         * chng from 0.02 to .6 */
                        dH2_abund = SDIV(dH2_abund);
                        drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
                }
                else
                {
                        drH2_abund = BigRadius;
                }

                /* check how molecular fractions of all heavy elements are chaning relative 
                 * to total gas phase abundance */
                dr_mole_abund = BigRadius;
                /* >>chng 04 jun 02, upper limit had been all species but now limit to real
                 * molecules since do not want this logic to work with the ions */
                for(i=0; i<mole.num_comole_calc; ++i)
                {
                        realnum abund,
                                abund_limit;
                        if( !dense.lgElmtOn[COmole[i]->nelem_hevmol] || COmole[i]->n_nuclei <= 1 )
                                continue;
                        /* >>chng 03 sep 21, add CO logic */
                        /* >>chng 04 mar 30, generalize to any molecule at all */
                        /* >>chng 04 mar 31 lower limit to abund had been 0.01, change
                         * to 0.001 to pick up approach to molecular fronts */
                        abund = COmole[i]->hevmol/dense.gas_phase[COmole[i]->nelem_hevmol];
                        /* is this an ice?  need special logic for ice since density increases
                         * exponentially when grain temp falls below sublimation temperature 
                         * >>chng 05 dec 06 - detect changes for smaller abundances for ices
                         * due to large changes in ice abundances */
                        if( COmole[i]->lgGas_Phase )
                        {
                                abund_limit = 1e-3f;
                        }
                        else
                        {
                                /* this is an ice - track its abundance at lower values so that
                                 * we resolve the sublimation transition region */
                                abund_limit = 1e-5f;
                        }

                        if( abund > abund_limit )
                        {
                                double drmole_one, relative_change, relative_change_denom;
                                /* >>chng 05 dec 08, use smaller abundance for the denominator since just taking
                                 * current abundance will overlook case where current density is vastly large
                                 * than old density */
                                if( struc.CO_molec[i][nzone-3]>SMALLFLOAT )
                                {
                                        relative_change_denom = MIN2( COmole[i]->hevmol , struc.CO_molec[i][nzone-3] );
                                }
                                else
                                {
                                        relative_change_denom = COmole[i]->hevmol;
                                }
                                /* the relative change in the abundance */
                                relative_change = 
                                        /*fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / COmole[i]->hevmol;*/
                                        /* >>chng 05 dec 08, use smaller abundance */
                                        fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / relative_change_denom;
                                /* >>chng 03 jun 16, change from 0.05 to 0.1, fine resolution of H2/H exposed
                                 * small numerical oscillations in Solomon process */
                                /* >>chng 04 jun 02, from 0.1 back to 0.05, more extensive CO etc network
                                 * caused oscillations in SiO abundance and Si Si+ density. */
                                /* >>chng 04 aug 03, from 0.05 to 0.035, leiden pdr model v2 had major
                                 * jump in eden */
                                /* >>chng 04 oct 18, from 0.035 back to 0.05, leiden pdr v2 actually due to having
                                 * PAHs in fully molecular limit (??), this caused cool flow pdr grid to trip on
                                 * too small dr */
                                relative_change = SDIV(relative_change);
                                /*>>chng 05 dec 08, relative_change must be less than one - with
                                 * revised logic above can be bigger than one */
                                if( relative_change > 1. )
                                        relative_change = 1./relative_change;
                                /*drmole_one = radius.drad*MAX2(0.2,0.035/relative_change );*/
                                /* >>chng 05 aug 23, thermal front allowed dr to become much too small
                                 * chng from 0.02 to .6 */
                                /*drmole_one = radius.drad*MAX2(0.2,0.05/relative_change );*/
                                drmole_one = radius.drad*MAX2(0.6,0.05/relative_change );
                                /* final dr will be the smallest we encounter */
                                if( drmole_one < dr_mole_abund )
                                {
                                        /* this is the dr used to set next dr - keep track of which moe was changing */
                                        dr_mole_abund = drmole_one;
                                        mole_dr_change = i;
                                        dCO_abund = relative_change;
                                }
                        }
                }
        }

        /* some consideration due to big H2 molecule */
        drSolomon_BigH2 = H2_DR();

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* can't make drmax large deep in model since causes feedback
         * oscillations with changes in heating or destruction rates
         * >>chng 96 oct 15, change from 5 to 10 */
        if( nzone < 5 )
        {
                /* >>chng 96 nov 11, had been 4 * drad up to 11, change to following
                 * to be similar to c90.01, max of 1.3 between 5 and 11 
                 * >>chng 04 oct 29  geometry.DirectionalCosin was ioncorrect applied to this factor  */
                drmax = 4.*radius.drad;
        }
        else
        {
                drmax = 1.3*radius.drad;
        }

        /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this
         * point if it occurs */
#       if 0
        /* look for oscillations in electron density or tempeature - freeze dr if these occur */
        dr_ne_oscil = BigRadius;
        dr_te_oscil = BigRadius;
        if( nzone >= 11 )
        {
                /* >>chng 96 oct 15, do not let zones increase if oscillations present */
                /* >>chng 96 oct 31, error to declare oscillation propto toler, the 
                 *heating cooling tolerance */
                realnum errorHC = POW2(conv.HeatCoolRelErrorAllowed);
                realnum errorNe = (realnum)POW2(conv.EdenErrorAllowed );
                limit = nzone -2;
                ASSERT( limit < struc.nzlim );
                for( k=nzone - 10; k < limit; k++ )
                {
                        /* check that square of change both chng sign and is
                         * greater than square of heat-tool error */
                        if( (struc.testr[k-1] - struc.testr[k])/struc.testr[k]*
                          (struc.testr[k] - struc.testr[k+1])/struc.testr[k] < 
                          -errorHC )
                        {
                                dr_te_oscil = radius.drad;
                        }
                        /* small residiual is to allow 0.01 rel error */
                        if( (struc.ednstr[k-1] - struc.ednstr[k])/struc.ednstr[k]*
                          (struc.ednstr[k] - struc.ednstr[k+1])/struc.ednstr[k] < 
                          -errorNe )
                        {
                                /* >>chng 96 oct 15, do not let zones increase if oscillations present */
                                /* radius_next keys from electron density oscillation*/
                                dr_ne_oscil = radius.drad;
                        }
                }
        }
        /* >>chng 05 apr 05, do not sense temp oscillation, so that we can move past this
         * point if it occurs */
        dr_ne_oscil = BigRadius;
        dr_te_oscil = BigRadius;
#       endif

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* check on several convergence criteria */

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        if( !conv.lgConvTemp )
        {
                drFail = radius.drad/2.;
        }
        else
        {
                drFail = BigRadius;
        }

        /* time dependent recombination - conditions become very homogeneous and 
         * crash into I fronts - use last iteration's  radius as guess of current case */
        recom_dr_last_iter = BigRadius;
        if( dynamics.lgStatic && dynamics.lgRecom )
        {
                /* first time through nzone == 1 */
                static long int nzone_recom = -1;
                if( nzone<=1 )
                {
                        /* initialization case */
                        nzone_recom = 0;
                }
                else if( radius.depth < struc.depth_last[struc.nzone_last-1] && 
                        nzone_recom < struc.nzone_last )
                {
                        ASSERT( nzone_recom>=0 && nzone_recom<struc.nzone_last );
                        /* case where we are within previous computed structure 
                         * first possibly increase nzone_recom so that it points deeper
                         * than this zone */
                        while( struc.depth_last[nzone_recom] < radius.depth &&
                                nzone_recom < struc.nzone_last-1 )
                                ++nzone_recom;
                        recom_dr_last_iter = struc.drad_last[nzone_recom]*3.;
                }
                else
                {
                        /* case where we have overrun the previous iteration's saved structure */
                        recom_dr_last_iter = BigRadius;
                }
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* change in electron density - radius_next keys from change in elec den,
         * remember old electron density during first call */
        /* this is low ionization solution */
        if( nzone > 2 )
        {
                /* next is-2 since nzone is on physics not c scale, and we want previous zone */
                Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3];
                Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
                dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;

                if( dEfrac > SMALLFLOAT )
                {
                        double fac = 0.04;
                        /* >>chng 03 dec 25, use smaller rel change in elec frac when most elec in ipMole or grains */
                        /* >>chng 04 sep 14, change to from from metals but comment out */
                        /* >>chng 04 sep 17, change to from from metals - uncomment */
                        if( dense.eden_from_metals > 0.5 )
                        {
                                /* >>chng 04 sep 18, change 0.02 from 0.01 */
                                /* >>chng 04 sep 18, change 0.02 from 0.04 */
                                fac = 0.04;
                        }
                        /* >>chng 04 feb 23, add test on hydrogen being predomintly
                         * recombined due to three-body recom, which is very sensitive
                         * to the electron density - but only do this in partially ionized medium */
                        else if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 && 
                                dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]>0.1 &&
                                dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]<0.8 )

                        {
                                fac = 0.02;
                        }
                        /* >>chng 04 sep 17, change to 0.1 from 0.2 */
                        drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
                }
                else
                {
                        drEfrac = BigRadius;
                }
        }
        else
        {
                dEfrac = 0.;
                drEfrac = BigRadius;
                Efrac_old = 0.;
                Efrac_new = 0.;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* do not let thickness get too large
         *  1    '' radius_next keys from relative depth'',e11.3)') */
        if( nzone > 20 )
        {
                /*drDepth = radius.depth/20.;*/
                /* >>chng 02 nov 05, change from 1/20 to 1/10 wasted zones early on */
                drDepth = radius.depth/10.;
        }
        else
        {
                drDepth = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* case where stopping thickness or edge specified, need to approach slowly */
        thickness_total = BigRadius;
        DepthToGo = BigRadius;
        if( StopCalc.HColStop < 5e29 )
        {
                double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
                /* >>chng 03 oct 28, forgot to div col den above by eff density */
                thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
        }

        if( StopCalc.colpls < 5e29 )
        {
                double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff );
                thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
        }

        if( StopCalc.col_h2 < 5e29 )
        {
                /* >>chng 03 apr 15, add this molecular hydrogen */
                double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.col_h2-colden.colden[ipCOL_H2g]-colden.colden[ipCOL_H2s]) / coleff );
                thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
        }

        if( StopCalc.col_h2_nut < 5e29 )
        {
                /* >>chng 03 apr 15, add this molecular hydrogen */
                double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.col_h2_nut-(2*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])+dense.xIonDense[ipHYDROGEN][1])) / coleff  );
                thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
        }

        if( StopCalc.col_H0_ov_Tspin < 5e29 )
        {
                /* >>chng 05 jan 09, add n(H0)/Tspin */
                double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff  );
                thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
        }

        if( StopCalc.col_monoxco < 5e29 )
        {
                /* >>chng 03 apr 15, add this, CO */
                double coleff = (double)SDIV( (findspecies("CO")->hevmol)*geometry.FillFac);
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.col_monoxco-findspecies("CO")->hevcol) / coleff );
                thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
        }

        if( StopCalc.colnut < 5e29 )
        {
                double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff );
                thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
        }

        /* this is case where outer radius is set */
        if( radius.router[iteration-1] < 5e29 )
        {
                thickness_total = MIN2(thickness_total , radius.router[iteration-1] );
                DepthToGo = MIN2(DepthToGo ,
                        radius.router[iteration-1] - radius.depth );
        }

        /* this is case where stopping optical depth was specified */
        if( StopCalc.iptnu != rfield.nupper )
        {
                /* end optical depth has been specified */
                double dt = SDIV(opac.opacity_abs[StopCalc.iptnu-1]*geometry.FillFac);
                DepthToGo = MIN2(DepthToGo ,
                        (StopCalc.tauend-opac.TauAbsGeo[0][StopCalc.iptnu-1])/dt);
        }
        /* stop AV - usually this is dust, but we consider all opacity sources,
         * so always include this */
        /* compute some average grain properties */
        AV_to_go = BigRadius;
        if( rfield.opac_mag_V_extended > SMALLFLOAT && rfield.opac_mag_V_point>SMALLFLOAT )
        {
                /* by default stop av is very large, and opacity can be very small, so ratio
                 * goes to inf - work with logs to see how big the number is */
                double ave = log10(StopCalc.AV_extended - rfield.extin_mag_V_extended) - 
                        log10(rfield.opac_mag_V_extended);
                double avp = log10(StopCalc.AV_point - rfield.extin_mag_V_point) - 
                        log10(rfield.opac_mag_V_point);
                AV_to_go = MIN2( ave , avp );
                if( AV_to_go < 37. )
                {
                        AV_to_go = pow(10., AV_to_go );
                        /* this is to make sure that we go slightly deeper than AV so that
                         * we ttrigger this stop */
                        AV_to_go *= 1.0001;
                }
                else
                        AV_to_go = BigRadius;
                /*fprintf(ioQQQ,"DEBUG next dr %.3e %.3e %.3e\n", AV_to_go , ave, avp );*/
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* set dr if one of above tests have triggered */
        if( DepthToGo <= 0. )
        {
                TotalInsanity();
        }
        else if( DepthToGo < BigRadius )
        {
                /* want to approach the outer edge slowly,
                 * the need for this logic is most evident in brems.in - 
                 * HI fraction varies across coronal model */
                drThickness = MIN2( thickness_total/10. , DepthToGo );
        }
        else
        {
                drThickness = BigRadius;
        }

        /*fprintf(ioQQQ,
                "DEBUG depth2go z%li drThickness2 = %e %e %e\n",
                nzone , drThickness , thickness_total , DepthToGo );*/

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* spherical models, do not want delta R/R big */
        drSphere = radius.Radius*0.04;

        /* optical depth to electron scattering */
        /* >>chng 04 jun 16, add filling factor, was missing */
        dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
        /* >>chng 02 oct 06, increase dr when constant temperature,
         * to prevent some ct models from taking too many cells */
        if( thermal.lgTemperatureConstant ) 
                dRTaue *= 3.;

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        if( dense.flong != 0. )
        {
                drFluc = 0.628/2./dense.flong;
                /* >>chng 04 sep 18, caused cautions that ionization jumps occurred.
                 * set to have the value */
                drFluc /= 2.;
        }
        else
        {
                drFluc = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /* if density fluctuations in place then override change in heat
         * for dr set */
        if( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn )
        {
                drdHdStep = BigRadius;
        }

        /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
        /*active dr sets */
        /* we are deep into model, use logic since we have several zones
         * of old data */
        radius.drNext = MIN2( MIN4( drmax, drInter, drLineHeat, winddr ),
                              MIN4( drFluc, GlobDr, DrGrainHeat, dr_change_heavy ) );
        radius.drNext = MIN2( MIN4( radius.drNext, SpecDr, drFail, WindAccelDR ),
                              MIN3( drSphere, radius.sdrmax, dRTaue ) );
        radius.drNext = MIN3( MIN3( radius.drNext, drDest, drdTdStep ),
                              MIN3( drdHdStep, drConPres, drTab ),
                              MIN3( drSolomon_BigH2, drLeiden_hack, recom_dr_last_iter ) );
        radius.drNext = MIN3( MIN4( radius.drNext, drHion, drDepth, dr_mole_abund ),
                              MIN4( AV_to_go, drEfrac, drHMase, drThickness ),
                              MIN3( drHe1ovHe2, drH2_heat_cool, drH2_abund ) );

        /*fprintf(ioQQQ,
                "DEBUG depth2go drNext = %e \n",                radius.drNext );*/

        /* keep dr constant in first two zones, if it wants to increase,
         * but always allow it to decrease.
         * to guard against large increases in efrac for balmer cont photo dominated models,
         */
        if( nzone <= 1 && radius.drNext > OldDR )
        {
                radius.drNext = OldDR;
        }

        /* a pressure failure has occurred - keep zone the same time, hoping to pass through
         * troubled region */
        if( !conv.lgConvPres || !conv.lgConvTemp )
        {
                radius.drNext = radius.drad;
        }

        /* >>chng 05 aug 05, in case of thermal front, where temp is falling quickly and
         * conditions change very fast, the zone thickness does not really affect the change
         * in conditions and can cause zoning to become very very thin, which causes
         * an abort.  this occurs between 200 and 1000K.  if we are doing temp soln,
         * temp is between these values, and temp is changing rapidly, do not make sone
         * thickness much smaller */
        drThermalFront = BigRadius;
        /* do not use thermal front logic in case of recombination time dep cloud */
        if( nzone >=5 && !dynamics.lgStatic )
        {
                /* >>chng 05 aug 23, upper bound of thermal from from 1000K to 4000K */
                /*if( phycon.te > 200. && phycon.te < 1000. && */
                if( phycon.te > 200. && phycon.te < 3000. && 
                        /* >>chng 05 aug 23, from > 10% in zone to to two zones > 5%,
                         * to fix leiden v3 with large H2 */
                        (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
                        (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
                        (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
                {
                        /* the 0.91 is to make dr unique, so that print statement that
                         * follows will identify this reason */
                        drThermalFront = radius.drad * 0.91;
                        radius.drNext = drThermalFront;
                }
        }

        /* dr = zero is a logical mistake */
        if( radius.drNext <= 0. )
        {
                fprintf( ioQQQ, " radius_next finds insane drNext:%10.2e\n", 
                  radius.drNext );
                fprintf( ioQQQ, " all drs follow:%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e\n", 
                  drmax, drInter, drLineHeat, winddr, drFluc, GlobDr, SpecDr, 
                  drFail, drSphere, radius.sdrmax, dRTaue, 
                  OldH2Abund, drDepth );
                cdEXIT(EXIT_FAILURE);
        }

        /* option to force min drad */
        if( radius.drNext < radius.sdrmin )
        {
                radius.drNext = radius.sdrmin;
        }

        /* this is general code that prevents zone thickness drNext from
         * becoming too thin, something that can happen for various bad reasons
         * HOWEVER we do not want to do this test for certain density laws,
         * for which very small zone thicknesses are unavoidable
         * the special cases are:
         * special density law,
         * globule density law,
         * carbon +-0 i front
         * the fluctuations command
         * drMinimum was set in radius_first to either sdrmin (set drmin) or
         * some fraction of the initial radius - it is always set
         * to something non-trivial.  
         * sdrmin is set with the "set dr" command - is SMALLFLOAT by default */
        if( ((strcmp(dense.chDenseLaw,"DLW1") != 0 && 
                strcmp(dense.chDenseLaw ,"GLOB") != 0) )&& 
                /* >>chng 04 feb 19, do not use this test - errors can still happen
                 * when all C is atomic! */
                /*dense.xIonDense[ipCARBON][0]/dense.gas_phase[ipCARBON] < 0.05) && */
                (dense.flong == 0.) &&
                /* >>chng 01 aug 11, add check against stopping on depth to go */
                radius.drNext != DepthToGo )
        {
                /* don't let dr get smaller than drMinimum, if this resets drNext
                 * then code stops with warning that zones got too thin
                 * this can happen due to numerical oscillations, which the code
                 * tries to damp out by making the zones thinner.
                 * scale by radius.Radius/radius.rinner to stop very spherical 
                 * simulations from false trigger on dr too small */
                if( radius.drNext* radius.Radius/radius.rinner  * 
                        /* drMinimum is drad * hden, to make it proportional to optical depth.
                         * This avoids false trigger across thermal fronts. */
                        dense.gas_phase[ipHYDROGEN] < 
                        radius.drMinimum )
                {
                        radius.drNext = radius.drMinimum/ dense.gas_phase[ipHYDROGEN];
                        /* leaving this at true will cause the model to stop with a warning
                         * that the zone thickness is too small */
                        radius.lgDrMinUsed = true;
                        /* set abort handler */
                        lgAbort = true;
                        /* must decrement nzone, since we will not complete this zone, and will not have
                         * valid structure data for it */
                        --nzone;
                        fprintf( ioQQQ, 
                                "\n DISASTER PROBLEM radius_next finds dr too small and aborts.  "
                                "This is zone %ld iteration %ld.  The thickness was %.2e\n", 
                                nzone, 
                                iteration,
                                radius.drNext);
                        fprintf( ioQQQ, 
                                " If this simulation seems reasonable you can override this limit "
                                "with the command SET DRMIN %.2e\n\n", 
                                radius.drNext*2);
                        /* set abort flag */
                        lgAbort = true;
                        return 1;
                }
        }

        /* factor to allow for slop in floating numbers */
        const double Z = 1.0001;

        /* following is to make thickness of model exact
         * n.b., on last zone, drNext can be NEGATIVE!!
         * DEPTH was incremented at start of zone calc in newrad,
         * has been outer edge of zone all throughout */
        double drOuterRadius = (radius.router[iteration-1]-radius.depth)*Z;
        radius.drNext = min(radius.drNext,drOuterRadius);

        /* this means outer limit exceeded */
        fixit(); /* is this useful? There is an assert below that drNext > 0. */
        if( radius.drNext < 0. )
        {
                radius.lgDrNeg = true;
        }
        else
        {
                radius.lgDrNeg = false;
        }

        /* set flag if dr set by maser */
        if( fp_equal( radius.drNext, drHMase ) )
        {
                rt.lgMaserSetDR = true;
        }

        /* all this is to only punch on last iteration
         * the punch dr command is not really a punch command, making this necessary
         * lgDRon is set true if "punch dr" entered */
        if( punch.lgDROn )
        {
                if( punch.lgDRPLst )
                {
                        /* lgDRPLst was set true if "punch" had "last" on it */
                        if( iterations.lgLastIt )
                        {
                                lgDoPun = true;
                        }
                        else
                        {
                                lgDoPun = false;
                        }
                }
                else
                {
                        lgDoPun = true;
                }
        }
        else
        {
                lgDoPun = false;
        }
        if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
        {
                if( !conv.lgConvTemp && nzone > 0 )
                {
                        fprintf( punch.ipDRout, "#>>>> A temperature failure occured.\n" );
                }
                if( !conv.lgConvPres && nzone > 0 )
                {
                        fprintf( punch.ipDRout, "#>>>> A pressure failure occured.\n" );
                }

                /* this is common part of each line, the zone count, depth, chosen dr, and depth2go */
                fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drNext, radius.Depth2Go );

                /*=======begin active dr sets */
                if( fp_equal( radius.drNext, drLineHeat ) )
                {
                        if( level == 1 )
                        {
                                strcpy( chLbl, chLineLbl(&TauLines[ipStrong]) );
                                fprintf( punch.ipDRout, "level 1 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", 
                                  chLbl, TauInwd, TauLines[ipStrong].Emis->pump, 
                                  TauLines[ipStrong].Emis->Pesc );
                        }
                        else if( level == 2 )
                        {
                                strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]));
                                fprintf( punch.ipDRout, "level 2 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n", 
                                  chLbl, TauInwd, TauLine2[ipStrong].Emis->pump, 
                                  TauLine2[ipStrong].Emis->Pesc );
                        }
                        else
                        {
                                fprintf( ioQQQ, " insanity pr line heat\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( fp_equal( radius.drNext, drDepth ) )
                {
                        fprintf( punch.ipDRout, "relative depth\n");
                }

                else if( fp_equal( radius.drNext, drThermalFront ) )
                {
                        fprintf( punch.ipDRout, "thermal front logic\n");
                }

                else if( fp_equal( radius.drNext, dr_change_heavy ) )
                {
                        ASSERT( ichange_heavy_nelem < LIMELM+3 );
                        fprintf( punch.ipDRout, 
                                "change in ionization, element %s ion (C scale) %li rel change %.2e ion frac %.2e\n",
                                elementnames.chElementName[ichange_heavy_nelem],
                                ichange_heavy_ion , 
                                change_heavy_frac_big ,
                                frac_change_heavy_frac_big);
                }

                else if( fp_equal( radius.drNext, drThickness ) )
                {
                        fprintf( punch.ipDRout, "depth to go\n");
                }

                else if( fp_equal( radius.drNext, AV_to_go ) )
                {
                        fprintf( punch.ipDRout, "A_V to go\n");
                }

                else if( fp_equal( radius.drNext, drTab ) )
                {
                        fprintf( punch.ipDRout, "spec den law, new old den%10.2e%10.2e\n", 
                          hdnew, dense.gas_phase[ipHYDROGEN] );
                }

                else if( fp_equal( radius.drNext, drHMase ) )
                {
                        fprintf( punch.ipDRout, "H maser dTauMase=%10.2e %li %li %li %li\n", 
                                rt.dTauMase,
                                rt.mas_species,
                                rt.mas_ion,
                                rt.mas_hi,
                                rt.mas_lo );
                }

                else if( fp_equal( radius.drNext, drHe1ovHe2 ) )
                {
                        /* radius_next keys from change in He2/He1 ionization, old=%11.3e sv=%11.3e FrLya*/
                        fprintf( punch.ipDRout, "change in He0/He+ ionization, ratio %.2e\n", 
                          Old_He_atom_ov_ion );
                }

                else if( fp_equal( radius.drNext, drH2_heat_cool ) )
                {
                        /* radius_next keys from change in H2 heating, old=%11.3e sv=%11.3e FrLya*/
                        fprintf( punch.ipDRout, "change in H2 heating/cooling, d(c,h)/H %.2e\n", 
                          dH2_heat_cool );
                }

                else if( fp_equal( radius.drNext, drH2_abund ) )
                {
                        /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/
                        fprintf( punch.ipDRout, "change in H2 abundance, d(H2)/H %.2e\n", 
                          dH2_abund );
                }

                else if( fp_equal( radius.drNext, dr_mole_abund ) )
                {
                        /* radius_next keys from change in CO mole abundance */
                        fprintf( punch.ipDRout, "change in heav ele mole abundance, d(mole)/elem %.2e mole=%i=%s\n", 
                          dCO_abund , mole_dr_change , COmole[mole_dr_change]->label);
                }

                else if( fp_equal( radius.drNext, drSolomon_BigH2 ) )
                {
                        /* radius_next keys from change in H2 abundance, old=%11.3e sv=%11.3e FrLya*/
                        fprintf( punch.ipDRout, "change in big H2 Solomon rate line opt depth\n");
                }

                else if( fp_equal( radius.drNext, recom_dr_last_iter ) )
                {
                        /* radius_next keys from previous iteration recom logic */
                        fprintf( punch.ipDRout, "previous iteration recom logic\n");
                }

                else if( fp_equal( radius.drNext, drHion ) )
                {
                        fprintf( punch.ipDRout, "change in H ionization fm to%11.3e%11.3e\n", 
                          SaveOHIIoHI, OHIIoHI );
                }

                else if( fp_equal( radius.drNext, drEfrac ) )
                {
                        fprintf( punch.ipDRout, 
                                "change in elec den, rel chng:%11.3e, cur %11.3e old=%11.3e\n", 
                          dEfrac , Efrac_old , Efrac_new );
                }

                else if( fp_equal( radius.drNext, drdHdStep ) )
                {
                        fprintf( punch.ipDRout, 
                                "change in heating; current %10.3e delta=%10.3e\n", 
                          thermal.htot, dHdStep );
                }

                else if( fp_equal( radius.drNext, drLeiden_hack ) )
                {
                        fprintf( punch.ipDRout, 
                                "Leiden hack\n" );
                }

                else if( fp_equal( radius.drNext, drConPres ) )
                {
                        fprintf( punch.ipDRout, "change in con accel\n"  );
                }

                else if( fp_equal( radius.drNext, drdTdStep ) )
                {
                        fprintf( punch.ipDRout, 
                                "change in temperature; current= %10.3e, dT/T= %.3f\n", 
                          phycon.te, dTdStep );
                }

                else if( fp_equal( radius.drNext, radius.sdrmin ) )
                {
                        fprintf( punch.ipDRout, "sdrmin\n"  );
                }

                else if( fp_equal( radius.drNext, radius.sdrmax ) )
                {
                        fprintf( punch.ipDRout, "sdrmax\n" );
                }

                else if( fp_equal( radius.drNext, drSphere ) )
                {
                        fprintf( punch.ipDRout, "sphericity\n" );
                }

                else if( fp_equal( radius.drNext, dRTaue ) )
                {
                        fprintf( punch.ipDRout, 
                                "optical depth to electron scattering\n" );
                }

                else if( fp_equal( radius.drNext, drFail ) )
                {
                        fprintf( punch.ipDRout, 
                                "temperature failure.\n" );
                }

                else if( fp_equal( radius.drNext, drmax ) )
                {
                        fprintf( punch.ipDRout, 
                                "DRMAX; nu opc dr=%10.2e%10.2e%10.2e\n", 
                          freqm, opacm, radius.drChange/
                          SDIV(opacm) );
                }

                else if( fp_equal( radius.drNext, drInter ) )
                {
                        fprintf( punch.ipDRout, 
                                "cont inter nu=%10.2e opac=%10.2e\n", 
                          freqm, opacm );
                }

                else if( fp_equal( radius.drNext, DrGrainHeat ) )
                {

                        fprintf( punch.ipDRout, 
                                "grain heating nu=%10.2e opac=%10.2e\n", 
                          grfreqm, gropacm );
                }

                else if( fp_equal( radius.drNext, winddr ) )
                {
                        fprintf( punch.ipDRout, 
                                "Wind, dVelRelative=%.3e windv=%.3e\n", 
                           dVelRelative ,
                           wind.windv);
                }

                else if( fp_equal( radius.drNext, drFluc ) )
                {
                        fprintf( punch.ipDRout, 
                                "density fluctuations\n"  );
                }

                else if( fp_equal( radius.drNext, GlobDr ) )
                {
                        fprintf( punch.ipDRout, 
                                "GLOB law new dr=%10.2e HDEN=%10.2e\n", 
                                GlobDr,
                          dense.gas_phase[ipHYDROGEN] );
                }

                else if( fp_equal( radius.drNext, SpecDr ) )
                {
                        fprintf( punch.ipDRout, 
                                "special law new dr=%10.2e HDEN=%10.2e\n", 
                                SpecDr,
                          dense.gas_phase[ipHYDROGEN] );
                }

                else if( fp_equal( radius.drNext, OldDR ) )
                {
                        fprintf( punch.ipDRout, "old DR.\n" );
                }

                else if( fp_equal( radius.drNext, drOuterRadius ) )
                {
                        fprintf( punch.ipDRout, "outer radius.\n" );
                }

                else
                {
                        fprintf( punch.ipDRout, 
                                " %4ld radius_next keys from insanity %10.2e\n", 
                          nzone, radius.drNext );

                        fprintf( ioQQQ, 
                                " %4ld radius_next keys from insanity %10.2e\n", 
                          nzone, radius.drNext );
                        cdEXIT(EXIT_FAILURE);
                }

                /*======end active dr sets */
        }

        ASSERT( radius.drNext > 0. );

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was%12.4e\n", 
                  radius.drNext, radius.drad );
        }

        return 0;
}

/*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
STATIC void ContRate(double *freqm, 
  double *opacm)
{
        long int i, 
          ipHi,
          iplow, 
          limit;
        double FreqH, 
          Freq_nonIonizing, 
          Opac_Hion, 
          Opac_nonIonizing, 
          Rate_max_Hion, 
          Rate_max_nonIonizing;

        DEBUG_ENTRY( "ContRate()" );

        /* 
         * find maximum continuum interaction rate,
         * these should be reset in following logic without exception,
         * if they are still zero at the end we have a logical error 
         */
        Rate_max_nonIonizing = 0.;
        Freq_nonIonizing = 0.;
        Opac_nonIonizing = 0.;

        /* this must be reset to val >= 0 */
        *opacm = -1.;
        *freqm = -1.;

        /* do up to carbon photo edge if carbon is turned on */
        /* >>>chng 00 apr 07, add test for whether element is turned on */
        if( dense.lgElmtOn[ipCARBON] )
        {
                /* carbon is turned on, use carbon 1 edge */
                ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
        }
        else
        {
                /* carbon turned off, use hydrogen Balmer continuum */
                ipHi = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]-1;
        }

        /* start at plasma frequency */
        for( i=rfield.ipPlasma; i < ipHi; i++ )
        {
                /* this does not have grain opacity since grains totally passive
                 * at energies smaller than CI edge */
                if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 
                  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
                {
                        Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
                          (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
                        Freq_nonIonizing = rfield.anu[i];
                        Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
                }
        }

        /* not every continuum extends beyond C edge-this whole loop can add to zero
         * use total opacity here
         * test is to put in fir continuum if free free heating is important */
        if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
        {
                /* this is index for energy where cloud free free optical depth is unity,
                 * is zero if no freq are opt thin */
                iplow = MAX2(1 , rfield.ipEnergyBremsThin );
        }
        else
        {
                /* >>>chng 00 apr 07, from Heavy.ipHeavy[0][5] to ipHi defined above, since
                 * would crash if element not defined */
                iplow = ipHi;
        }
        /* do not go below plasma frequency */
        iplow = MAX2( rfield.ipPlasma , iplow );

        /* this energy range from carbon edge to hydrogen edge */
        limit = MIN2(rfield.nflux,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1);
        for( i=iplow; i < limit; i++ )
        {
                if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 
                  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
                {
                        Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
                          (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
                        Freq_nonIonizing = rfield.anu[i];
                        Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
                }
        }

        /* variables to check continuum interactions over Lyman continuum */
        Rate_max_Hion = 0.;
        Opac_Hion = 0.;
        FreqH = 0.;

        /* not every continuum extends beyond 1 Ryd-this whole loop can add to zero */
        for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ )
        {
                if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 
                  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
                {
                        /* Rate_max_Hion = anu(i)*flux(i)/widflx(i)*opac(i) */
                        Rate_max_Hion = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
                          (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
                        FreqH = rfield.anu[i];
                        Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
                }
        }

        /* use Lyman continuum if its opacity is larger than non-h ion */
        if( Rate_max_nonIonizing < 1e-30  && Opac_Hion > SMALLFLOAT )
        {
                /* this happens for laser source - use Lyman continuum */
                *opacm = Opac_Hion;
                *freqm = FreqH;
        }
        /* >>chng 05 aug 03, add last test on Opac_Hion for case where we go very
         * deep and very little radiation is left */
        else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/Rate_max_nonIonizing > 1e-10 && Opac_Hion > SMALLFLOAT )
        {
                /* use Lyman continuum */
                *opacm = Opac_Hion;
                *freqm = FreqH;
        }
        else
        {
                /* not much rate in the Lyman continuum, stick with low energy */
                *opacm = Opac_nonIonizing;
                *freqm = Freq_nonIonizing;
        }

#       if 0
        /*>>chng 06 aug 12, do not let zones become optically thick to
        * peak of dust emission - one of NA's vastly optically thick dust clouds
        * did not conserve total energy - very opticallly thick to ir dust emission
        * so ir was absorbed and reemitted many times
        * this makes sure the cells remain optically thin to dust emission */
        if( gv.lgDustOn && gv.lgGrainPhysicsOn )
        {
                double fluxGrainPeak = -1.;
                long int ipGrainPeak = -1;
                long int ipGrainPeak2 = -1;

                i = 0;
                while( rfield.anu[i] < 0.0912  && i < (rfield.nflux-2) )
                {
                        /* >>chng 06 aug 23.  Only want to do this test for the IR dust continuum, therefore only 
                         * check on optical depth for wavelengths greater than 1 micron */
                        if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak )
                        {
                                ipGrainPeak = i;
                                fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i];
                        }
                        ++i;
                }
                ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );

                /* >>chng 06 aug 23.  We have just found the wavelength and flux at the peak of the IR emission.
                 * Now we want to find the wavelength, short of the peak, which corresponds to 5% of the 
                 * peak emission.  This wavelength will be where the code checks to make sure the zone has
                 * not become optically thick.  Since dust opacity generally decreases with wavelength,  this assures that
                 * the optical depths are small for all wavelengths where the flux is appreciable.  Tests show that
                 * this allows for flux/luminosity conservation to within ~5% for an AV of 1e4 mag and at least 2 iterations*/
                i = ipGrainPeak;
                /* >>chng 06 aug 23.  Only want to do this test for the IR dust continuum, therefore only 
                 * check on opacities for wavelengths shortward of the peak and greater than 1 micron 
                 * this routine can be called with the dust emission totally zero - it is only evaluated 
                 * late to save time.  don't do the following tests if peak dust emission is zero */
                if( fluxGrainPeak > 0. )
                {
                        while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
                        {
                                /* find wavelength where flux = 5% of peak flux, shortward of the peak */
                                if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
                                {
                                        /* This will be the array number and flux at 10% of the peak */
                                        ipGrainPeak2 = i;
                                }
                                ++i;
                        }
                        ASSERT( ipGrainPeak2>=0 );

                        /* use this as a limit to the dr if opacity is larger */
                        if( opac.opacity_abs[ipGrainPeak2] > *opacm )
                        {
                                /* scale opacity by factor of 10, which further decreases the zone thickness*/
                                *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
                                *freqm = rfield.anu[ipGrainPeak2];
                                /*fprintf(ioQQQ,"DEBUT opac peak %.2e %.2e \n",
                                        *opacm , *freqm );*/
                        }
                }
        }
#       endif

        {
                /* following should be set true to print contributors */
                enum {DEBUG_LOC=false};
                if( DEBUG_LOC )
                {
                        fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
                        Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
                        Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
                        );

                }
        }

        /* these were set to -1 at start, and must have been reset in one of the
         * two loops.  Logic error if still <0. */
        /* >>chng 05 aug 03, change logic to -1 on entry and check at least zero
         * here - will be zero if NO radiation field exists, perhaps since totally
         * attenuated */
        ASSERT( *opacm >= 0. && *freqm >= 0. );
        return;
}

/*GrainRateDr called by radius_next to find grain heating rate dr */
STATIC void GrainRateDr(double *grfreqm, 
  double *gropacm)
{
        long int i, 
          iplow;
        double xMax;

        DEBUG_ENTRY( "GrainRateDr()" );

        /* in all following changed from anu2 to anu  july 25 95
         *
         * find maximum continuum interaction rate */

        /* not every continuum extends beyond C edge-this whole loop can add to zero
         * use total opacity here
         * test is to put in fir continuum if free free heating is important */
        if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
        {
                /* this is pointer to energy where cloud free free optical depth is unity,
                 * is zero if no freq are opt thin */
                iplow = MAX2(1 , rfield.ipEnergyBremsThin );
        }
        else
        {
                /* do up to carbon photo edge if carbon is turned on */
                /* >>>chng 00 apr 07, add test for whether element is turned on */
                if( dense.lgElmtOn[ipCARBON] )
                {
                        /* carbon is turned on, use carbon 1 edge */
                        iplow = Heavy.ipHeavy[ipCARBON][0];
                }
                else
                {
                        /* carbon truned off, use hydrogen balmer continuum */
                        iplow = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s];
                }

        }

        xMax = -1.;
        /* integrate up to H edge */
        for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
        {
                if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
                {
                        xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
                          opac.opacity_abs[i];
                        *grfreqm = rfield.anu[i];
                        *gropacm = opac.opacity_abs[i];
                }
        }
        /* integrate up to heii edge if he is turned on,
         * this logic will not make sense if grains on but he off, which in itself makes no sense*/
        if( dense.lgElmtOn[ipHELIUM] )
        {
                for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
                {
                        if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
                        {
                                xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
                                  opac.opacity_abs[i];
                                *grfreqm = rfield.anu[i];
                                *gropacm = opac.opacity_abs[i];
                        }
                }
        }

        /* possible that there is NO ionizing radiation, in extreme cases,
         * if so then xMax will still be negative */
        if( xMax <= 0. )
        {
                *gropacm = 0.;
                *grfreqm = 0.;
        }
        return;
}

#undef WIND_CHNG_VELOCITY_RELATIVE

#if 0
/*ChkRate called by radius_next to check how rates of destruction of various species changes */
STATIC void ChkRate(
          /* element number on C scale */
          long int nelem, 
          /* change in destruction rate */
          double *dDestRate, 
          /* old and new destruction rates */
          double *DestRateOld,
          double *DestRateNew,
          /* stage of ionization on the physical scale */
          long int *istage)
{
        long int i;

        double average, 
          dDest;

        static double OldDest[LIMELM][LIMELM];

        DEBUG_ENTRY( "ChkRate()" );

        /* return if this element is not turned on */
        if( !dense.lgElmtOn[nelem] )
        {
                *dDestRate = 1e-3;
                *DestRateOld = 1e-3;
                *DestRateNew = 1e-3;
                *istage = 0;
                return;
        }

        /* for first zone, and during search, we will do nothing but
         * still must return finite numbers */
        *istage = 1;
        *dDestRate = 0.;
        *DestRateOld = 0.;
        *DestRateNew = 0.;
        *dDestRate = 0.;

        if( nzone <= 1 )
        {
                for( i=0; i < nelem+1; i++ )
                {
                        OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
                }
        }

        else if( dense.xIonDense[nelem][0]/dense.gas_phase[nelem] <  0.9 )
        {
                /* do not use this method if everything is atomic */
                for( i=0; i < (nelem); i++ )
                {
                        /* last check below, .5 chosen so that do not key off
                         * predominantly neutral species where self-opacity
                         * could cause oscillations */
                        if( ((dense.xIonDense[nelem][i]/dense.gas_phase[nelem] > 0.01 && 
                                dense.xIonDense[nelem][i]/dense.gas_phase[nelem] < 0.9) && 
                                dense.xIonDense[nelem][i+1]/dense.gas_phase[nelem] > .05) && 
                                OldDest[nelem][i] > 0. )
                        {
                                /* last check on old dest in case just lowered ionization
                                 * stage, so no history */
                                /* check that rate is positive */
                                if( ionbal.RateIonizTot[nelem][i] <= 0. )
                                {
                                        fprintf( ioQQQ, " ChkRate gets insane destruction rate for ion%4ld%4ld%10.2e\n", 
                                          nelem+1, i, ionbal.RateIonizTot[nelem][i] );
                                        cdEXIT(EXIT_FAILURE);
                                }

                                /* do not consider unless of middling ionization, and
                                 * rate is going down (to prevent dr osciallating)
                                 * no absolute value in following since do not want to
                                 * consider cases where ionization rate increases */
                                average = (OldDest[nelem][i] + ionbal.RateIonizTot[nelem][i])* 0.5;

                                dDest = (OldDest[nelem][i] - ionbal.RateIonizTot[nelem][i])/ average;
                                /* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ term + if rate going down */

                                if( *dDestRate < dDest )
                                {
                                        /* biggest rate so far, remember change in rates and ionization stage */
                                        *dDestRate = dDest;
                                        *istage = i+1;
                                        *DestRateOld = OldDest[nelem][i];
                                        *DestRateNew = ionbal.RateIonizTot[nelem][i];
                                }

                        }
                        OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
                }
        }
        return;
}
#endif

---