/home66/gary/public_html/cloudy/c08_branch/source/radius_increment.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_increment do work associated with geometry increments of this zone, called before RT_tau_inc */
/*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "hydrogenic.h"
#include "dynamics.h"
#include "rfield.h"
#include "hextra.h"
#include "colden.h"
#include "geometry.h"
#include "opacity.h"
#include "thermal.h"
#include "doppvel.h"
#include "dense.h"
#include "mole.h"
#include "h2.h"
#include "timesc.h"
#include "hmi.h"
#include "lines_service.h"
#include "taulines.h"
#include "trace.h"
#include "wind.h"
#include "phycon.h"
#include "pressure.h"
#include "grainvar.h"
#include "molcol.h"
#include "conv.h"
#include "hyperfine.h"
#include "mean.h"
#include "struc.h"
#include "radius.h"
/*cmshft - so Compton scattering shift of spectrum */
STATIC void cmshft(void);

#if !defined(NDEBUG)
/*pnegopc punch negative opacities on io unit, iff 'set negopc' command was 
 * given.  This function is only used in debug mode */
STATIC void pnegopc(void);
#endif

void radius_increment(void)
{
#       if !defined(NDEBUG)
        bool lgFlxNeg;
#       endif

        long int nelem,
          i, 
          j ,
          ion,
          nzone_minus_1 , 
          mol;
        double 
          avWeight,
          DilutionHere ,
          escatc, 
          fmol,
          opfac, 
          relec, 
          rforce, 
          t;

        double ajmass, 
                Error,
          rjeans;

        realnum dradfac;

        DEBUG_ENTRY( "radius_increment()" );

        /* when this sub is called radius is the outer edge of zone */

        if( lgAbort )
        {
                return;
        }

        /* following block only set of asserts to check that iso populations are sane */
#       if !defined(NDEBUG)
        if( !dynamics.lgAdvection )
        {
                long int ipISO;
                for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
                {
                        for( nelem=ipISO; nelem<LIMELM; ++nelem )
                        {
                                /* >>chng 05 feb 05, rm div by gas_phase for element, to eden, 
                                 * see explain for this date below */
                                if( dense.lgElmtOn[nelem] && 
                                        dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10. &&
                                        !(iso.chTypeAtomUsed[ipISO][nelem][0]=='L' && 
                                        iso.xIonSimple[ipISO][nelem]<1e-10) )
                                {
                                        /* we will check that ground pops are within this factor of the xIonDense value
                                         * these are only converged down to some extent - probably this should be related
                                         * to the ionization convergence error */
                                        /* >>chng 05 feb 05, for very highly ionized sims, where N or O is He-like,
                                         * the error can approach 1%, and depends on the convergence criteria.
                                         * the code does not try to converge these quantities in an absolute sense, 
                                         * only the quantities relative to the electron or hydrogen density.  so it is not
                                         * safe to do the assert for small values 
                                         * change is to only do this branch if abundance is large enough to be detected by
                                         * the ionization convergence solvers */
#                                       define ERR_CHK  1.001
                                        double err_chk = ERR_CHK;
                                        /* >>chng 05 sep 02, when low-T solver used solution is approximate,
                                         * and it must not matter (lot-T solver should not be used if it
                                         * does matter - so use more liberal expectation */
                                        if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' )
                                                err_chk *= 10.;
                                        /* make sure that populations are valid, first check Pop2Ion  */
                                        if( StatesElem[ipISO][nelem][0].Pop >
                                                dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk &&
                                                /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
                                                iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 )
                                        {
                                                fprintf(ioQQQ,
                                                        " DISASTER radius_increment finds inconsistent populations, Pop2Ion zn %.2f ipISO %li nelem %li %.4e %.4e solver:%s\n", 
                                                        fnzone,
                                                        ipISO , nelem ,
                                                        StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO] ,
                                                        dense.xIonDense[nelem][nelem-ipISO] ,
                                                        iso.chTypeAtomUsed[ipISO][nelem]);
                                                fprintf(ioQQQ," level solver %s found pop_ion_ov_neut of %.5e",
                                                        iso.chTypeAtomUsed[ipISO][nelem] ,
                                                        iso.pop_ion_ov_neut[ipISO][nelem] );
                                                fprintf(ioQQQ," ion_solver found pop_ion_ov_neut of %.5e\n",
                                                        dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) );
                                                fprintf(ioQQQ,
                                                        "simple %.3e Test shows Pop2Ion (%.6e) > xIonDense[nelem]/[nelem+1] (%.6e) \n", 
                                                        iso.xIonSimple[ipISO][nelem],
                                                        StatesElem[ipISO][nelem][0].Pop ,
                                                        dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO]) );
                                        }

                                        /* >>chng 06 jul 24, add assert that this is not search 
                                         * phase to confirm that removing other asserts on this 
                                         * are ok - we cannot be in search phase in this routine 
                                         * if ok, rm this and all ref to lgSearch, perhaps conv.h header */
                                        ASSERT( !conv.lgSearch );
                                        ASSERT( StatesElem[ipISO][nelem][0].Pop <=
                                                dense.xIonDense[nelem][nelem-ipISO]/
                                                (double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk ||
                                                /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
                                                iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);

                                        /* make sure that populations are valid, first check Pop2Ion  */
                                        if( StatesElem[ipISO][nelem][0].Pop*
                                                dense.xIonDense[nelem][nelem+1-ipISO] >
                                                dense.xIonDense[nelem][nelem-ipISO]*err_chk &&
                                                /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
                                                iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15)
                                        {
                                                fprintf(ioQQQ," DISASTER PopLo fnz %.2f ipISO %li nelem %li pop_ion_ov_neut %.2e 2pops %.6e %.6e solver:%s\n", 
                                                        fnzone,
                                                        ipISO , nelem ,
                                                        iso.pop_ion_ov_neut[ipISO][nelem] , 
                                                        StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO],
                                                        dense.xIonDense[nelem][nelem-ipISO]*err_chk ,
                                                        iso.chTypeAtomUsed[ipISO][nelem]);
                                        }
                                        ASSERT( 
                                                /*(conv.lgSearch || !conv.nPres2Ioniz) || */
                                                (StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO]<=
                                                        dense.xIonDense[nelem][nelem-ipISO]*err_chk) ||
                                                /*>>chng 06 jul 22, only do check if pops are within tolerance of double */
                                                iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);
#                                       undef ERR_CHK
                                }
                        }
                }
        }
#       endif

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, 
                        " radius_increment called; radius=%10.3e rinner=%10.3e DRAD=%10.3e drNext=%10.3e ROUTER=%10.3e DEPTH=%10.3e\n", 
                  radius.Radius, radius.rinner, radius.drad, radius.drNext, 
                  radius.router[iteration-1], radius.depth );
        }

        /* remember mean and largest errors on electron density */
        Error = fabs(dense.eden - dense.EdenTrue)/SDIV(dense.EdenTrue);
        if( Error > conv.BigEdenError )
        {
                conv.BigEdenError = (realnum)Error;
                dense.nzEdenBad = nzone;
        }
        conv.AverEdenError += (realnum)Error;

        /* remember mean and largest errors between heating and cooling */
        Error = fabs(thermal.ctot - thermal.htot) / thermal.ctot;
        conv.BigHeatCoolError = MAX2((realnum)Error , conv.BigHeatCoolError );
        conv.AverHeatCoolError += (realnum)Error;

        /* remember mean and largest pressure errors */
        Error = fabs(pressure.PresTotlCurr - pressure.PresTotlCorrect) / pressure.PresTotlCorrect;
        conv.BigPressError = MAX2((realnum)Error , conv.BigPressError );
        conv.AverPressError += (realnum)Error;

        /* integrate total mass over model */
        dense.xMassTotal += dense.xMassDensity * (realnum)(radius.drad_x_fillfac*radius.r1r0sq);

        /* check cooling time for this zone, remember longest */
        timesc.time_therm_long = MAX2(timesc.time_therm_long,1.5*dense.pden*BOLTZMANN*phycon.te/
          thermal.ctot);

        /* fmol is fraction of hydrogen in form of any molecule or ion */
        /* if N_H_MOLEC != 8, we probably need to change this line... */
        ASSERT( N_H_MOLEC == 8);
        fmol = (hmi.Hmolec[ipMHm] + 2.*(hmi.H2_total + hmi.Hmolec[ipMH2p]))/dense.gas_phase[ipHYDROGEN];

        /* remember the largest H2 fraction that occurs in the model */
        hmi.BiggestH2 = MAX2(hmi.BiggestH2,(realnum)fmol);

        /* largest fraction of atoms in molecules */
        for( i=0; i<mole.num_comole_calc; ++i )
        {
                if( dense.lgElmtOn[COmole[i]->nelem_hevmol] )
                {
                        realnum frac = COmole[i]->hevmol / SDIV(dense.gas_phase[COmole[i]->nelem_hevmol]);
                        frac *= (realnum) COmole[i]->nElem[COmole[i]->nelem_hevmol];
                        COmole[i]->xMoleFracMax = MAX2(frac,COmole[i]->xMoleFracMax);
                }
        }

        /* H 21 cm equilibrium timescale, H21cm returns H (not e) collisional 
         * deexcitation rate (not cs) */
        t = H21cm_H_atom( phycon.te )* dense.xIonDense[ipHYDROGEN][0] +
                /* >>chng 02 feb 14, add electron term as per discussion in */
                /* >>refer      H1      21cm    Liszt, H., 2001, A&A, 371, 698 */
                H21cm_electron( phycon.te )*dense.eden;

        /* only update time scale if t is significant */
        if( t > SMALLFLOAT )
                timesc.TimeH21cm = MAX2( 1./t, timesc.TimeH21cm );

        /* remember longest CO timescale */
        if( dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0] > SMALLFLOAT )
        {
                double timemin;
                double product = SDIV(dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0]);
                struct molecule *co_molecule;
                co_molecule = findspecies("CO");
                timemin = MIN2(timesc.AgeCOMoleDest[co_molecule->index] ,
                        timesc.AgeCOMoleDest[co_molecule->index] * co_molecule->hevmol/ product );
                /* this is rate CO is destroyed, equal to formation rate in equilibrium */
                timesc.BigCOMoleForm = MAX2(  timesc.BigCOMoleForm , timemin );
        }

        /* remember longest H2  destruction timescale timescale */
        timesc.time_H2_Dest_longest = MAX2(timesc.time_H2_Dest_longest, timesc.time_H2_Dest_here );

        /* remember longest H2 formation timescale timescale */
        timesc.time_H2_Form_longest = MAX2( timesc.time_H2_Form_longest , timesc.time_H2_Form_here );

        /* increment counter if this zone possibly thermally unstable
         * this flag was set in conv_temp_eden_ioniz.cpp, 
         * derivative of heating and cooling negative */
        if( thermal.lgUnstable )
                thermal.nUnstable += 1;

        /* remember Stromgren radius - where hydrogen ionization falls below half */
        if( !rfield.lgUSphON && dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN] > 0.49 )
        {
                rfield.rstrom = (realnum)radius.Radius;
                rfield.lgUSphON = true;
        }

        /* ratio of inner to outer radii, at this point
         * radius is the outer radius of this zone */
        DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/
          radius.Radius);

        if( trace.lgTrace && trace.lgConBug )
        {
                fprintf( ioQQQ, " Energy, flux, OTS:\n" );
                for( i=0; i < rfield.nflux; i++ )
                {
                        fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i], 
                          rfield.flux[i] + rfield.outlin[i] + rfield.ConInterOut[i], 
                          rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]);
                }
                fprintf( ioQQQ, "\n" );
        }

        /* diffuse continua factor
         * if lgSphere then all diffuse radiation is outward (COVRT=1)
         * lgSphere not set then COVRT=0.0 */

        /* begin sanity check */
        /* only do this test in debug mode */
#       if !defined(NDEBUG)
        lgFlxNeg = false;
        for( i=0; i < rfield.nflux; i++ )
        {
                if( rfield.flux[i] < 0. )
                {
                        fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" );
                        fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n", 
                          rfield.flux[i], rfield.anu[i], i );
                        lgFlxNeg = true;
                }
                if( rfield.otscon[i] < 0. )
                {
                        fprintf( ioQQQ, " radius_increment finds negative intensity in otscon.\n" );
                        fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n", 
                          rfield.otscon[i], rfield.anu[i], i );
                        lgFlxNeg = true;
                }
                if( opac.tmn[i] < 0. )
                {
                        fprintf( ioQQQ, " radius_increment finds negative tmn.\n" );
                        fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
                          opac.tmn[i], rfield.anu[i], i, rfield.chLineLabel[i] );
                        lgFlxNeg = true;
                }
                if( rfield.otslin[i] < 0. )
                {
                        fprintf( ioQQQ, " radius_increment finds negative intensity in otslin.\n" );
                        fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
                          rfield.otslin[i], rfield.anu[i], i, rfield.chLineLabel[i]  );
                        lgFlxNeg = true;
                }
                if( rfield.outlin[i] < 0. )
                {
                        fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" );
                        fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
                          rfield.outlin[i], rfield.anu[i], i, rfield.chLineLabel[i]  );
                        lgFlxNeg = true;
                }
                if( rfield.ConInterOut[i] < 0. )
                {
                        fprintf( ioQQQ, " radius_increment finds negative intensity in ConInterOut.\n" );
                        fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
                          rfield.ConInterOut[i], rfield.anu[i], i, rfield.chContLabel[i]  );
                        lgFlxNeg = true;
                }
                if( opac.opacity_abs[i] < 0. )
                {
                        opac.lgOpacNeg = true;
                        /* this sub will punch negative opacities on io unit,
                         * iff 'set negopc' command was given */
                        pnegopc();
                }
        }
        if( lgFlxNeg )
        {
                fprintf( ioQQQ, " Insanity has occurred, this is zone%4ld\n", 
                  nzone );
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }
        /*end sanity check*/
#       endif

        /* rfield.lgOpacityFine flag set false with no fine opacities command 
         * tests show that always evaluating this changes fast run of
         * pn_paris from 26.7 sec to 35.1 sec 
         * but must always update fine opacities since used for transmission */
        if( rfield.lgOpacityFine )
        {
                dradfac = (realnum)radius.drad_x_fillfac;
                /* increment the fine opacity array */
                for( i=0; i<rfield.nfine; ++i )
                {
                        rfield.fine_opt_depth[i] += 
                                rfield.fine_opac_zone[i]*dradfac;
                }

                /* evaluate fine opacity transmission coefficient for transmission
                 * through this zone. This is assumed to be a scattering opacity, 
                 * only do this if including scattering */
                if( opac.lgScatON )
                {
                        /* sum over coarse continuum */
                        for( i=0; i < rfield.nflux-1; i++ )
                        {
                                /* find transmission coefficient if lower and upper bounds of coarse continuum is within
                                 * boundaries of fine continuum 
                                 * unity is default */
                                if( rfield.ipnt_coarse_2_fine[i] && rfield.ipnt_coarse_2_fine[i+1] )
                                {
                                        rfield.trans_coef_zone[i] = 0.;
                                        for( j=rfield.ipnt_coarse_2_fine[i]; j<rfield.ipnt_coarse_2_fine[i+1]; ++j )
                                        {
                                                /* find transmission coefficient for this zone */
                                                rfield.trans_coef_zone[i] += 1.f / (1.f + rfield.fine_opac_zone[j]*dradfac );
                                        }
                                        /* first branch is normal case, where fine continuum is finer than
                                         * coarse continuum.  But, when end temp is very high, fine continuum is
                                         * very coarse, so may be just one cell, and following will not pass */
                                        if( rfield.ipnt_coarse_2_fine[i+1]>rfield.ipnt_coarse_2_fine[i] )
                                        {
                                                rfield.trans_coef_zone[i] /= (rfield.ipnt_coarse_2_fine[i+1]-rfield.ipnt_coarse_2_fine[i]);
                                        }
                                        else
                                        {
                                                /* in case where fine is coarser than coarse, just use firs cell */
                                                rfield.trans_coef_zone[i] = 
                                                        1.f / (1.f + rfield.fine_opac_zone[rfield.ipnt_coarse_2_fine[i]]*dradfac );
                                        }
                                }
                                else
                                {
                                        /* this branch, not within find continuum energy range - transmission unity */
                                        rfield.trans_coef_zone[i] = 1.f;
                                }

                                /* the total is just the current total times the transmission coefficient in this zone. */
                                rfield.trans_coef_total[i] *= rfield.trans_coef_zone[i];
                        }
                }
        }

        /* electron scattering opacity */
        escatc = 6.65e-25*dense.eden;

        /* this loop should not be to <= nflux since we not want to clobber the
         * continuum unit integration */
        relec = 0.;
        rforce = 0.;
        for( i=0; i < rfield.nflux; i++ )
        {
                /* sum total continuous optical depths */
                opac.TauAbsGeo[0][i] += (realnum)(opac.opacity_abs[i]*radius.drad_x_fillfac);
                opac.TauScatGeo[0][i] += (realnum)(opac.opacity_sct[i]*radius.drad_x_fillfac);

                /* following only optical depth to illuminated face */
                opac.TauAbsFace[i] += (realnum)(opac.opacity_abs[i]*radius.drad_x_fillfac);
                opac.TauScatFace[i] += (realnum)(opac.opacity_sct[i]*radius.drad_x_fillfac);

                /* these are total in inward direction, large if spherical */
                opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i];

                /* TMN is array of scale factors which account for attenuation
                 * of continuum across the zone (necessary to conserve energy
                 * at the 1% - 5% level.) sphere effects in
                 * drNext was set by NEXTDR and will be next dr */
                /* compute both total and Thomson scat rad acceleration */
                rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ 
                        rfield.ConInterOut[i])* rfield.anu[i]*(opac.opacity_abs[i] + 
                  opac.opacity_sct[i]);

                relec += ((rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ 
                        rfield.ConInterOut[i])* escatc)*rfield.anu[i];

                /* attenuation of flux by optical depths IN THIS ZONE 
                 * AngleIllumRadian is 1/COS(theta), is usually 1, reset with illuminate command,
                 * option for illumination of slab at an angle */
                opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac*
                        geometry.DirectionalCosin);

                /* e(-tau) in inward direction, up to illuminated face */
                opac.ExpmTau[i] *= (realnum)opac.ExpZone[i];

                /* e2(tau) in inward direction, up to illuminated face */
                opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
                ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. );

                /* on second and later iterations define outward E2 */
                if( iteration > 1 )
                {
                        /* e2 from current position to outer edge of shell */
                        realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
                        opac.E2TauAbsOut[i] = (realnum)e2( tau );
                        ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
                }

                /* DilutionHere is square of ratio of inner to outer radius */
                opfac = opac.ExpZone[i]*DilutionHere;
                ASSERT( opfac <= 1.0 );

                /* >>chng 07 may 22, break continuum into three parts */
                rfield.flux_beam_const[i] *= (realnum)opfac;
                rfield.flux_beam_time[i] *= (realnum)opfac;
                rfield.flux_isotropic[i] *= (realnum)opfac;
                rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
                        rfield.flux_isotropic[i];

                /* update SummedCon here since flux changed */
                rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];

                /* outward lines and continua */
                rfield.ConInterOut[i] *= (realnum)opfac;
                rfield.outlin[i] *= (realnum)opfac;
                rfield.outlin_noplot[i] *= (realnum)opfac;

                rfield.ConEmitOut[i] *= (realnum)opfac;
                rfield.ConEmitOut[i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i];

                /* set occupation numbers, first attenuated incident continuum */
                rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i];

                /* >>chng 00 oct 03, add diffuse continua */
                /* local diffuse continua */
                rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[nzone][i]*rfield.convoc[i];

                /* >>chng 05 feb 16, occupation number of outward diffuse continuum */
                rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[i]*rfield.convoc[i];
        }

        /* begin sanity check, compare total Lyman continuum optical depth 
         * with amount of extinction there */

        /* this is amount continuum attenuated to illuminated face, 
         * but only do test if flux positive, not counting scattering opacity,
         * and correction for spherical dilution not important */
        /* >>chng 02 dec 05, add test for small float, protect against models 
         * where we have gone below smallfloat, and so float is ragged */
        if( rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]>SMALLFLOAT &&
                (rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
                SDIV(rfield.flux_total_incident[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) ) > SMALLFLOAT && 
                !opac.lgScatON &&
                radius.depth/radius.Radius < 1e-4 )
        {
                /* ratio of current to incident continuum, converted to optical depth */
                /* >>chng 99 apr 23, this crashed on alpha due to underflow to zero in argy
                 * to log.  defended two ways - above checks that ratio of fluxes is large enough, 
                 * and here convert to double.
                 * error found by Peter van Hoof */
                double tau_effec = -log((double)rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
                        (double)opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
                        (double)rfield.flux_total_incident[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]);

                /* this is computed absorption optical depth to illuminated face */
                /* >>chng 06 jul 11, add geometry factor for illumination at an angle - bug fix - should have
                 * always been there */
                double tau_true = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]*geometry.DirectionalCosin;

                /* first test is relative error, second is to absolute error and comes
                 * in for very small tau, where differences are in the round-off */
                if( fabs( tau_effec - tau_true ) / t > 0.01 && 
                        /* for very small inner optical depths, the tmn correction is major,
                         * and this test is not precise */
                        fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) )
                {
                        /* in print below must add extra HI col den since this will later be
                         * incremented in RT_tau_inc */
                        fprintf( ioQQQ,
                                " PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n",
                                nzone, 
                                tau_effec , 
                                tau_true ,
                                6.327e-18*(dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac+colden.colden[ipCOL_H0]) );
                        TotalInsanity();
                }
        }
        /* end sanity check */

        /* do scattering opacity, intensity goes as 1/(1+tau) 
         * scattering opacity is turned off when sphere is set */
        if( opac.lgScatON )
        {
                for( i=0; i < rfield.nflux; i++ )
                {
                        /* assume half forward scattering
                         * opfac = 1. / ( 1. + dReff*0.5 * scatop(i) )
                         * >>chng 97 apr 25, remove .5 to get agreement with
                         * Lightman and White equation 11 in small epsilon limit,
                         * >>refer      continuum       RT      Lightman, A.P., & White, T.R. 1988, ApJ, 335, 57 */
                        opfac = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
                        ASSERT( opfac <= 1.0 );
                        rfield.flux_beam_const[i] *= (realnum)opfac;
                        rfield.flux_beam_time[i] *= (realnum)opfac;
                        rfield.flux_isotropic[i] *= (realnum)opfac;
                        rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
                                rfield.flux_isotropic[i];

                        rfield.ConInterOut[i] *= (realnum)opfac;
                        rfield.ConEmitOut[i] *= (realnum)opfac;
                        rfield.outlin[i] *= (realnum)opfac;
                        rfield.outlin_noplot[i] *= (realnum)opfac;
                }
        }

        /* now do slight reshuffling of energy due to Compton scattering */
        cmshft();

        /* remember the largest value */
        wind.AccelMax = (realnum)MAX2(wind.AccelMax,wind.AccelTot);

        /* force multiplier; relec can be zero for very low densities - so use double
         * form of safe_div - fmul is of order unity - wind.AccelLine and wind.AccelCont 
         * are both floats to will underflow long before relec will - fmul is only used
         * in output, not in any physics */
        relec *= (EN1RYD/SPEEDLIGHT/dense.xMassDensity);
        if( relec > SMALLFLOAT )
                wind.fmul = (realnum)( (wind.AccelLine + wind.AccelCont) / relec);
        else
                wind.fmul = 0.;

        /* keep track of average acceleration */
        wind.AccelAver += wind.AccelTot*(realnum)radius.drad_x_fillfac;
        wind.acldr += (realnum)radius.drad_x_fillfac;

        /* following is integral of radiative force */
        pressure.pinzon = (realnum)(wind.AccelTot*dense.xMassDensity*radius.drad_x_fillfac*geometry.DirectionalCosin);
        /*fprintf(ioQQQ," debuggg pinzon %.2f %.2e %.2e %.2e\n", 
                fnzone,pressure.pinzon,dense.xMassDensity,wind.AccelTot);*/
        pressure.PresInteg += pressure.pinzon;

        /* sound is sound travel time, sqrt term is sound speed */
        timesc.sound_speed_isothermal = sqrt(pressure.PresGasCurr/dense.xMassDensity);
        /* adiabatic sound speed assuming mono-atomic gas - gamma is 5/3*/
        timesc.sound_speed_adiabatic = sqrt(5.*pressure.PresGasCurr/(3.*dense.xMassDensity) );
        timesc.sound += radius.drad_x_fillfac / timesc.sound_speed_isothermal;

        /* attenuate neutrons if they are present */
        if( hextra.lgNeutrnHeatOn )
        {
                /* correct for optical depth effects */
                hextra.totneu *= (realnum)sexp(hextra.CrsSecNeutron*dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac*geometry.DirectionalCosin);
                /* correct for spherical effects */
                hextra.totneu *= (realnum)DilutionHere;
        }

        /* following radiation factors are extinguished by 1/r**2 and e- opacity
         * also bound electrons */

        /* do all emergent spectrum from illuminated face if model is NOT spherical */
        if( !geometry.lgSphere )
        {
                double Reflec_Diffuse_Cont;

                /* >>chng 01 jul 14, from lower limit of 0 to plasma frequency -
                 * never should have added diffuse emission from below the plasma frequency */
                for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
                {
                        if( opac.TauAbsGeo[0][i] < 30. )
                        {
                                /* ConEmitLocal is diffuse emission per unit vol, fill factor
                                 * the 1/2 comes from isotropic emission */
                                Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.*
                                        radius.drad_x_fillfac * opac.E2TauAbsFace[i]*radius.r1r0sq;

                                /* ConEmitReflec - reflected diffuse continuum */
                                rfield.ConEmitReflec[i] += (realnum)(Reflec_Diffuse_Cont);

                                /* the reflected part of the incident continuum */
                                rfield.ConRefIncid[i] += (realnum)(rfield.flux[i]*opac.opacity_sct[i]*
                                  radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);

                                /* reflected line emission */
                                rfield.reflin[i] += (realnum)(rfield.outlin[i]*opac.opacity_sct[i]*
                                  radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
                        }
                }
        }

        /* following is general method to find means weighted by various functions
         * called in IterStart to initialize to zero, called here to put in numbers
         * results will be weighted by radius and volume
         * this is the only place things must be entered to create averages 
         * code is in mean.c */
        aver("zone",1.,1.,"          ");
        aver("doit",phycon.te,1.,"    Te    ");
        aver("doit",phycon.te,dense.eden,"  Te(Ne)  ");
        aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHYDROGEN][1]," Te(NeNp) ");
        aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][1]  ," Te(NeHe+)");
        aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][2]  ,"Te(NeHe2+)");
        aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][1]  ," Te(NeO+) " );
        aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][2]  ," Te(NeO2+)");
        /*aver("doit",phycon.te,hmi.Hmolec[ipMH2g],"  Te(H2)  ");*/
        aver("doit",phycon.te,hmi.H2_total,"  Te(H2)  ");
        aver("doit",dense.gas_phase[ipHYDROGEN],1.,"   N(H)   ");
        aver("doit",dense.eden,dense.xIonDense[ipOXYGEN][2],"  Ne(O2+) ");
        aver("doit",dense.eden,dense.xIonDense[ipHYDROGEN][1],"  Ne(Np)  ");

        /* save information about structure of model, now used to get t^2 */
        /* max because if program aborts during search phase, will get to here
         * with nzone = -1 */
        nzone_minus_1 = MAX2( 0, nzone-1 );
        /* this is number of struc.xx zones with valid data */
        struc.nzone = nzone_minus_1+1;
        ASSERT(nzone_minus_1>=0 && nzone_minus_1 < struc.nzlim );
        struc.testr[nzone_minus_1] = (realnum)phycon.te;
        /* number of particles per unit vol */
        struc.DenParticles[nzone_minus_1] = dense.pden;
        /* >>chng 02 May 2001 rjrw: add hden for dilution */
        struc.hden[nzone_minus_1] = (realnum)dense.gas_phase[ipHYDROGEN];
        /* total grams per unit vol */
        struc.DenMass[nzone_minus_1] = dense.xMassDensity;
        struc.heatstr[nzone_minus_1] = thermal.htot;
        struc.coolstr[nzone_minus_1] = thermal.ctot;
        struc.volstr[nzone_minus_1] = (realnum)radius.dVeff;
        struc.drad[nzone_minus_1] = (realnum)radius.drad;
        struc.drad_x_fillfac[nzone_minus_1] = (realnum)radius.drad_x_fillfac;
        struc.histr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][0];
        struc.hiistr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][1];
        struc.ednstr[nzone_minus_1] = (realnum)dense.eden;
        struc.o3str[nzone_minus_1] = dense.xIonDense[ipOXYGEN][2];
        struc.pressure[nzone_minus_1] = (realnum)pressure.PresTotlCurr;
        struc.windv[nzone_minus_1] = (realnum)wind.windv;
        struc.AccelTot[nzone_minus_1] = wind.AccelTot;
        struc.AccelGravity[nzone_minus_1] = wind.AccelGravity;
        struc.pres_radiation_lines_curr[nzone_minus_1] = (realnum)pressure.pres_radiation_lines_curr;
        struc.GasPressure[nzone_minus_1] = (realnum)pressure.PresGasCurr;
        struc.depth[nzone_minus_1] = (realnum)radius.depth;
        /* save absorption optical depth from illuminated face to current position */
        struc.xLyman_depth[nzone_minus_1] = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]];
        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
        {
                struc.gas_phase[nelem][nzone_minus_1] = dense.gas_phase[nelem];
                for( ion=0; ion<nelem+2; ++ion )
                {
                        struc.xIonDense[nelem][ion][nzone_minus_1] = dense.xIonDense[nelem][ion];
                }
        }

        /* the hydrogen molecules */
        for(mol=0;mol<N_H_MOLEC;mol++) 
        {
                struc.H2_molec[mol][nzone_minus_1] = hmi.Hmolec[mol];
        }

        /* the heavy element molecules */
        for(mol=0;mol<mole.num_comole_calc;mol++) 
        {
                struc.CO_molec[mol][nzone_minus_1] = COmole[mol]->hevmol;
        }

        colden.dlnenp += dense.eden*(double)(dense.xIonDense[ipHYDROGEN][1])*radius.drad_x_fillfac;
        colden.dlnenHep += dense.eden*(double)(dense.xIonDense[ipHELIUM][1])*radius.drad_x_fillfac;
        colden.dlnenHepp += dense.eden*(double)(dense.xIonDense[ipHELIUM][2])*radius.drad_x_fillfac;
        colden.dlnenCp += dense.eden*(double)(dense.xIonDense[ipCARBON][1])*radius.drad_x_fillfac;

        /* >>chng 05 jan 09, add integral of n(H0) / Tspin */
        colden.H0_ov_Tspin += (double)(dense.xIonDense[ipHYDROGEN][0]) / hyperfine.Tspin21cm*radius.drad_x_fillfac;

        /* >>chng 05 Mar 07, add integral of n(OH) / Tspin */
        colden.OH_ov_Tspin += (double)(findspecies("OH")->hevmol) / phycon.te*radius.drad_x_fillfac;

        /* this is Lya excitation temperature */
        hydro.TexcLya = (realnum)TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );

        /* count number of times Lya excitation temp hotter than gas */
        if( hydro.TexcLya > phycon.te )
        {
                hydro.nLyaHot += 1;
                if( hydro.TexcLya > hydro.TLyaMax )
                {
                        hydro.TLyaMax = hydro.TexcLya;
                        hydro.TeLyaMax = (realnum)phycon.te;
                        hydro.nZTLaMax = nzone;
                }
        }

        /* column densities in various species */
        colden.colden[ipCOL_HTOT] += (realnum)(dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac);
        ASSERT( colden.colden[ipCOL_HTOT] >SMALLFLOAT );
        colden.colden[ipCOL_HMIN] += hmi.Hmolec[ipMHm]*(realnum)radius.drad_x_fillfac;
        /* >>chng 02 sep 20, from htwo to H2_total */
        /* >>chng 05 mar 14, rather than H2_total, give H2g and H2s */
        colden.colden[ipCOL_H2g] += hmi.Hmolec[ipMH2g]*(realnum)radius.drad_x_fillfac;
        colden.colden[ipCOL_H2s] += hmi.Hmolec[ipMH2s]*(realnum)radius.drad_x_fillfac;
        /* this is a special form of column density - should be proportional to total shielding */
        colden.coldenH2_ov_vel += hmi.H2_total*(realnum)radius.drad_x_fillfac / DoppVel.doppler[LIMELM+2];
        colden.colden[ipCOL_HeHp] += hmi.Hmolec[ipMHeHp]*(realnum)radius.drad_x_fillfac;
        colden.colden[ipCOL_H2p] += hmi.Hmolec[ipMH2p]*(realnum)radius.drad_x_fillfac;
        colden.colden[ipCOL_H3p] += hmi.Hmolec[ipMH3p]*(realnum)radius.drad_x_fillfac;
        colden.colden[ipCOL_Hp] += dense.xIonDense[ipHYDROGEN][1]*(realnum)radius.drad_x_fillfac;
        colden.colden[ipCOL_H0] += dense.xIonDense[ipHYDROGEN][0]*(realnum)radius.drad_x_fillfac;

        colden.colden[ipCOL_elec] += (realnum)(dense.eden*radius.drad_x_fillfac);
        /* He I t2S column density*/
        colden.He123S += dense.xIonDense[ipHELIUM][1]*
                StatesElem[ipHE_LIKE][ipHELIUM][ipHe2s3S].Pop*(realnum)radius.drad_x_fillfac;

        /* the ortho and para column densities */
        h2.ortho_colden += h2.ortho_density*radius.drad_x_fillfac;
        h2.para_colden += h2.para_density*radius.drad_x_fillfac;
        ASSERT( fabs( h2.ortho_density + h2.para_density - hmi.H2_total ) / SDIV( hmi.H2_total ) < 1e-4 );
        /*fprintf(ioQQQ,"DEBUG ortho para\t%.3e\t%.3e\ttot\t%.3e\t or pa colden\t%.3e\t%.3e\n", 
                h2.ortho_density, h2.para_density,hmi.H2_total,
                h2.ortho_colden , h2.para_colden);*/

        /*>>chng 27mar, GS, Column density of F=0 and F=1 levels of H0*/
        colden.H0_21cm_upper += (HFLines[0].Hi->Pop*radius.drad_x_fillfac);
    colden.H0_21cm_lower +=(HFLines[0].Lo->Pop*radius.drad_x_fillfac);
        /*fprintf(ioQQQ,"DEBUG pophi-poplo\t%.3e\t%.3e\radius\t%.3e\t col_hi\t%.3e\t%.3e\n", 
                HFLines[0].PopHi, HFLines[0].PopLo, radius.drad_x_fillfac,
                 HFLines[0].PopHi*radius.drad_x_fillfac,colden.H0_21cm_upper );*/

        /* the CII and SiII atoms are resolved */
        for( i=0; i < 5; i++ )
        {
                /* pops and column density for C II atom */
                colden.C2Colden[i] += colden.C2Pops[i]*(realnum)radius.drad_x_fillfac;
                /* pops and column density for SiII atom */
                colden.Si2Colden[i] += colden.Si2Pops[i]*(realnum)radius.drad_x_fillfac;
        }
        for( i=0; i < 3; i++ )
        {
                /* pops and column density for CI atom */
                colden.C1Colden[i] += colden.C1Pops[i]*(realnum)radius.drad_x_fillfac;
                /* pops and column density for OI atom */
                colden.O1Colden[i] += colden.O1Pops[i]*(realnum)radius.drad_x_fillfac;
        }
        for( i=0; i < 4; i++ )
        {
                /* pops and column density for CIII atom */
                colden.C3Colden[i] += colden.C3Pops[i]*(realnum)radius.drad_x_fillfac;
        }

        /* now add total molecular column densities */
        molcol("ADD ",ioQQQ);

        /* increment forming the mean ionization and temperature */
        MeanInc();

        /*-----------------------------------------------------------------------*/

        /* calculate average atomic weight per hydrogen of the plasma */
        avWeight = 0.;
        for( nelem=0; nelem < LIMELM; nelem++ ) 
        {
                avWeight += dense.gas_phase[nelem]*dense.AtomicWeight[nelem];
        }
        avWeight /= dense.gas_phase[ipHYDROGEN];

        /* compute some average grain properties */
        rfield.opac_mag_B_point = 0.;
        rfield.opac_mag_V_point = 0.;
        rfield.opac_mag_B_extended = 0.;
        rfield.opac_mag_V_extended = 0.;
        long int nd;
        for( nd=0; nd < gv.nBin; nd++ )
        {

                /* this is total extinction in magnitudes at V and B, for a point source 
                 * total absorption and scattering,
                 * does not discount forward scattering to be similar to stellar extinction
                 * measurements made within ism */
                rfield.extin_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
                        radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
                rfield.opac_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
                        dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

                rfield.extin_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
                        radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

                rfield.opac_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
                        dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

                /* this is total extinction in magnitudes at V and B, for an extended source 
                 * total absorption and scattering,
                 * DOES discount forward scattering to apply for extended source like Orion */
                rfield.extin_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund*
                        radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
                rfield.opac_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund*
                        dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

                rfield.extin_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund*
                        radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
                rfield.opac_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
                        gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund*
                        dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

                gv.bin[nd]->avdust += gv.bin[nd]->tedust*(realnum)radius.drad_x_fillfac;
                gv.bin[nd]->avdft += gv.bin[nd]->DustDftVel*(realnum)radius.drad_x_fillfac;
                gv.bin[nd]->avdpot += (realnum)(gv.bin[nd]->dstpot*EVRYD*radius.drad_x_fillfac);
                gv.bin[nd]->avDGRatio += (realnum)(gv.bin[nd]->dustp[1]*gv.bin[nd]->dustp[2]*
                        gv.bin[nd]->dustp[3]*gv.bin[nd]->dustp[4]*gv.bin[nd]->dstAbund/avWeight*
                        radius.drad_x_fillfac);
        }

        /* there are some quantities needed to calculation the Jeans mass and radius */
        colden.TotMassColl += dense.xMassDensity*(realnum)radius.drad_x_fillfac;
        colden.tmas += (realnum)phycon.te*dense.xMassDensity*(realnum)radius.drad_x_fillfac;
        colden.wmas += dense.wmole*dense.xMassDensity*(realnum)radius.drad_x_fillfac;

        /* now find minimum Jeans length and mass; length in cm */
        rjeans = 7.79637 + (phycon.alogte - log10(dense.wmole) - log10(dense.xMassDensity*
          geometry.FillFac))/2.;

        /* minimum Jeans mass in gm */
        ajmass = 3.*(rjeans - 0.30103) + log10(4.*PI/3.*dense.xMassDensity*
          geometry.FillFac);

        /* now remember smallest */
        colden.rjnmin = MIN2(colden.rjnmin,(realnum)rjeans);
        colden.ajmmin = MIN2(colden.ajmmin,(realnum)ajmass);

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, " radius_increment returns\n" );
        }
        return;
}

#if !defined(NDEBUG)
/*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */
STATIC void pnegopc(void)
{
        long int i;
        FILE *ioFile;

        DEBUG_ENTRY( "pnegopc()" );

        if( opac.lgNegOpacIO )
        {
                /* option to punch negative opacities */
                ioFile = open_data( "negopc.txt", "w", AS_LOCAL_ONLY );
                for( i=0; i < rfield.nflux; i++ )
                {
                        fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i], 
                          opac.opacity_abs[i] );
                }
                fclose( ioFile);
        }
        return;
}
#endif
/* Note - when this file is compiled in fast optimized mode,
 * a good compiler will complain that the routine pnegopc is a
 * static, local routine, but that it has not been used.  It is
 * indeed used , but only when NDEBUG is not set, so it is only
 * used when the code is compiled in debug mode.  So this is
 * not a problem. */


/*cmshft - so Compton scattering shift of spectrum */
STATIC void cmshft(void)
{
        long int i;

        DEBUG_ENTRY( "cmshft()" );

        /* first check whether Compton scattering is in as heat/cool */
        if( rfield.comoff == 0. )
        { 
                return;
        }

        if( rfield.comoff != 0. )
        { 
                return;
        }

        /* do reshuffle */
        for( i=0; i < rfield.nflux; i++ )
        {
                continue;
                /* watch this space for some really great code!!!!
                 * COMUP needs factor of TE to be Compton cooling */
        }
        return;
}

---