/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * 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 */
/*TODO	2	- this routine is very important since it sets the pace for the calculation,
 * and directly affects the convergence of the code.  Most of the logic is very old and
 * messy.  
 * 1) make sure all test cases have punch dr
 * 2) cat all these reasons together into one file and sort on the reason
 * 3) discover what logic is the main pacesetter for the code
 * 4) which are never triggered and so can be removed
 */
#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 "tabden.h"
#include "fabden.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);
/*lint -e777 floating tests for equality are ok and frequent here */

/*radius_next use adaptive logic to find next zone thickness 
 * return 0 if ok, 1 for abort */
int radius_next(void)
{
	char chLbl[11];
	int 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, 
	  dVeldRad,
	  fac2, 
	  factor, 
	  freqm, 
	  grfreqm=0., 
	  gropacm=0., 
	  hdnew, 
	  opacm, 
	  OldDR ,
	  winddr, 
	  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;
	int lgFirstCall;

#	ifdef DEBUG_FUN
	fputs( "<+>radius_next()\n", debug_fp );
#	endif


	/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	 *
	 * 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( n=0; n<N_ELEM_CHECK; ++n )*/
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
		{
			/*nelem = n_elem_check[n];*/
			if( dense.lgElmtOn[nelem] )
			{
				float change;
				/* this is the limit to the ionization we will check -
				 * also used in prt_comment to check on whether oscillations occurred */
				float frac_limit;
				if( nelem<=ipHELIUM )
				{
					frac_limit = 1e-4f;
					change = CHANGE_ION_HHE;
				}
				else
				{
					/* >>chng 04 feb 28, from limit to this / 2,
					 * 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 )
				{
					float abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
					float 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) */
						/*float abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);*/
						float abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]);
						float 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.lgTSetOn )
	{
		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 */
	if( !thermal.lgTSetOn )
	{
		/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
		/* 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.AngleIllum);
	}
	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 = vfmax(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen,FEND);

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

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

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

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

	else
	{
		fprintf( ioQQQ, " insanity following ChkRate\n" );
		ShowMe();
		puts( "[Stop in radius_next]" );
		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 */
	/*>>chng 04 may 25,
	 *  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 fractional accel over length */
		dVeldRad = fabs(wind.windv-OldWindVelocity)/
			MAX2(v,0.1*timesc.sound_speed_isothermal)/radius.drad;

		if( 1.1*dVeldRad*radius.drad > 0.03  )
		{
			/* dVeldRad is D(vel)/vel / DRAD, computed in convpres */
			winddr = 0.03/dVeldRad;
		}
		else
		{
			winddr = 1.1*radius.drad;
		}

		/* >>chng 02 nov 05, add dr from advective term,
		 * the 1/500 came from looking at one set of structure plots */
		if( dynamics.lgAdvection )
		{
			/* >>chng 02 dec 11, from 5 to 20 */
			winddr = MIN2( winddr , dynamics.dRad / 20. );
			/*>>chng 04 oct 06, set dVeldRad to dynamics.dRad since dVeldRad is printed as part
			 * of reason for choosing this criteria, want to reflect valid reason. */
			dVeldRad = dynamics.dRad;
		}
	}
	else
	{
		winddr = BigRadius;
		dVeldRad = 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 );
			puts( "[Stop in radius_next]" );
			cdEXIT(EXIT_FAILURE);
		}
		factor = radius.glbden*pow(radius.glbrad/radius.glbdst,radius.glbpow);
		fac2 = radius.glbden*pow(radius.glbrad/(radius.glbdst - (float)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 = fabden(radius.Radius+radius.drad,radius.depth+
			  radius.drad);
		}
		else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
		{
			hdnew = tabden(radius.Radius+radius.drad,radius.depth+
			  radius.drad);
		}
		else
		{
			fprintf( ioQQQ, " dlw insanity in radius_next\n" );
			puts( "[Stop in radius_next]" );
			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(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.AngleIllum);
	}
	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].TauIn;
				TauDTau = TauLines[ipStrong].PopOpc * TauLines[ipStrong].opacity / 
					DoppVel.doppler[TauLines[ipStrong].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].TauIn;
				TauDTau = TauLine2[ipStrong].PopOpc * TauLine2[ipStrong].opacity / 
					DoppVel.doppler[TauLine2[ipStrong].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].TauIn;
				TauDTau = C12O16Rotate[ipStrong].PopOpc * C12O16Rotate[ipStrong].opacity / 
					DoppVel.doppler[C12O16Rotate[ipStrong].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].TauIn;
				TauDTau = C13O16Rotate[ipStrong].PopOpc * C13O16Rotate[ipStrong].opacity / 
					DoppVel.doppler[C13O16Rotate[ipStrong].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].TauIn;
				TauDTau = HFLines[ipStrong].PopOpc * HFLines[ipStrong].opacity / 
					DoppVel.doppler[HFLines[ipStrong].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.lgTSetOn )
	{
		/* 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 ipMole;
		Old_H2_abund = struc.Hmolec[ipMH2g][nzone-3] + struc.Hmolec[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(ipMole=0; ipMole<NUM_HEAVY_MOLEC; ++ipMole)
		{
			float abund,
				abund_limit;
			if( !dense.lgElmtOn[co.nelem_hevmol[ipMole]] )
				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 = co.hevmol[ipMole]/dense.gas_phase[co.nelem_hevmol[ipMole]];
			/* 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( co.lgGas_Phase[ipMole] )
			{
				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.COmolec[ipMole][nzone-3]>SMALLFLOAT )
				{
					relative_change_denom = MIN2( co.hevmol[ipMole] , struc.COmolec[ipMole][nzone-3] );
				}
				else
				{
					relative_change_denom = co.hevmol[ipMole];
				}
				/* the relative change in the abundance */
				relative_change = 
					/*fabs( co.hevmol[ipMole] - struc.COmolec[ipMole][nzone-3] ) / co.hevmol[ipMole];*/
					/* >>chng 05 dec 08, use smaller abundance */
					fabs( co.hevmol[ipMole] - struc.COmolec[ipMole][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 = ipMole;
					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.AngleIllum 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 */
		float errorHC = POW2(conv.HeatCoolRelErrorAllowed);
		float errorNe = (float)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;
	}

	/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
	/* 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( (co.hevmol[ipCO])*geometry.FillFac);
		DepthToGo = MIN2(DepthToGo ,
			(StopCalc.col_monoxco-co.hevcol[ipCO]) / 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.lgTSetOn ) 
		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 = vfmin(drmax,drInter,drLineHeat,winddr,drFluc,GlobDr,
	  DrGrainHeat,dr_change_heavy,FEND);
	radius.drNext = vfmin(radius.drNext,SpecDr,drFail,drSphere,radius.sdrmax,
	  dRTaue,FEND);
	radius.drNext = vfmin(radius.drNext,drDest,drdTdStep,drdHdStep,
	  drConPres,drTab,drSolomon_BigH2,drLeiden_hack,FEND);
	radius.drNext = vfmin(radius.drNext,drHion,drDepth, dr_mole_abund , AV_to_go ,
	  drEfrac , drHMase , drThickness , drHe1ovHe2 , drH2_heat_cool , drH2_abund , FEND);

	/*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;
	}

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

	/*TODO	1	rm logic detecting temp and pres failure - these normally occur
	 * near troubled parts of cooling curve */
	/* 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;
	if( nzone >=5 )
	{
		/* >>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 );
		puts( "[Stop in radius_next]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* set flag if dr set by maser */
	if( 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( 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].pump, 
				  TauLines[ipStrong].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].pump, 
				  TauLine2[ipStrong].Pesc );
			}
			else
			{
				fprintf( ioQQQ, " insanity pr line heat\n" );
				puts( "[Stop in radius_next]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

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

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

		else if( radius.drNext == dr_change_heavy )
		{
			fprintf( punch.ipDRout, 
				"change in ionization, element %s ion %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( radius.drNext == drThickness )
		{
			fprintf( punch.ipDRout, "depth to go\n");
		}

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

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

		else if( 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( 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( 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( 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( 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 , co.chLab[mole_dr_change]);
		}

		else if( 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( radius.drNext == drHion )
		{
			fprintf( punch.ipDRout, "change in H ionization fm to%11.3e%11.3e\n", 
			  SaveOHIIoHI, OHIIoHI );
		}

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

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

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

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

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

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

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

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

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

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

		else if( 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( radius.drNext == drInter )
		{
			fprintf( punch.ipDRout, 
				"cont inter nu=%10.2e opac=%10.2e\n", 
			  freqm, opacm );
		}

		else if( radius.drNext == DrGrainHeat )
		{

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

		else if( radius.drNext == winddr )
		{
			fprintf( punch.ipDRout, 
				"Wind, dVeldRad=%10.3e\n", 
			   dVeldRad );
		}

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

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

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

		else if( radius.drNext == OldDR )
		{
			fprintf( punch.ipDRout, "old DR.\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 );
			puts( "[Stop in radius_next]" );
			cdEXIT(EXIT_FAILURE);
		}

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

	/* 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,
	 * not at carbon +-0 i front
	 * not flucuations 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 only set wih the "set dr" command */
	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 */
		/* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth
		 * avoid false trigger across thermal fronts 
		 * add * dense.gas_phase */
		if( radius.drNext * 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\n\n", 
				nzone, 
				iteration);
			return(1);
		}
	}

	/* factor to allow for slop in floating numbers */
#	define	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 */
	radius.drNext = MIN2(radius.drNext,(radius.router[iteration-1]-
	  radius.depth)*Z);

	/* this means outer limit exceeded */
	if( radius.drNext < 0. )
	{
		radius.lgDrNeg = TRUE;
	}
	else
	{
		radius.lgDrNeg = FALSE;
	}

	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 );
	}

#	ifdef DEBUG_FUN
	fputs( " <->radius_next()\n", debug_fp );
#	endif

	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;

#	ifdef DEBUG_FUN
	fputs( "<+>ContRate()\n", debug_fp );
#	endif

	/* 
	 * 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 truned off, use hydrogen balmer continuum */
		ipHi = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]-1;
	}

	for( i=1; 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;
	}

	/* 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;
	}

	/* >>chng 05 aug 03, i dont see what this logic is trying to do - remove it
	 * it must be incorrect */
#	if 0
	if( Rate_max_Hion > Rate_max_nonIonizing )
	{
	}
	else
	{
		*opacm = Opac_nonIonizing;
		*freqm = Freq_nonIonizing;
	}
#	endif

	{
		/* following should be set true to print contributors */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		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. );

#	ifdef DEBUG_FUN
	fputs( " <->ContRate()\n", debug_fp );
#	endif
	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;

#	ifdef DEBUG_FUN
	fputs( "<+>GrainRateDr()\n", debug_fp );
#	endif

	/* 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.;
	}

#	ifdef DEBUG_FUN
	fputs( " <->GrainRateDr()\n", debug_fp );
#	endif
	return;
}
#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];

#	ifdef DEBUG_FUN
	fputs( "<+>ChkRate()\n", debug_fp );
#	endif

	/* return if this element is not turned on */
	if( !dense.lgElmtOn[nelem] )
	{
		*dDestRate = 1e-3;
		*DestRateOld = 1e-3;
		*DestRateNew = 1e-3;
		*istage = 0;
		
#		ifdef DEBUG_FUN
		fputs( " <->ChkRate()\n", debug_fp );
#		endif
		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] );
					puts( "[Stop in chkrate]" );
					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];
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->ChkRate()\n", debug_fp );
#	endif
	return;
}
#endif
/*lint +e777 floating tests for equality are ok and frequent here */
