/* This file is part of Cloudy and is copyright (C) 1978-2004 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 nextdr to find energy of maximum continuum-gas interaction */
/*GrainRateDr called by nextdr 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 "co.h"
#include "opacity.h"
#include "dense.h"
#include "heavy.h"
#include "grainvar.h"
#include "elementnames.h"
#include "converge.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 nextdr 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 nextdr to find energy of maximum continuum-gas interaction */
static void ContRate(double *freqm, 
  double *opacm);

/*GrainRateDr called by nextdr 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 , limit , k;

	double thickness_total , drThickness , DepthToGo ;
	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, 
	  drLineHeat, 
	  drPressure, 
	  drSphere, 
	  drTab, 
	  drdHdStep, 
	  drdTdStep, 
	  drmax, 
	  dVeldRad,
	  fac2, 
	  factor, 
	  freqm, 
	  grfreqm, 
	  gropacm, 
	  hdnew, 
	  opacm, 
	  OldDR ,
	  winddr, 
	  x,
	  dr_ne_oscil,
	  dr_te_oscil;

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

	/* >>chng 04 feb 15, kill this block - not used */
#	if 0
	/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
	/* check on how much outward flux is being added relative to
	 * attenuated incident beam
	 *  1    '' radius_next keys from ratio out2 tot continuua'',e9.1)') */
	if( (rfield.SumIncCon > 0. && rfield.SumOutCon > 0.) && 
	  !thermal.lgTSetOn )
	{
		/* evaluated in metdif */
		Out2Tot = rfield.SumOutCon/rfield.SumIncCon;
		/* Out2HTot = SumHOutCon / SumHIncCon
		 * try to keep outward frac less than 5% of total */
		drOutFrac = radius.drad*MAX2(0.7,0.05/MAX2(1e-10,Out2Tot));
	}
	else
	{
		Out2Tot = FLT_MAX;
		drOutFrac = BigRadius;
	}
#	endif

	/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
	/* 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 */
			drConPres = 0.05*pressure.PresTotlInit/(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 nextdr]" );
		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 nextdr]" );
			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 nextdr]" );
			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);
			drH2_heat_cool = radius.drad*MAX2(0.2,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 */
	drH2_abund = BigRadius;
	dr_mole_abund = BigRadius;
	mole_dr_change = -1;
	if( nzone>=4 )
	{
		/* first do H2 abundance */
		double Old_H2_abund,
			Older_H2_abund;
		int mole;
		Old_H2_abund = struc.Hmolec[ipMH2g][nzone-3] + struc.Hmolec[ipMH2s][nzone-3];
		Older_H2_abund = struc.Hmolec[ipMH2g][nzone-4] + struc.Hmolec[ipMH2s][nzone-4];
		/* >>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;
			/* >>chng 04 feb 27, add logic on existance of big H2 molecule*/
			if( h2.lgH2ON  && hmi.lgBigH2_evaluated && 
				/* >>chng 04 mar 13, do not make dr very small if H2 oscillating */
				(Old_H2_abund-hmi.H2_total)*(Older_H2_abund-Old_H2_abund)>0)
			{
				if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-3 )
				{
					/* lots of H2 */
					fac = 0.05;
				}
				else
					fac = 0.1;
			}
			/* this is percentage change in H2 density */
			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 */
			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(mole=0; mole<NUM_HEAVY_MOLEC; ++mole)
		{
			float abund;
			if( !dense.lgElmtOn[co.nelem_hevmol[mole]] )
				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[mole]/dense.gas_phase[co.nelem_hevmol[mole]];
			if( abund > 0.001 )
			{
				double drmole_one;
				/* the relative change in the abundance */
				double relative_change = 
					fabs( co.hevmol[mole] - struc.COmolec[mole][nzone-3] ) / co.hevmol[mole];
				/* >>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 */
				relative_change = SDIV(relative_change);
				drmole_one = radius.drad*MAX2(0.2,0.035/relative_change );
				/* final dr will be the smallest we encounter */
				if( drmole_one < dr_mole_abund )
				{
					/* this is the dr set used to pace changes in molecules */
					dr_mole_abund = drmole_one;
					mole_dr_change = mole;
					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 */
		drmax = 4.*radius.drad/geometry.AngleIllum;
	}
	else
	{
		drmax = 1.3*radius.drad/geometry.AngleIllum;
	}

	/* 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 = (float)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;
			}
		}
	}

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

	/* is current pressure ok? */
	if( conv.lgConvPres )
	{
		drPressure = BigRadius;
	}
	/* >>chng 04 may 25, add following test after nzone > 0 */
	else if( nzone > 0 && !pressure.lgSonicPoint && !pressure.lgStrongDLimbo )
	{
		/* WJH 22 May 2004: if we are trying to get through the sonic point, then 
		 * we may get lots of pressure failures. We would like to soldier on regardless */

		drPressure = radius.drad/2.;
	}
	else
	{
		drPressure = BigRadius;
	}

	/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
	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 mole 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[ipCHII]) / 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[ipCOLH2]) / 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[ipCOLH2]+dense.xIonDense[ipHYDROGEN][1])) / coleff  );
		thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / 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[ipCHI]) / 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);
	}

	/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
	/* 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,
		" 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,
	  drPressure,dRTaue,FEND);
	radius.drNext = vfmin(radius.drNext,drDest,drdTdStep,drdHdStep,
	  drConPres,drTab,drSolomon_BigH2,FEND);
	radius.drNext = vfmin(radius.drNext,drHion,drDepth, dr_mole_abund ,
	  drEfrac , drHMase , drThickness , drHe1ovHe2 , drH2_heat_cool , drH2_abund , dr_ne_oscil,
	  dr_te_oscil , FEND);

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

	/* 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%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, drPressure, dRTaue, 
		  OldH2Abund, drDepth );
		puts( "[Stop in nextdr]" );
		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 nextdr]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( radius.drNext == drDepth )
		{
			fprintf( punch.ipDRout, "relative depth\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 == 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 == dr_ne_oscil )
		{
			/* radius_next keys from electron density oscillation*/
			fprintf( punch.ipDRout, "electron density oscillation\n");
		}

		else if( radius.drNext == dr_te_oscil )
		{
			/* radius_next keys from electron density oscillation*/
			fprintf( punch.ipDRout, "temperature oscillation\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 );
		}
#		if 0
		else if( radius.drNext == drDest )
		{
			fprintf( punch.ipDRout, 
				"change in dest rates, atom=%s%3ld old rate %.2e new %.2e\n", 
			  chDestAtom, istage ,DestOld , DestNew );
		}

		else if( radius.drNext == drOutFrac )
		{
			fprintf( punch.ipDRout, 
				"ratio out2 tot continuua%10.2e%10.2e%7ld %4.4s %4.4s\n", 
			  Out2Tot, 
			  rfield.SumOutMax, 
			  rfield.ipSumOutMax, 
			  rfield.chLineLabel[rfield.ipSumOutMax-1], 
			  rfield.chContLabel[rfield.ipSumOutMax-1] );
			/* in above SumOutMax is anu where peak occured */
		}
#		endif

		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 == 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 == drPressure )
		{
			fprintf( punch.ipDRout, 
				"pressure failure\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/
			  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 nextdr]" );
			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 */
		if( radius.drNext < radius.drMinimum )
		{
			radius.drNext = radius.drMinimum;
			/* 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 */
			conv.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 nextdr to find energy of maximum continuum-gas interaction */
static void ContRate(double *freqm, 
  double *opacm)
{
	long int i, 
	  ipHi,
	  iplow, 
	  limit;
	double FreqH, 
	  FreqSub, 
	  OpacH, 
	  OpacSub, 
	  xMaxH, 
	  xMaxSub;

#	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 
	 */
	xMaxSub = 0.;
	FreqSub = 0.;
	OpacSub = 0.;

	/* 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 below CI edge */
		if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] - 
		  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > xMaxSub )
		{
			xMaxSub = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
			  (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
			FreqSub = rfield.anu[i];
			OpacSub = 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]) > xMaxSub )
		{
			xMaxSub = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
			  (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
			FreqSub = rfield.anu[i];
			OpacSub = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
		}
	}

	/* variables to check continuum interactions over lyman continuum */
	xMaxH = 0.;
	OpacH = 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]) > xMaxH )
		{
			/* xMaxH = anu(i)*flux(i)/widflx(i)*opac(i) */
			xMaxH = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
			  (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
			FreqH = rfield.anu[i];
			OpacH = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
		}
	}


	/* use lyman continuum if its opacity is larger than non-h ion */
	if( xMaxSub < 1e-30 )
	{
		/* this happens for laser source - use Lyman continuum */
		*opacm = OpacH;
		*freqm = FreqH;
	}
	else if( OpacH > OpacSub && xMaxH/xMaxSub > 1e-10 )
	{
		/* use Lyman continuum */
		*opacm = OpacH;
		*freqm = FreqH;
	}
	else
	{
		/* not much rate in the lyman continuum, stick with low energy */
		*opacm = OpacSub;
		*freqm = FreqSub;
	}

	if( xMaxH > xMaxSub )
	{
	}
	else
	{
		*opacm = OpacSub;
		*freqm = FreqSub;
	}

	{
		/* 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", 
			xMaxSub,FreqSub,OpacSub,
			xMaxH,FreqH ,OpacH,*freqm,*opacm
			);
	
		}
	}

	/* these were set to zero at start, and must have been reset if one of the
	 * two loops.  Logic error if still zero. */
	ASSERT( *opacm > 0. && *freqm > 0. );

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

/*GrainRateDr called by nextdr 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 nextdr 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 */
