/* 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 */
/*zero zero out or initialize variables, called by cdInit, but also by optimize_func during optimization,
 * this is called before any commands are parsed, called one time per model, at very start */
/*rfield_optac_zero zero out rfield arrays between certain limits */
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "punch.h"
#include "oxy.h"
#include "doppvel.h"
#include "dense.h"
#include "grains.h"
#include "rt.h"
#include "he.h"
#include "struc.h"
#include "h2.h"
#include "coolheavy.h"
#include "lines.h"
#include "dynamics.h"
#include "carb.h"
#include "sil.h"
#include "mean.h"
#include "helike.h"
#include "atomfeii.h"
#include "iso.h"
#include "conv.h"
#include "geometry.h"
#include "timesc.h"
#include "peimbt.h"
#include "ionbal.h"
#include "continuum.h"
#include "atmdat.h"
#include "mole.h"
#include "ca.h"
#include "input.h"
#include "atoms.h"
#include "pressure.h"
#include "stopcalc.h"
#include "numderiv.h"
#include "wind.h"
#include "negdrg.h"
#include "colden.h"
#include "yield.h"
#include "hmi.h"
#include "rfield.h"
#include "abund.h"
#include "radius.h"
#include "opacity.h"
#include "neutrn.h"
#include "broke.h"
#include "cpu.h"
#include "secondaries.h"
#include "called.h"
#include "phycon.h"
#include "warnings.h"
#include "thermal.h"
#include "cooling.h"
#include "fe.h"
#include "hyperfine.h"
#include "zero.h"

/*zero zero out or initialize variables, called by cdInit, but also by optimize_func during optimization,
 * this is called before any commands are parsed, called one time per model, at very start */
void zero(void)
{
	long int i,
	  ion,
	  ipISO ,
	  nelem,
	  ns;

	/* this is used to signify the first call to this routine.  At that
	 * stage some memory has not been allocated so must not initialize,
	 * set false at very end of this routine */
	static int lgFirstCall = TRUE;

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

	/* this routine is called exactly one time at the start of
	 * the calculation of a single model.  When the code is used as a subroutine
	 * this routine is called one time for each model.  It is called before
	 * the boundary conditions are read in, and is never called again
	 * during that calculation of the one model.  
	 * All default variables that must be initialized before a calculation starts
	 * must appear in the routine.  In a grid they are reset for each model
	 */

	/* these used to be in block data logic */
	Zerologic();

	/* set all initial abundances */
	AbundancesZero();

	/* zero out parameters needed by large FeII atom */
	FeIIZero();

	/* zero out warnings, cautions, notes, etc */
	wcnint();

	/* limits for highest and lowest stages of ionization in TrimStage */
	ionbal.trimhi = 1e-6;
	ionbal.lgTrimhiOn = TRUE;
	ionbal.trimlo = 1e-10;

	hyperfine.lgLya_pump_21cm = TRUE;

	/* variable to do with geometry */
	geometry.nprint = 1000;
	geometry.lgZoneSet = FALSE;
	geometry.lgZoneTrp = FALSE;
	geometry.lgEndDflt = TRUE;

	/* this is default number of zones
	 * >>chng 96 jun 5, from 400 to 500 for thickes corners4 grid */
	/* >>chng 04 jan 30, from 600 to 800, code uses finer zoning today */
	/* >>chng 04 dec 24, from 800 to 1400, so that HII region - molecular cloud
	 * sims do not need set nend - all sims in test suite will run ok without set nend */
	geometry.nEndDflt = 1400;

	for( i=0; i < ITRDIM; i++ )
	{
		geometry.nend[i] = geometry.nEndDflt;
		/*>>chng 03 nov 13, from 1e30 to 1e31, because default inner radius raised to 1e30 */
		radius.router[i] = 1e31;
	}

	/* change amgle of illumination */
	geometry.AngleIllum = 1.;
	geometry.fiscal = 1.;
	geometry.FillFac = 1.;
	geometry.filpow = 0.;

	/* default is open geometry, not sphere */
	geometry.lgSphere = FALSE;
	/* the radiative transport covering factor */
	geometry.covrt = 0.;
	/* the geometric covering factor */
	geometry.covgeo = 1.;
	/* default is expanding when geometry set */
	geometry.lgStatic = FALSE;
	/* option to tell code not to complain when geometry static done without iterating,
	 * set with (OK) option on geometry command */
	geometry.lgStaticNoIt = FALSE;
	/* this is exponent for emissivity contributing to observed luminosity, r^2.
	 * set to 1 with aperture slit, to 0 with aperture beam command */
	geometry.iEmissPower = 2;

	carb.p1909 = 0.;
	carb.p2326 = 0.;
	sil.p1895 = 0.;
	oxy.p1666 = 0.;
	oxy.p1401 = 0.;

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

	/* clear flag indicating the ionization convergence failures 
	 * have ever occured in current zone */
	conv.lgConvIonizThisZone = FALSE;

	/* general abort flag */
	lgAbort = FALSE;

	/* cooling tolerance heating tolerance - allowed error in heating -  cooling balance */
	/*conv.HeatCoolRelErrorAllowed = 0.02f;*/
	/* >>chng 04 sep 25, change te toler from 0.02 to 4x smaller, 0.005, drove instabilities
	 * in chemistry */
	conv.HeatCoolRelErrorAllowed = 0.005f;

	/* this is the default allowed relative error in the electron density */
	conv.EdenErrorAllowed = 0.01;

	conv.LimFail = 20;
	conv.lgMap = FALSE;

	/* true if level 2 lines were contributors to the ots rates, set in dimacool */
	conv.lgLevel2_OTS_Imp = TRUE;
	/* true if level 2 lines were contributors to the cooling, set in dimacool */
	conv.lgLevel2_Cool_Imp = TRUE;

	/* this counts how many times ionize is called in this model after startr,
	 * and is flag used by ionize to understand it is being called the first time*/
	conv.nTotalIoniz = 0;
	/* these are variables to remember the biggest error in the
	 * electron density, and the zone in qhich it occurred */
	conv.BigEdenError = 0.;
	conv.AverEdenError = 0.;
	conv.BigHeatCoolError = 0.;
	conv.AverHeatCoolError = 0.;
	conv.BigPressError = 0.;
	conv.AverPressError = 0.;
	strcpy( conv.chSolverEden , "simple" );
	/* >>chng 03 mar 24, turning this on resulted in host of errors,
	 * nearly all due to more iterations per zone, test suite took longer time
	 * see edensolver.txt in tsuite directory for full list.  Of these, only
	 * orion_hii_pdr, hizqso, and globule, actually had eden faults 
	 * make inital steps to soln larger, not dependent on tolerance */
	/* this did not work - conserv.in failed 
	 * this should work - used for dynamics 
	strcpy( conv.chSolverEden , "new" );*/

	strcpy( conv.chSolverTemp , "simple" );
	/* most models failed with impossible bracket to solver with this on,
	 * check blister.in for one, this does not really work
	strcpy( conv.chSolverTemp , "new" );*/

	/* iterate to convergence flag */
	conv.lgAutoIt = FALSE;
	/* convergence criteria */
	conv.autocv = 0.20f;
	conv.lgConvTemp = TRUE;
	conv.lgConvPres = TRUE;
	conv.lgConvEden = TRUE;
	/* >>chng 04 jan 25, only set lgConvIoniz TRUE where used in ConvXXX path */
	/*conv.lgConvIoniz = TRUE;*/

	/* this option, use the new atmdat_rad_rec recombination rates */
	PhFitOn.lgPhFit = FALSE;

	/* age of the cloud, to check for time-steady */
	timesc.CloudAgeSet = -1.f;
	/* some timescale for CO and H2 */
	timesc.time_H2_Dest_longest = 0.;
	timesc.time_H2_Form_longest = 0.;
	/* remains neg if not evaluated */
	timesc.time_H2_Dest_here = -1.;
	timesc.time_H2_Form_here = 0.;

	timesc.AgeCOMoleDest = 0.;
	timesc.BigCOMoleForm = 0.;

	timesc.TimeH21cm = 0.;
	timesc.sound_speed_isothermal = 0.;

	peimbt.tsqden = 1e7;

	/* CO related variables */
	co.codfrc = 0.;
	co.codtot = 0.;
	co.CODissHeat = 0.;
	/* the CO molecule */
	co.COCoolBigFrac = 0.;
	co.lgCOCoolCaped = FALSE;
	/* largest fraction of atoms in molecules */
	for( i=0; i<NUM_COMOLE_CALC; ++i )
		co.xMoleFracMax[i] = 0.;

	/* flag to turn off molecular network */
	co.lgNoCOMole = FALSE;

	NumDeriv.lgNumDeriv = FALSE;
	/* nsum is line pointer within large stack of line intensities */
	LineSave.nsum = 0;

	/* index within the line in the line stack 
	 * default is Hbeta total - the third line in the stack
	 * 0th is a zero for sanity, 1st is unit, 2nd is a comment */
	/* >>chng 02 apr 22 from 2 to 3 since added unit at 1 */
	LineSave.ipNormWavL = 3;
	LineSave.ipEmerNormWavL = 0;
	LineSave.WavLNorm = 4861.f;
	LineSave.lgNormFound = FALSE;
	LineSave.lgNormSet = FALSE;
	LineSave.sig_figs = 4;

	/* the label for the normalization line */
	strcpy( LineSave.chNormLab, "    " );

	/* the scale factor for the normalization line */
	LineSave.ScaleNormLine = 1.;

	/* this is scale factor, reset with set resoluton command, for setting
	 * the continuum resolution.  Setting to 0.1 will increase resolution by 10x.
	 * this multiplies the resolution contained in the continuum_mesh.dat file */
	continuum.ResolutionScaleFactor = 1.;

	continuum.lgCoStarInterpolationCaution = FALSE;
	continuum.lgCon0 = FALSE;

	/* threshold for faintest heating cooling to punch with 
	 * punch heating or punch cooling commands,
	 * can be reset with PUNCH WEAKHEATCOOL command */
	punch.WeakHeatCool = 0.05f;

	/* set of variables used to control punch command */
	punch.npunch = 0;
	punch.cp_npun = 0;
	punch.lgPunContinuum = FALSE;
	punch.lgHashEndIter = TRUE;
	punch.lgDRPLst = FALSE;
	/* this is the string that will appear after one model in the punch output,
	 * reset with the "set punch hash" command */
	strcpy( punch.chHashString , " ###########################"  );
	/* this is the string that will be added as a prefix to filenames for punch output,
	 * reset with the "set punch prefix" command */
	strcpy( punch.chFilenamePrefix , ""  );
	/* punch every one continuum point */
	punch.ncPunchSkip = 1;
	/* flush file after every iteration - set punch flush */
	punch.lgFLUSH = FALSE;

	/* default line width for plots, in units of speed of light
	 * PunchLWidth = 1.
	 * >>chng 97 jul 10, from unity to 1000 km/sec */
	punch.PunchLWidth = (float)(SPEEDLIGHT/1000.e5);

	/* free free heating, cooling, net */
	CoolHeavy.lgFreeOn = TRUE;
	CoolHeavy.brems_cool_h = 0.;
	CoolHeavy.colmet = 0.;

	CoolHeavy.brems_cool_net = 0.;
	hydro.cintot = 0.;

	/* option to print emissivity instead of intensity/luminosity */
	hydro.lgHydEmiss = FALSE;
	hydro.lgHiPop2 = FALSE;
	hydro.pop2mx = 0.;
	hydro.lgHLyaMased = FALSE;
	hydro.lgReevalRecom = FALSE;

	/* flag for Lya masing */
	hydro.xLaMase = 0.;
	hydro.HCollIonMax = 0.;

	/* this is flag indicating which type of model atom to use, 
	 * set with atom hydrogen populations, departure, lowt */
	/* >>chng 02 mar 13, this had been set to POPU which had the effect
	 * of forcing the populations solution even for very extreme low ionization.
	 * a few models did have neg pops as a result.  There was no reason given for the
	 * change nor was a time stamp present.  So, the code had been running without
	 * this protection for some time.  Changed back to none here, to enable the
	 * test in hydrolevel, and also made limits there much more stringent (since
	 * they had (unintentionally?) been very stringent for some time now). */
	/*strcpy( iso.chTypeAtomSet , "POPU" );*/
	for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		strcpy( iso.chTypeAtomSet[ipISO] , "none" );
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
		{
			/* what was actually used */
			strcpy( iso.chTypeAtomUsed[ipISO][nelem] , "none" );
		}
	}

	/* type of hydrogen atom topoff, options are " add" and "scal" 
	 * in versions 90 this was " add", but was "scal" in 91
	 * >>chng 99 jan 16, changed back to " add"*/
	/*strcpy( hydro.chHTopType, "scal" );*/
	strcpy( hydro.chHTopType, " add" );

	/* Lya excitation temperature, counter for hotter than gas */
	hydro.TexcLya = 0.;
	hydro.TLyaMax = 0.;
	hydro.nLyaHot = 0;

	/* option to kill damping wings of Lya */
	hydro.DampOnFac = 1.;

	/* is continuum pumping of H lyman lines included?  yes, but turned off
	 * with atom h-like lyman pumping off command */
	hydro.lgLymanPumping = TRUE;

	/* mutilplicative factor for all continuum pumping of H I Lyman lines,
	 * account for possible emission in the line */
	hydro.xLymanPumpingScaleFactor = 1.f;

	/* >>refer	abundance	D/H	Pettini, M., & Bowen, D.V., 2001, ApJ, 560, 41 */
	/* quoted error is +/- 0.35 */
	hydro.D2H_ratio = 1.65e-5;

	/* zero fracs of He0 dest due to 23S */
	he.nzone = 0;
	he.frac_he0dest_23S = 0.;
	he.frac_he0dest_23S_photo = 0.;

	 /* flag set by compile he-like command, says to regenerate table of recombination coef */
	helike.lgCompileRecomb = FALSE;
	helike.lgNoRecombInterp = FALSE;
	/* flag saying whether to use a set of physics meant to 
	 * replicate the Benjamin et al. HeI results,
	 * set with ATOM HE-LIKE BENJAMIN command */
	helike.lgSetBenjamin = FALSE;

	/* how the gbar cs will be treated - set with atom he-like gbar command */
	/*TODO	2	change this to CS_new */
	helike.lgCS_Vriens = TRUE;
	helike.lgCS_Vrinceanu = TRUE;
	helike.lgCS_therm_ave = FALSE;
	helike.lgCS_None = FALSE;
	/* when set try actually set to 1 or 2, depending on which fit is to be used,
	 * 1 is the broken power law fit */
	/* >>chng 02 dec 21, change to broken power law fit */
	helike.lgCS_new = 1;
	/* This flag says whether the density is high enough that helium is sufficiently l-mixed. */
	helike.lgCritDensLMix = TRUE;
	/* flag saying whether to include fine-structure mixing in spontaneous decays	
	 * set with ATOM HE-LIKE FSM command */
	helike.lgFSM = 0;
	/* This is the flag saying whether to generate errors for helike rates.  FALSE
	 * means don't generate errors.	*/
	helike.lgRandErrGen = FALSE;
	/* this is the flag saying whether we have done the random error generation yet.  FALSE
	 * means not done yet.	*/
	helike.lgErrGenDone = FALSE;
	/* this is the flag saying whether we should include excess recombination in the
	 * helike sequence.  Should only be off if testing effect of topoff approximations. */
	helike.lgTopoff = TRUE;
	helike.lgHugeCaseB = FALSE;
	/* This sets the default dielectronic recombination for helike ions to the 
	 * Gu 2003 treatment.	*/
	helike.lgDielRecom = 1;
	
	/* the h2 molecule */
	H2_Zero();

	/**********************************************************************
	 * all parameters having to do with secondary ionization 
	 * by suprathermal electrons 
	 **********************************************************************/
	secondaries.SetCsupra = 0.;
	secondaries.lgCSetOn = FALSE;
	secondaries.lgSecOFF = FALSE;
	secondaries.SecHIonMax = 0.;

	secondaries.heatef = 1.;
	secondaries.efionz = 0.;
	secondaries.exctef = 0.;
	secondaries.x12tot = 0.;
	secondaries.sec2total = 0.;

	for( nelem=0; nelem< LIMELM; ++nelem )
	{
		secondaries.Hx12[0][nelem] = 0.;
		secondaries.Hx12[1][nelem] = 0.;
	}
	secondaries.Hx12[0][LIMELM] = 0.;
	secondaries.Hx12[1][LIMELM] = 0.;

	if( lgFirstCall )
	{
		/* malloc space for supra[nelem][ion] */
		secondaries.csupra = (float **)MALLOC( (unsigned)LIMELM*sizeof(float *) );
		secondaries.csupra_effic = (float **)MALLOC( (unsigned)LIMELM*sizeof(float *) );
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
		{
			secondaries.csupra[nelem] = (float *)MALLOC( (unsigned)(nelem+1)*sizeof(float) );
			secondaries.csupra_effic[nelem] = (float *)MALLOC( (unsigned)(nelem+1)*sizeof(float) );
		}
	}
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		for( ion=0; ion<nelem+1; ++ion )
		{
			/* secondary ionization rate for each species */
			secondaries.csupra[nelem][ion] = 0.;
			/* the rate of each species relative to H0 */
			secondaries.csupra_effic[nelem][ion] = 1.f;
		}
	}
	/* this scale factor is from table 10 of Tielens & Hollenbach 1985 */
	secondaries.csupra_effic[ipHELIUM][0] = 1.08f;

	/* on first call, these arrays do not exist, only zero here on 
	 * second and later calls, on first call, create them */
	if( lgFirstCall )
	{
		/* these will save bound electron recoil information data */
		if( (ionbal.ipCompRecoil = 
			(long**)MALLOC(sizeof(long*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (ionbal.CompRecoilIonRate = 
			(double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (ionbal.CompRecoilIonRateSave = 
			(double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (ionbal.CompRecoilHeatRate = 
			(double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (ionbal.CompRecoilHeatRateSave = 
			(double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (ionbal.PhotoRate_Shell = 
			(double****)MALLOC(sizeof(double***)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (ionbal.CollIonRate_Ground = 
			(double***)MALLOC(sizeof(double**)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (ionbal.xInnerShellIonize = 
			(double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();

		/* these are source and sink terms for heavy element ionization balance from the
		 * chemistry */
		if( (mole.source = 
			(double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (mole.sink = 
			(double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ))==NULL )
			BadMalloc();
		if( (mole.xMoleChTrRate = 
			(float***)MALLOC(sizeof(float**)*(unsigned)LIMELM ))==NULL )
			BadMalloc();

		/* space for ionization recombination arrays */
		if( (ionbal.RateIonizTot = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		if( (ionbal.RateRecomTot = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		if( (ionbal.RR_rate_coef_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		if( (ionbal.DR_rate_coef_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		if( (ionbal.RR_Verner_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();

		/* rate coefficients [cm3 s-1] for Badnell DR recombination */
		if( (ionbal.DR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		if( (ionbal.DR_Badnell_rate_coef_mean_ion = (double *)MALLOC(sizeof(double)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		/* do these rate coefficients exist? */
		if( (ionbal.lgDR_Badnell_rate_coef_exist = (int **)MALLOC(sizeof(int *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		if( (ionbal.RR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		/* do these rate coefficients exist? */
		if( (ionbal.lgRR_Badnell_rate_coef_exist = (int **)MALLOC(sizeof(int *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
		/* rate coefficients [cm3 s-1] for older DR recombination */
		if( (ionbal.DR_old_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();

		/* create arrays for ions */
		for(nelem=0; nelem<LIMELM; ++nelem )
		{
			if( (ionbal.DR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.lgDR_Badnell_rate_coef_exist[nelem] = (int *)MALLOC(sizeof(int)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.RR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.lgRR_Badnell_rate_coef_exist[nelem] = (int *)MALLOC(sizeof(int)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.DR_old_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();

			if( (ionbal.RateIonizTot[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.RateRecomTot[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.RR_rate_coef_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.DR_rate_coef_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.RR_Verner_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ) )==NULL )
				BadMalloc();
			if( (ionbal.xInnerShellIonize[nelem] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			if( (ionbal.ipCompRecoil[nelem] = 
				(long*)MALLOC(sizeof(long)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			if( (ionbal.CompRecoilIonRate[nelem] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			if( (ionbal.CompRecoilIonRateSave[nelem] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			if( (ionbal.CompRecoilHeatRate[nelem] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			if( (ionbal.CompRecoilHeatRateSave[nelem] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			if( (ionbal.PhotoRate_Shell[nelem] = 
				(double***)MALLOC(sizeof(double**)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			if( (ionbal.CollIonRate_Ground[nelem] = 
				(double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) ))==NULL )
				BadMalloc();
			/* chemistry source and sink terms for ionization ladders */
			if( (mole.source[nelem] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) ))==NULL )
				BadMalloc();
			if( (mole.sink[nelem] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) ))==NULL )
				BadMalloc();
			if( (mole.xMoleChTrRate[nelem] = 
				(float**)MALLOC(sizeof(float*)*(unsigned)(nelem+2) ))==NULL )
				BadMalloc();

			for( ion=0; ion<nelem+1; ++ion )
			{
				/* >>chng 03 aug 09, set these to impossible values */
				ionbal.RateIonizTot[nelem][ion] = -1.;
				ionbal.RateRecomTot[nelem][ion] = -1.;
				ionbal.xInnerShellIonize[nelem][ion] = -1.;
				ionbal.ipCompRecoil[nelem][ion] = -99;
				ionbal.CompRecoilIonRate[nelem][ion] = -1.;
				ionbal.CompRecoilIonRateSave[nelem][ion] = -1.;
				ionbal.CompRecoilHeatRate[nelem][ion] = -1.;
				ionbal.CompRecoilHeatRateSave[nelem][ion] = -1.;

				/* finish mallocing space */
				if( (ionbal.PhotoRate_Shell[nelem][ion] = 
					(double**)MALLOC(sizeof(double*)*(unsigned)NSHELLS ))==NULL )
					BadMalloc();
				if( (ionbal.CollIonRate_Ground[nelem][ion] = 
					(double*)MALLOC(sizeof(double)*(unsigned)2 ))==NULL )
					BadMalloc();
				for( ns=0; ns<NSHELLS; ++ns )
				{
					if( (ionbal.PhotoRate_Shell[nelem][ion][ns] = 
						(double*)MALLOC(sizeof(double)*(unsigned)3 ))==NULL )
						BadMalloc();
				}

				/* now set to impossible values */
				ionbal.ipCompRecoil[nelem][ion] = -100000;
				ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.;
				ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.;
				ionbal.DR_old_rate_coef[nelem][ion] = 0.;
				SET_DBL_SNAN(ionbal.DR_rate_coef_used[nelem][ion] );
				SET_DBL_SNAN(ionbal.RR_rate_coef_used[nelem][ion] );
				SET_DBL_SNAN(ionbal.RR_Verner_rate_coef[nelem][ion] );
			}
		}
	}
	for(ion=0; ion<LIMELM; ++ion )
	{
		ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 0.;
	}

	/* now zero out these arrays */
	for( nelem=0; nelem< LIMELM; ++nelem )
	{
		for( ion=0; ion<nelem+1; ++ion )
		{

			ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
			ionbal.CompRecoilIonRate[nelem][ion] = 0.;
			ionbal.xInnerShellIonize[nelem][ion] = 0.;
			ionbal.CollIonRate_Ground[nelem][ion][0] = 0.;
			ionbal.CollIonRate_Ground[nelem][ion][1] = 0.;
			ionbal.RateRecomTot[nelem][ion] = 0.;
			for( ns=0; ns < NSHELLS; ++ns )
			{
				/* must be zero since ion routines use these when
				 * not yet defined */
				ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 0.;
				ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.;
				ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.;
			}
		}
		/* these have one more ion than above */
		for( ion=0; ion<nelem+2; ++ion )
		{
			long int ion2;
			/* zero out the source and sink arrays */
			mole.source[nelem][ion] = 0.;
			mole.sink[nelem][ion] = 0.;
			if( (mole.xMoleChTrRate[nelem][ion] = 
				(float*)MALLOC(sizeof(float)*(unsigned)(nelem+2) ))==NULL )
				BadMalloc();
			for( ion2=0; ion2<nelem+2; ++ion2 )
			{
				mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
			}
		}
	}

	/* capture of molecules onto grain surfaces - formation of ices
	 * flag says to include this process - turned off with the
	 * command NO GRAIN MOLECULES */
	mole.lgGrain_mole_deplete = TRUE;

	ionbal.lgPhotoIoniz_On = TRUE;
	ionbal.lgCompRecoil = TRUE;
	ionbal.lgInnerShellLine_on = TRUE;
	ionbal.lgInnerShell_Romas = TRUE;

	for( i=0; i<4; ++i )
	{
		ionbal.GuessDiel[i] = 1.;
	}

	/* default condition is burgess supressed, nusbaumer and storye not */
	ionbal.lgSupDie[0] = TRUE;
	ionbal.lgSupDie[1] = FALSE;

	ionbal.lgNoCota = FALSE;
	for( i=0; i < LIMELM; i++ )
	{
		ionbal.CotaRate[i] = 0.;
	}
	ionbal.ilt = 0;
	ionbal.iltln = 0;
	ionbal.ilthn = 0;
	ionbal.ihthn = 0;
	ionbal.ifail = 0;
	ionbal.lgGrainIonRecom = TRUE;

	/* turn on Badnell DR by default */
	ionbal.lgDR_recom_Badnell_use = TRUE;
	/* do not turn on Badnell DR by default */
	ionbal.lgDR_recom_Badnell_use = FALSE;
	/* turn on Badnell RR by default */
	ionbal.lgRR_recom_Badnell_use = TRUE;
	/* do not turn on Badnell RR by default */
	ionbal.lgRR_recom_Badnell_use = FALSE;
	/* option to print recom coef then exit */
	ionbal.lgRecom_Badnell_print = FALSE;
	/* set true to use Badnell mean in place of existing hacks */
	ionbal.lg_use_DR_Badnell_rate_coef_mean_ion = FALSE;
	ionbal.lg_use_DR_Badnell_rate_coef_mean_ion = TRUE;

	/* setting this non-zero will turn on guess for dr for entire range of ions */
	/* >>chng 03 nov 22, change from false to true - turn on process */
	/* >>refer	dr	guess	Kraemer, S.B., Ferland, G.J., & Gabel,  J.R. 2004, ApJ in press */
	ionbal.lg_guess_coef = TRUE;
	ionbal.guess_noise = 0.;

	/* do O H charge transfer in chemistry by default, changes with
	 * set HO charge transfer */
	ionbal.lgHO_ct_chem = TRUE;

	/**********************************************************************
	 * these are options to print errors to special window,
	 * set with print errors command, 
	 * output will go to standard error
	 * defined in cdInit 
	 **********************************************************************/
	lgPrnErr = FALSE;
	ioPrnErr = stderr;
	lgAssertFPE = FALSE;

	/* the code is ok at startup */
	broke.lgBroke = FALSE;
	broke.lgFixit = FALSE;
	broke.lgCheckit = FALSE;

	/* the ipSolar array is a set of pointers used for reading
	 * in abundances with the "abundances" command
	 * the hydrogen abundance of unity is not read in */
	abund.npSolar = LIMELM - 1;
	for( i=0; i < abund.npSolar; i++ )
	{
		abund.ipSolar[i] = i + 2;
	}

	/* turn element on if this is first call to this routine,
	 * but do not reset otherwise since would clobber how elements were set up */
	{
		static int nCalled = 0;
		if( nCalled == 0 )
		{
			for( nelem=0; nelem < LIMELM; nelem++ )
			{
				/* never turn on element if turned off on first iteration */
				dense.lgElmtOn[nelem] = TRUE;

				/* option to set ionization with element ioniz cmnd,
				* default (FALSE) is to solve for ionization */
				dense.lgSetIoniz[nelem] = FALSE;
			}
		}
		++nCalled;
	}

	/* main arrays to save ionization fractions*/
	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		dense.gas_phase[nelem] = 0.;
		dense.xMolecules[nelem] = 0.;
		for( ion=0; ion < LIMELM+1; ion++ )
		{
			dense.xIonDense[nelem][ion] = 0.;
		}
	}
	dense.xMassTotal = 0.;
	/* list of atomic weights, mass in amu, used for thermal line widths */
	/* >>refer	all	atomic weight	Coplen, T.B. 2001, J. Phys. Chem REf Data, 30, 701 */
	dense.AtomicWeight[0] = 1.00794f;
	dense.AtomicWeight[1] = 4.0026022f;
	dense.AtomicWeight[2] = 6.9412f;
	dense.AtomicWeight[3] = 9.0121823f;
	dense.AtomicWeight[4] = 10.8117f;
	dense.AtomicWeight[5] = 12.01078f;
	dense.AtomicWeight[6] = 14.00672f;
	dense.AtomicWeight[7] = 15.99943f;
	dense.AtomicWeight[8] = 18.9984032f;
	dense.AtomicWeight[9] = 20.17976f;
	dense.AtomicWeight[10] = 22.989770f;
	dense.AtomicWeight[11] = 24.30506f;
	dense.AtomicWeight[12] = 26.9815382f;
	dense.AtomicWeight[13] = 28.08553f;
	dense.AtomicWeight[14] = 30.9737612f;
	dense.AtomicWeight[15] = 32.0655f;
	dense.AtomicWeight[16] = 35.4532f;
	dense.AtomicWeight[17] = 39.9481f;
	dense.AtomicWeight[18] = 39.09831f;
	dense.AtomicWeight[19] = 40.0784f;
	dense.AtomicWeight[20] = 44.9559108f;
	dense.AtomicWeight[21] = 47.8671f;
	dense.AtomicWeight[22] = 50.94151f;
	dense.AtomicWeight[23] = 51.99616f;
	dense.AtomicWeight[24] = 54.9380499f;
	dense.AtomicWeight[25] = 55.8472f;
	dense.AtomicWeight[26] = 58.9332009f;
	dense.AtomicWeight[27] = 58.69342f;
	dense.AtomicWeight[28] = 63.5463f;
	dense.AtomicWeight[29] = 65.392f;

	/* this is the simple fred hamann feii atom */
	for( i=0; i < NFEII; i++ )
	{
		fe2ovr_la.Fe2PopLte[i] = 0.;
		fe2ovr_la.feopc[i] = 0.;
		fe2ovr_la.Fe2TauLte[i] = opac.taumin;
	}

	/* zero out volume and column density save arrays */
	MeanZero();

	/* parameters to do with wind */
	wind.lgWindOK = TRUE;
	wind.windv0 = 0.;
	wind.comass = 0.;
	wind.windv = 0.;
	wind.emdot = 0.;
	wind.AccelAver = 0.;
	wind.acldr = 0.;
	wind.agrav = 0.;
	wind.AccelTot = 0.;
	wind.AccelCont = 0.;
	wind.AccelLine = 0.;
	wind.AccelPres = 0.;
	wind.AccelMax = 0.;
	wind.fmul = 0.;
	wind.lgVelPos = TRUE;

	/* zero out heating rates */
	HeatZero();

	/* some parameters dealing with calcium */
	ca.Ca2RmLya = 0.;
	ca.popca2ex = 0.;
	ca.oldcla = 0.;
	ca.Ca3d = 0.;
	ca.Ca4p = 0.;
	ca.dstCala = 0.;

	/* zero out molecular species */
	for( i=0; i < NUM_COMOLE_CALC; i++ )
	{
		co.hevmol[i] = 0.;
		co.hevcol[i] = 0.;
	}
	/* ratio of C12O16 to C13O16 */
	co.RatioC12O16_2_C13O16 = 30.;

	/* flag saying that H2O water dest rate went to zero */
	co.lgH2Ozer = FALSE;

	/* extra electron density, set with eden command */
	dense.EdenExtra = 0.;

	/* forced electron density, set with set eden command */
	dense.EdenSet = 0.;

	/* this is the default allowed relative error in the pressure */
	conv.PressureErrorAllowed = 0.01f;

	/* this is abort option set with SET PRESIONIZ command */
	conv.limPres2Ioniz = 100000000;

	/* some titles and line images */
	for( i=0; i<74; ++i)
	{
		input.chTitle[i] = ' ';
	}
	input.chTitle[75] = '\0';

	/* velocity field information */
	/* the turbulent velocity at illuminated face, internall in cm/s,
	 * but entered with turbulence command in km/s */
	DoppVel.TurbVel = 0.;
	/* The scale in cm over which the turbulence is dissipated.  Normally 0,
	 * only set if dissipate keyword appears on turbulence command */
	DoppVel.DispScale = 0.;
	/* equipartition option on turbulence command, to set turbulence from B */
	DoppVel.lgTurbEquiMag = FALSE;

	/* flag for constant pressure, include continuum too */
	pressure.lgContRadPresOn = TRUE;
	/* pressure related variables */
	pressure.PresInteg = 0.;
	pressure.pinzon = 0.;

	/* individual terms for the pressure equation of state */
	/* >>chng 05 jan 01, all three are set true at start *
	 * code default is constant density, so all this is ignored
	 * is turned to constant pressure then these flags must be adjusted so
	 * that we get the pressure we expect.  with all true, all three
	 * contributors to pressure will be counted - with constant gas pressure
	 * these are turned off */
	pressure.lgPres_radiation_ON = TRUE;
	pressure.lgPres_magnetic_ON = TRUE;
	pressure.lgPres_ram_ON = TRUE;

	pressure.PresRamCurr = 0.;
	pressure.PresRadCurr = 0.;
	pressure.lgPradCap = FALSE;
	pressure.lgPradDen = FALSE;
	pressure.lgLineRadPresOn = TRUE;
	/* normally abort when radiation pressure exceeds gas pressure in const pres mod,
	 * this is option to keep going, set with NO ABORT on constant pressure command */
	pressure.lgRadPresAbortOK = TRUE;
	/* Ditto for whether to stop at sonic point, this gets set to false
	 * for some of the dynamics pressure modes (strongd, shock, antishock)*/
	pressure.lgSonicPointAbortOK = TRUE;
	/* constant density is the default */
	strcpy( dense.chDenseLaw, "CDEN" );
	/* this flag will say we hit the sonic point */
	pressure.lgSonicPoint = FALSE;
	/* True when no physical solution for desired pressure in strong D fronts */
	pressure.lgStrongDLimbo = FALSE; 


	pressure.RadBetaMax = 0.;
	pressure.ipPradMax_nzone = -1;
	pressure.PresMax = 0.;

	/* zero out some dynamics variables */
	DynaZero();

	/* check whether this is a 32-bit machine
	 * use te which is a float.
	 * then set te after test */
	{
		float flt_test;
		flt_test = 1e-37f;
		flt_test = (float)(flt_test /1e40);
		if( flt_test == 0. )
		{
			lgBit32 = TRUE;
		}
		else
		{
			lgBit32 = FALSE;
		}
	}

	/* now set physical conditions array */
	phycon.te = 1e4;
	phycon.sqrte = 100.;
	phycon.lgPhysOK = TRUE;
	/* largest relative changes in Te, ne, H+, H2, and CO in structure
	 * this is computed as part of prtcomment so does not exist when code not talking,
	 * set to zero in zero and still zero if prtcomment not called */
	phycon.BigJumpTe = 0.;
	phycon.BigJumpne = 0.;
	phycon.BigJumpH2 = 0.;
	phycon.BigJumpCO = 0.;

	dense.xNucleiTotal = 1.;
	/* WJH */
	dense.xMassDensity0 = -1.0f; 
	
	dense.eden = 1.;
	/* >>chng 02 feb 28, this was set to 1 here, and to -99 in zerologic, so that
	 * check whether hden was specified was defeated.  Reset to correct -99. here.*/
	dense.gas_phase[ipHYDROGEN] = -99.;

	/* controls density fluctuations, when true variations are due to density changing,
	 * when false are due to abundances changing if changes are on */
	dense.lgDenFlucOn = TRUE;
	dense.lgDenFlucRadius = TRUE;
	dense.flong = 0.;
	dense.flcPhase = 0.;
	dense.HCorrFac = 1.f;

	/* following will force updaing all temperature - density variables */
	tfidle(TRUE);

	atoms.nNegOI = 0;
	for( i=0; i< N_OI_LEVELS; ++i )
	{
		atoms.popoi[i] = 0.;
	}

	/* do we want to punch negative opacities */
	opac.lgNegOpacIO = FALSE;

	/* this is flag for turning on case b */
	opac.lgCaseB = FALSE;

	/* this is separate flag for turning off collisions from n=2 */
	opac.lgCaseB_HummerStorey = FALSE;

	/* this is separate flag for turning off excited state photoionization */
	opac.lgCaseB_no_photo = FALSE;
	/* another case b option, turn off background opacities, no Pdest */
	opac.lgCaseB_no_pdest = FALSE;

	strcpy( opac.chOpcTyp, "    " );

	/* smallest allowed line and Lya optical depths, reset with
	 * caseb command */
	opac.taumin = 1e-20f;
	opac.tlamin = 1e-20f;

	opac.otsmin = 0.;

	/* this flag says to use the static opacities,
	 * only evaluate them at start of each new zone.
	 * when set false with 
	 * no static opacities
	 * command, always reevalute them */
	opac.lgOpacStatic =  TRUE;

	/* set true in radinc if negative opacites ever occur */
	opac.lgOpacNeg = FALSE;

	/* can turn of scattering opacities for some geometries */
	opac.lgScatON = TRUE;

	/* variables having to do with compiling and/or using the
	 * ancillary file of stored opacities */
	opac.lgCompileOpac = FALSE;
	/* "no file opacity" command sets following var false, says not to use file */
	/*TODO	2	file opacities are disabled for now - reinstate this when arrays settle down */
	opac.lgUseFileOpac = FALSE;

	opac.telec = opac.taumin;
	opac.thmin = opac.taumin;

	opac.tpcah[0] = opac.taumin;
	opac.tpcah[1] = 1e20f;

	opac.lgTauOutOn = FALSE;

	dynamics.dDensityDT = 0.;
	negdrg.lgNegGrnDrg = FALSE;

	/* various criteria for stopping model */
	/* >>chng 04 dec 21, remove from here and init to 1e30 in zero */
	/*StopCalc.tauend = 0.;*/
	StopCalc.tauend = 1e30f;

	/* >>chng 05 nov 22 - NPA.  Stop calculation when fraction of oxygen frozen
	 * out on grains ices gets too high - formation of ices */
	/*StopCalc.StopDepleteFrac = 0.99f;*/
	/* >>chng 05 dec 16, with revised ion solver logic, code should be able to
	 * converge away from situation where ices have disturbed the chemistry and
	 * net negative atomic abundances result.  now we say soln not converged and
	 * soldier on 
	 * this test should not be necessary */
	StopCalc.StopDepleteFrac = 1.02f;

	StopCalc.xMass = 0.;
	StopCalc.taunu = 0.;
	StopCalc.iptnu = -1;
	/* stopping AV */
	StopCalc.AV_extended = 1e30f;
	StopCalc.AV_point = 1e30f;
	/* highest allowed temperature */
	StopCalc.T2High = 1e12f;
	/* highest and lowest temperatues to ever allow */
	/*TODO	3	next two should be constants since not reset - are hardwired limits to 
	 * code's range of temperature validity */
	StopCalc.TeHighest = 1e10f;
	StopCalc.TeLowest = 2.8f;
	/* the floor sets a limit to the temperature in the calculation -
	 * if te falls below this, we do a constant temperature cloud at
	 * this temperature */
	StopCalc.TeFloor = 0.;

	/* this is in cddefines.h and is 4000 */
	StopCalc.tend = TEND;
	/* ending column densities */
	StopCalc.HColStop = COLUMN_INIT;
	StopCalc.colpls = COLUMN_INIT;
	StopCalc.colnut = COLUMN_INIT;
	StopCalc.col_h2 = COLUMN_INIT;
	StopCalc.col_h2_nut = COLUMN_INIT;
	StopCalc.col_H0_ov_Tspin = COLUMN_INIT;
	StopCalc.col_monoxco = COLUMN_INIT;

	/* stopping electron density */
	StopCalc.StopElecDensity = -COLUMN_INIT;

	/* stopping electron and molecular fractions */
	StopCalc.StopElecFrac = -FLT_MAX;
	StopCalc.StopHPlusFrac = -FLT_MAX;
	/* stopping molecular fraction has opposite sign - want to stop when 2H_2/NH gt this */
	StopCalc.StopH2MoleFrac = FLT_MAX;
	/* this flag says that 21cm line optical depth is the stop quantity */
	StopCalc.lgStop21cm = FALSE;

	/* limits to stop line command */
	StopCalc.nstpl = 0;
	for( i=0; i < MXSTPL; i++ )
	{
		StopCalc.stpint[i] = FLT_MAX;
		StopCalc.ipStopLin1[i] = 0;
		StopCalc.ipStopLin2[i] = 0;
		StopCalc.StopLineWl[i] = 0;
		StopCalc.StopLineWl2[i] = 0;
		strcpy(StopCalc.chStopLabel[i] , "none" );
		strcpy(StopCalc.chStopLabel2[i] , "none" );
	}

	/* effects of fast neutrons */
	neutrn.frcneu = 0.;
	neutrn.effneu = 1.;
	neutrn.totneu = 0.;
	neutrn.lgNeutrnHeatOn = FALSE;
	neutrn.CrsSecNeutron = 4e-26;

	opac.stimax[0] = 0.;
	opac.stimax[1] = 0.;
	
	/* this is limit on ionization we shall check for zoning
	 * and prtcomment */
	struc.dr_ionfrac_limit = 1e-3f;
	struc.nzone = -1;
	for (i=0;i<N_H_MOLEC;i++) 
	{
		/* number of protons in each molecule */
		int hmi_protons[N_H_MOLEC] = {1,1,1,2,2,3,2,1};
		/* >>chng 03 nov 28, add electrons to hmi struc */
		/* number of electrons in each molecule, THAT IS NOT COUNTED AS ION,
		 * used to get H mole electron sum */
		int hmi_electrons[N_H_MOLEC] = {0,0,-1,0,1,1,0,1};

		hmi.Hmolec[i] = 0.;
		hmi.nProton[i] = hmi_protons[i];
		/* number of free electrons contributed, -1 for H- */
		hmi.nElectron[i] = hmi_electrons[i];
	}

	hmi.H2_frac_abund_set = 0.;
	hmi.hmihet = 0.;
	hmi.h2plus_heat = 0.;
	hmi.HeatH2Dish_used = 0.;
	hmi.HeatH2Dexc_used = 0.;
	hmi.HeatH2Dish_TH85 = 0.;
	hmi.HeatH2Dexc_TH85 = 0.;
	hmi.HeatH2Dish_BigH2 = 0.;
	hmi.HeatH2Dexc_BigH2 = 0.;
	hmi.UV_Cont_rel2_Draine_DB96_face = 0.;
	hmi.UV_Cont_rel2_Draine_DB96_depth = 0.;
	hmi.UV_Cont_rel2_Habing_TH85_face = 0.;
	hmi.UV_Cont_rel2_Habing_TH85_depth = 0.;
	hmi.HeatH2DexcMax = 0.;
	hmi.CoolH2DexcMax = 0.;
	hmi.hmitot = 0.;
	hmi.H2Opacity = 0.;
	hmi.hmidep = 1.;
	hmi.h2dep = 1.;
	hmi.h2pdep = 1.;
	hmi.h3pdep = 1.;
	hmi.BiggestH2 = -1.f;
	/* option to turn off H molecules */
	hmi.lgNoH2Mole = FALSE;

	/* option to kill effects of H2 in CO chemistry - set with
	 * set leiden hack h2* off */
	hmi.lgLeiden_Keep_ipMH2s = TRUE;
	hmi.lgLeidenCRHack = TRUE;

	/* option to turn on the UMIST rates, naturally this will be 1, set to zero
	   with the set UMIST rates command */
	co.lgUMISTrates = TRUE;

	/* option to use diffuse cloud chemical rates from Table 8 of
	 * >> refer Federman, S. R. & Zsargo, J. 2003, ApJ, 589, 319
	 * By default, this is FALSE - changed with set chemistry command */
	co.lgFederman = TRUE;

	/* this says which guess of the rate of the solomon process to use,
	 * default is Tielens & Hollenbach 1985a, changed with 
	 * set h2 solomon command, options are TH85 and BD96,
	 * the second for the Bertoldi & Draine rates - they
	 * differ by 1 dex.   when large H2 turned on this is ignored */
	/* the Tielens & Hollenbach 1985 treatment */
	hmi.chSolomon = 'T';
	/* the Bertoldi & Draine 1996 treatment */
	/* >>chng 03 nov 15, change default to BD96 */
	hmi.chSolomon = 'B';
	/* >>chng 05 dec 08, use the Elwert et al. approximations as the default */
	hmi.chSolomon = 'E';

	/* default grain formation pumping - Takahashi 2001 */
	hmi.chGrainFormPump = 'T';

	/* set which approximation for Jura rate - Cazaux & Tielens 
	 * >>refer	H2	form	Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */
	hmi.chJura = 'C';

	/* scale factor to multiply Jura rate, set jura rate command */
	hmi.ScaleJura = 1.f;

	/* binding energy for change in H2 population while on grain surface,
	 * set with "set h2 Tad" command */
	hmi.Tad = 800.;

	/* zero out some column densities */
	for( i=0; i < NCOLD; i++ )
	{
		colden.colden[i] = 0.;
	}
	colden.He123S = 0.;
	colden.coldenH2_ov_vel = 0.;
	/* ortho and para column densities */
	h2.ortho_colden = 0.;
	h2.para_colden = 0.;
	for( i=0; i < 5; i++ )
	{
		colden.C2Pops[i] = 0.;
		colden.C2Colden[i] = 0.;
		/* pops and column density for SiII atom */
		colden.Si2Pops[i] = 0.;
		colden.Si2Colden[i] = 0.;
	}
	for( i=0; i < 3; i++ )
	{
		/* pops and column density for CI atom */
		colden.C1Pops[i] = 0.;
		colden.C1Colden[i] = 0.;
		/* pops and column density for OI atom */
		colden.O1Pops[i] = 0.;
		colden.O1Colden[i] = 0.;
		/* pops and column density for CIII atom */
		colden.C3Pops[i] = 0.;
	}
	for( i=0; i < 4; i++ )
	{
		/* pops and column density for CIII atom */
		colden.C3Colden[i] = 0.;
	}

	/* variables to do with Jeans mass and radius */
	colden.TotMassColl = 1.;
	colden.tmas = 1.;
	colden.wmas = 1.;
	colden.rjnmin = FLT_MAX;
	colden.ajmmin = FLT_MAX;

	/* >>chng 01 jul 26, moved next statement from below loop to avoid bug in gcc 2.95.3, PvH */
	/* default diffuse fields is outward only */
	strcpy( rfield.chDffTrns, "OU2" );

	/* >>chng 05 jun 30, malloc rather than declare these arrays
	 * three arrays that are needed to save tabulated continua, LIMSPC is limit to number
	 * of continua that can be entered */
	if( lgFirstCall )
	{
		if( (rfield.tNuRyd = (float**)MALLOC((size_t)(LIMSPC*sizeof(float*)) ) ) == NULL )
			BadMalloc();
		if( (rfield.tslop = (float**)MALLOC((size_t)(LIMSPC*sizeof(float*)) ) ) == NULL )
			BadMalloc();
		if( (rfield.tFluxLog = (float**)MALLOC((size_t)(LIMSPC*sizeof(float*)) ) ) == NULL )
			BadMalloc();
		if( (rfield.lgContMalloc = (int*)MALLOC((size_t)(LIMSPC*sizeof(int)) ) ) == NULL )
			BadMalloc();
	}
	/* this flag will say that the above vector has been malloced into the continuum array 
	 * this can be returned once the continuum is generated */
	for( i=0; i<LIMSPC; ++i )
	{
		rfield.lgContMalloc[i] = FALSE;
	}

	rfield_opac_zero( 0 , rfield.nupper );

	/* strcpy( rfield.chDffTrns, "OU2" ); */
	rfield.lgOutOnly = TRUE;
	rfield.lgUSphON = FALSE;

	/* flags for whether continuum is defined over all energies */
	rfield.lgMMok = TRUE;
	rfield.lgHPhtOK = TRUE;
	rfield.lgXRayOK = TRUE;
	rfield.lgGamrOK = TRUE;

	/* these are the low and high energy bounds of the continuum */
	rfield.emm = 1.001e-8f;
	rfield.egamry = 7.354e6f;

	rfield.nflux = rfield.nupper;

	/* where cloud is optically thick to brems */
	rfield.ipEnergyBremsThin = 0;
	rfield.EnergyBremsThin = 0.;

	/* this is the faintest the high-energy tail of the continuum be */
	rfield.FluxFaint = 1e-10f;

	for( i=0; i < LIMSPC; i++ )
	{
		/* range(i,1) = 1.
		 * >>chng 96 dec 18, from inf mass ryd to H mass ryg */
		rfield.range[i][0] = HIONPOT;
		rfield.range[i][1] = rfield.egamry;
		rfield.ioTableRead[i] = NULL;
	}
	rfield.comtot = 0.;
	rfield.comoff = 1.;
	rfield.cmcool = 0.;
	rfield.cinrat = 0.;
	rfield.extin_mag_B_point = 0.;
	rfield.extin_mag_V_point = 0.;
	rfield.extin_mag_B_extended = 0.;
	rfield.extin_mag_V_extended = 0.;
	rfield.EnerGammaRay = 7676.;

	/* variables dealing with the radius */
	radius.rinner = 0.;
	radius.distance = 0.;
	radius.Radius = 0.;
	radius.Radius_mid_zone = 0.;
	radius.depth = DEPTH_OFFSET;
	radius.depth_mid_zone = DEPTH_OFFSET/2.;
	radius.depth_x_fillfac = 0.;
	radius.lgRadiusKnown = FALSE;
	radius.drad = 0.;
	radius.r1r0sq = 1.;
	/* this is changed with the roberto command, to go from out to in */
	radius.dRadSign = 1.;
	radius.lgDrNeg = FALSE;

	/* RDFALT is log of default starting radius (cm) */
	/* >>chng 03 nov 12, from 25 to 30 for lya clouds */
	/*radius.rdfalt = 25.;*/
	radius.rdfalt = 30.;

	/* set default cylinder thickness */
	radius.CylindHigh = 1e35f;
	radius.lgCylnOn = FALSE;

	radius.drad_x_fillfac = 1.;
	radius.drad_x_fillfac_mean = 1.;
	radius.dVeff = 1.;
	radius.dRNeff = 1.;
	radius.lgdR2Small = FALSE;

	radius.glbdst = 0.;
	radius.glbrad = 0.;

	radius.sdrmin = SMALLFLOAT;
	radius.sdrmax = 1e30;
	radius.lgSMinON = FALSE;
	radius.lgDrMnOn = TRUE;

	radius.lgDrMinUsed = FALSE;

	/* drChange was reset to get orion flux in h-beta correct
	 * drChange is really tau of current zone */
	radius.drChange = 0.15f;

	rfield.lgHabing = FALSE;
	/* this flag to make sure more than one table read does not occur in same run */
	rfield.lgTableRead = FALSE;

	/* flag to turn off Lya ots */
	rfield.lgLyaOTS = TRUE;
	/* heii rec and lya ots */
	rfield.lgHeIIOTS = TRUE;
	rfield.lgKillOTSLine = FALSE;
	rfield.lgKillOutLine = FALSE;
	rfield.lgKillOutCont = FALSE;

	/* rfield.DiffPumpOn is unity unless process disabled by setting to 1
	 * with no diffuse line pumping command */
	rfield.DiffPumpOn = 1.;

	rfield.lgInducProcess = TRUE;
	/* 02 jun 13, by Ryan...added this line	*/
	rfield.lgCompileGauntFF = FALSE;

	/* >>chng 03 nov 28, add option to not do line transfer */
	rfield.lgDoLineTrans = TRUE;

	/* flag saying whether to constantly reevaluated opacities -
	 * set false with no opacity reevaluate command */
	rfield.lgOpacityReevaluate = TRUE;

	/* flag saying whether to constantly reevaluated ionization -
	 * set false with no ionization reevaluate command */
	rfield.lgIonizReevaluate = TRUE;

	/* this flag says that CMB has been set */
	rfield.lgCMB_set = FALSE;

	/* line overlap opacity, turn off with no fine opacity command */
	rfield.lgOpacityFine = TRUE;
	/*rfield.lgOpacityFine = FALSE;fixit(); rm this line */
	/* this element is default for choosing line width */
	rfield.fine_opac_nelem = ipIRON;
	/* there will be this many resolution elements in one fwhm for this element,
	 * at the lowest temperature to be considered */
	rfield.fine_opac_nresolv = 1;

	/* will fine optical depths be punched? */
	rfield.lgPunchOpacityFine = FALSE;

	/* reset some variable related to cooling */
	CoolZero();

	thermal.lgColNeg = TRUE;
	thermal.lgCNegChk = TRUE;
	thermal.CoolHeatMax = 0.;
	thermal.wlCoolHeatMax = 0;
	thermal.totcol = 0.;
	thermal.heatl = 0.;
	thermal.coolheat = 0.;
	thermal.lgCExtraOn = FALSE;
	thermal.CoolExtra = 0.;
	thermal.ctot = 1.;

	thermal.HeatNet = 0.;
	thermal.htot = 1.;
	thermal.power = 0.;
	thermal.FreeFreeTotHeat = 0.;
	thermal.char_tran_cool = 0.;
	thermal.char_tran_heat = 0.;

	fnzone = 0.;
	nzone = 0;
	/* save initial condition for talk in case PRINT LAST used */
	called.lgTalkSave = called.lgTalk;

	thermal.thist = 0.;
	thermal.tlowst = 1e20f;
	thermal.lgTSetOn = FALSE;
	thermal.ConstTemp = 0.;
	thermal.ConstGrainTemp = 0.;
	thermal.lgTeHigh = FALSE;
	thermal.nUnstable = 0;
	thermal.lgUnstable = FALSE;
	thermal.lgTeBD96 = FALSE;
	thermal.lgTLaw = FALSE;
	thermal.lgTeSN99 = FALSE;

	/* try toe predict next zone's temperature in constant density models,
	 * as done in ZoneStart.  Turned off with no tepred command */
	thermal.lgPredNextTe = TRUE;

	oxy.poiii2 = 0.;
	oxy.poiii3 = 0.;
	oxy.poiexc = 0.;

	oxy.d5007r = 0.;
	oxy.d5007t = 0.;
	oxy.d4363 = 0.;
	oxy.d6300 = 0.;

	atmdat.nsbig = 0;

	/* option to turn off collisional ionization with "no collisional ionization" cmmnd */
	atmdat.lgCollIonOn = TRUE;

	/***************************************************
	 * charge transfer ionization and recombination 
	 ***************************************************/
	/* HCharHeatMax, HCharCoolMax are largest fractions of heating in cur zone
	 * or cooling due to ct */
	atmdat.HCharHeatMax = 0.;
	atmdat.HCharCoolMax = 0.;

	atmdat.HIonFrac = 0.;
	atmdat.HCharExcIonTotal = 0.;
	atmdat.HIonFracMax = 0.;
	atmdat.HCharExcRecTotal = 0.;
	/* option to turn off all charge transfer, turned off with no charge transfer command */
	atmdat.lgCTOn = TRUE;

	/* flag saying that charge transfer heating should be included,
	 * turned off with no CTHeat commmand */
	atmdat.HCharHeatOn = 1.;
	for( nelem=0; nelem< LIMELM; ++nelem )
	{
		for( ion=0; ion<LIMELM; ++ion )
		{
			atmdat.HCharExcIonOf[nelem][ion] = 0.;
			atmdat.HCharExcRecTo[nelem][ion] = 0.;
		}
	}

	/* >>chng 97 jan 6, from 0 to 8.5e-10*q as per Alex Dalgarno e-mail
	 * >>chng 97 feb 6, from 8.5e-10*q 1.92e-9x as per Alex Dalgarno e-mail */
	atmdat.HCTAlex = 1.92e-9;

	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		abund.solar[nelem] = abund.SolarSave[nelem];
		/* these are depltion scale factors */
		abund.depset[nelem] = 1.;
		/*begin sanity check */
		if( abund.solar[nelem] == 0. || abund.depset[nelem] == 0. )
		{
			fprintf( ioQQQ, " ZERO finds insane abundance or depletion.\n" );
			fprintf( ioQQQ, " atomic number=%6ld abundance=%10.2e depletion=%10.2e\n", 
			  nelem, abund.solar[nelem], abund.depset[nelem] );
			ShowMe();
			puts( "[Stop in zero]" );
			cdEXIT(EXIT_FAILURE);
		}
		/*end sanity check */
	}

	abund.lgAbTaON = FALSE;
	abund.lgAbnSolar = FALSE;

	/* typical ISM depletion factors, subjective mean of Cowie and Songaila
	 * and Jenkins 
	 * */
	abund.Depletion[0] = 1.;
	abund.Depletion[1] = 1.;
	abund.Depletion[2] = .16f;
	abund.Depletion[3] = .6f;
	abund.Depletion[4] = .13f;
	abund.Depletion[5] = 0.4f;
	abund.Depletion[6] = 1.0f;
	abund.Depletion[7] = 0.6f;
	abund.Depletion[8] = .3f;
	abund.Depletion[9] = 1.f;
	abund.Depletion[10] = 0.2f;
	abund.Depletion[11] = 0.2f;
	abund.Depletion[12] = 0.01f;
	abund.Depletion[13] = 0.03f;
	abund.Depletion[14] = .25f;
	abund.Depletion[15] = 1.0f;
	abund.Depletion[16] = 0.4f;
	abund.Depletion[17] = 1.0f;
	abund.Depletion[18] = .3f;
	abund.Depletion[19] = 1e-4f;
	abund.Depletion[20] = 5e-3f;
	abund.Depletion[21] = 8e-3f;
	abund.Depletion[22] = 6e-3f;
	abund.Depletion[23] = 6e-3f;
	abund.Depletion[24] = 5e-2f;
	abund.Depletion[25] = 0.01f;
	abund.Depletion[26] = 0.01f;
	abund.Depletion[27] = 0.01f;
	abund.Depletion[28] = .1f;
	abund.Depletion[29] = .25f;

	abund.lgDepln = FALSE;
	abund.ScaleMetals = 1.;

	/* set true with no auger command, normally false */
	yield.lgKillAuger = FALSE;

	/* for totally forbidden He lines this will be the A */
	/* >>chng 04 may 17, change from 1e-20 to 1e-30 to do high n rydberg levels */
	iso.SmallA = 1e-30f;

	/* set parameter for atoms done in unified isoelectronic sequences */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* number of Lyman lines to include in opacities, this can be vastly larger
		 * than the number of actual levels in the model atom */
		iso.nLyman[ipISO] = 100;

		/* controls whether l-mixing and collisional ionization included */
		iso.lgColl_l_mixing[ipISO] = TRUE;
		iso.lgColl_excite[ipISO] = TRUE;
		iso.lgColl_ionize[ipISO] = TRUE;

		/* this will become a check on how well case b worked */
		for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
		{
			/* option to print departure coefficient */
			iso.lgPrtDepartCoef[ipISO][nelem] = FALSE;
			/* print hydrogenic level populations */
			iso.lgPrtLevelPops[ipISO][nelem] = FALSE;
			iso.CaseBCheck[ipISO][nelem] = -FLT_MAX;
			iso.TwoNu_induc_dn_max[ipISO][nelem] = 0.;
			/* a first guess at the recom coefficients */
			iso.RadRec_caseB[ipISO][nelem] = 1e-13;
			iso.lgFullSize[ipISO][nelem] = FALSE;
		}
	}
	/* reset with SET IND2 command, turns on/off induced two photon */
	iso.lgInd2nu_On = FALSE;

	/* hydrogen redistribution functions */
	iso.ipLyaRedist[ipH_LIKE] = ipLY_A;
	iso.ipResoRedist[ipH_LIKE] = ipCRD;
	iso.ipSubRedist[ipH_LIKE] = ipCRDW;

	/* this is the upper level for each Lya, which uses the special ipLY_A */
	iso.nLyaLevel[ipH_LIKE] = ipH2p;
	iso.nLyaLevel[ipHE_LIKE] = ipHe2p1P;

	/* he-like redistribution functions */
	/*TODO	2	why are higher lyman lines (ipResoRedist) not same as
	 * for H-like? */
	iso.ipLyaRedist[ipHE_LIKE] = ipPRD;
	iso.ipResoRedist[ipHE_LIKE] = ipPRD;
	iso.ipSubRedist[ipHE_LIKE] = ipCRDW;

	/* option to turn off an element */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		/* set of scale factors for changing abundances with elements command */
		abund.ScaleElement[nelem] = 1.;
		/* option to have abundances from table */
		abund.lgAbunTabl[nelem] = FALSE;
	}

	rt.dTauMase = 0.;
	rt.lgMaserCapHit = FALSE;
	rt.lgMaserSetDR = FALSE;

	rt.DoubleTau = 1.;
	rt.lgFstOn = TRUE;
	rt.lgElecScatEscape = TRUE;

	/* there was a call to TestCode */
	lgTestCodeCalled = FALSE;
	/* test code enabled with set test command */
	lgTestCodeEnabled = FALSE;

	/* zero out some grain variables */
	GrainZero();

	/* following should never be set non-zero  */
	FeII.fe21406 = 0.;
	FeII.fe21507 = 0.;
	FeII.fe21508 = 0.;
	FeII.fe21609 = 0.;
	FeII.fe21308 = 0.;
	FeII.fe21207 = 0.;
	FeII.fe21106 = 0.;
	FeII.fe21006 = 0.;
	FeII.fe21204 = 0.;
	FeII.fe21103 = 0.;
	FeII.fe21104 = 0.;
	FeII.fe21001 = 0.;
	FeII.fe21002 = 0.;
	FeII.fe20201 = 0.;
	FeII.fe20302 = 0.;
	FeII.fe20706 = 0.;
	FeII.fe20807 = 0.;
	FeII.fe20908 = 0.;
	FeII.fe21007 = 0.;
	FeII.fe21107 = 0.;
	FeII.fe21108 = 0.;
	FeII.fe21110 = 0.;
	FeII.fe21208 = 0.;
	FeII.fe21209 = 0.;
	FeII.fe21211 = 0.;

	fe.Fe4CoolTot = 0.;
	fe.fe40401 = 0.;
	fe.fe42836 = 0.;
	fe.fe42829 = 0.;
	fe.fe42567 = 0.;
	fe.fe41207 = 0.;
	fe.fe41206 = 0.;
	fe.fe41106 = 0.;
	fe.fe41007 = 0.;
	fe.fe41008 = 0.;
	fe.fe40906 = 0.;

	fe.Fe7CoolTot = 0.;
#	if 0
	fe.Fe7_6087 = 0.;
	fe.Fe7_5721 = 0.;
	fe.Fe7_6601 = 0.;
	fe.Fe7_3760 = 0.;
	fe.Fe7_3588 = 0.;
#	endif

	/* this is flag saying whether this is very first call,
	 * a time when space has not been allocated */
	lgFirstCall = FALSE;

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

/*rfield_opac_zero zero out rfield arrays between certain limits */
void rfield_opac_zero( long lo , long ihi )
{
	long int i;

	/* >>chng 01 aug 19, space not allocated yet,
	 * following code must also be present in contcreatemesh where
	 * space allocated for the first time */
	if( lgRfieldMalloced )
	{
		unsigned long n=(unsigned long)(ihi-lo+1);
		memset(&rfield.OccNumbDiffCont[lo]      , 0 , n*sizeof(float) );
		memset(&rfield.OccNumbContEmitOut[lo]   , 0 , n*sizeof(float) );
		memset(&rfield.ContBoltz[lo]            , 0 , n*sizeof(double) );
		memset(&rfield.ConEmitLocal[lo]         , 0 , n*sizeof(float) );
		memset(&rfield.ConEmitReflec[lo]        , 0 , n*sizeof(float) );
		memset(&rfield.ConEmitOut[lo]           , 0 , n*sizeof(float) );
		memset(&rfield.reflin[lo]               , 0 , n*sizeof(float) );
		memset(&rfield.ConRefIncid[lo]          , 0 , n*sizeof(float) );
		memset(&rfield.SummedCon[lo]            , 0 , n*sizeof(float) );
		memset(&rfield.OccNumbBremsCont[lo]     , 0 , n*sizeof(float) );
		memset(&rfield.convoc[lo]               , 0 , n*sizeof(float) );
		memset(&rfield.flux[lo]                 , 0 , n*sizeof(float) );
		memset(&rfield.SummedOcc[lo]            , 0 , n*sizeof(float) );
		memset(&rfield.SummedDif[lo]            , 0 , n*sizeof(float) );
		memset(&rfield.flux_accum[lo]           , 0 , n*sizeof(float) );
		memset(&rfield.otslin[lo]               , 0 , n*sizeof(float) );
		memset(&rfield.otscon[lo]               , 0 , n*sizeof(float) );
		memset(&rfield.ConInterOut[lo]          , 0 , n*sizeof(float) );
		memset(&rfield.outlin[lo]               , 0 , n*sizeof(float) );
		memset(&rfield.outlin_noplot[lo]        , 0 , n*sizeof(float) );
		memset(&rfield.ConOTS_local_OTS_rate[lo], 0 , n*sizeof(float) );
		memset(&rfield.ConOTS_local_photons[lo] , 0 , n*sizeof(float) );
		memset(&rfield.otsconNew[lo]            , 0 , n*sizeof(float) );
		memset(&rfield.otslinNew[lo]            , 0 , n*sizeof(float) );
		memset(&opac.OldOpacSave[lo]            , 0 , n*sizeof(double) );
		memset(&opac.opacity_abs[lo]            , 0 , n*sizeof(double) );
		memset(&opac.opacity_sct[lo]            , 0 , n*sizeof(double) );
		memset(&opac.albedo[lo]                 , 0 , n*sizeof(double) );
		memset(&opac.FreeFreeOpacity[lo]        , 0 , n*sizeof(double) );

		for( i=lo; i <= ihi; i++ )
		{
			opac.TauTotalGeo[0][i] = opac.taumin;
			opac.TauAbsGeo[0][i] = opac.taumin;
			opac.TauScatGeo[0][i] = opac.taumin;
			opac.tmn[i] = 1.;
			opac.ExpZone[i] = 1.;
			opac.e2TauAbs[i] = 1.;
			opac.ExpmTau[i] = 1.;
			opac.OpacStatic[i] = 1.;
		}
		/* also zero out fine opacity fine grid fine mesh array */
		memset(rfield.fine_opac , 0 , (unsigned long)rfield.nfine_malloc*sizeof(float) );
		/* also zero out fine opacity array */
		memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(float) );
	}
	return;
}


