/*zero zero out or initialize variables, called by cdInit, but also by func during optimization */
#include <string.h>
#include <float.h>
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "nfeii.h"
#include "incidspec.h"
#include "pnunit.h"
#include "o3exc.h"
#include "dtmase.h"
#include "plwidth.h"
#include "mgexc.h"
#include "testit.h"
#include "grainvar.h"
#include "hctmin.h"
#include "abuntabllg.h"
#include "punpoint.h"
#include "hhe2phtots.h"
#include "chmax.h"
#include "pmp2s.h"
#include "mean.h"
#include "vrrfit.h"
#include "pop371.h"
#include "elecden.h"
#include "showme.h"
#include "codish.h"
#include "zonecnt.h"
#include "timesc.h"
#include "heots.h"
#include "peimbt.h"
#include "struc.h"
#include "dclgas.h"
#include "htopoff.h"
#include "grainemis.h"
#include "lyaext.h"
#include "photrate.h"
#include "phfiton.h"
#include "fe2neg.h"
#include "punchskip.h"
#include "ionfracs.h"
#include "prnline.h"
#include "typmatrx.h"
#include "hmollte.h"
#include "lamase.h"
#include "holod.h"
#include "prtmaser.h"
#include "insanity.h"
#include "trimstg.h"
#include "fluct.h"
#include "ioreco.h"
#include "supdie.h"
#include "dtrans.h"
#include "dheton.h"
#include "ipsolar.h"
#include "elmton.h"
#include "diffheav.h"
#include "scalem.h"
#include "physok.h"
#include "h2opac.h"
#include "double.h"
#include "hmidep.h"
#include "ca2mly.h"
#include "pheat.h"
#include "hevmolec.h"
#include "lotsco.h"
#include "pca2ex.h"
#include "h2ozer.h"
#include "turb.h"
#include "gionrc.h"
#include "fntset.h"
#include "jmin.h"
#include "drneg.h"
#include "coatom.h"
#include "difher.h"
#include "comneg.h"
#include "cextra.h"
#include "hctonoff.h"
#include "globul.h"
#include "input.h"
#include "bit32.h"
#include "negoi.h"
#include "pressure.h"
#include "dndt.h"
#include "intalk.h"
#include "sumheat.h"
#include "tff.h"
#include "h2max.h"
#include "stopcalc.h"
#include "numderiv.h"
#include "wind.h"
#include "e2tau.h"
#include "exptau.h"
#include "grmetl.h"
#include "negdrg.h"
#include "cmfrac.h"
#include "gaunts.h"
#include "colden.h"
#include "stimax.h"
#include "hmi.h"
#include "hmihet.h"
#include "dcala.h"
#include "tcool.h"
#include "heat.h"
#include "hgamsc.h"
#include "bounds.h"
#include "pmpcah.h"
#include "ffsum.h"
#include "tehigh.h"
#include "hphoto.h"
#include "nsbig.h"
#include "dalerr.h"
#include "rfield.h"
#include "drminu.h"
#include "uspher.h"
#include "hemase.h"
#include "bndsok.h"
#include "habing.h"
#include "resion.h"
#include "thlo.h"
#include "the1lo.h"
#include "he3lines.h"
#include "charexc.h"
#include "taumin.h"
#include "taucon.h"
#include "he.h"
#include "touton.h"
#include "solars.h"
#include "solar.h"
#include "deplon.h"
#include "radius.h"
#include "opac.h"
#include "neutrn.h"
#include "fe2tau.h"
#include "tgff.h"
#include "broke.h"
#include "nustbl.h"
#include "fluors.h"
#include "max2pht.h"
#include "facexp.h"
#include "sumoutinc.h"
#include "angllum.h"
#include "secondaries.h"
#include "called.h"
#include "phycon.h"
#include "aaaa.h"
#include "warnings.h"
#include "colzro.h"
#include "nhe1lvl.h"
#include "nhe3lvl.h"
#include "he1bn.h"
#include "he3bn.h"
#include "he3n.h"
#include "he3pht.h"
#include "he1dmp.h"
#include "fe2cool.h"
#include "fe4cool.h"
#include "zero.h"

void zero(void)
{
	long int i,
	  j, 
	  nd;

#	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.  It is called before
	 * the boundary conditions are read in, and is never called again
	 * during that calculation.  When the code is used as a subroutine
	 * this routine is called first, when a specific model is initialized.
	 * All default variables that must be initialized before a calculation
	 * must appear in the routine.
	 */

	/* set version number */
	AAAA();

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

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

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

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

	dclgas.DustCoolGas = 0.;
	PunPoint.lgPunPoint = FALSE;

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

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

	/* change amgle of illumination */
	Angllum.AngleIllum = 1.;

	for( i=0; i < NZLIM; i++ )
	{
		struc.testr[i] = 0.;
		struc.volstr[i] = 0.;
		struc.radstr[i] = 0.;
		struc.histr[i] = 0.;
		struc.hiistr[i] = 0.;
		struc.ednstr[i] = 0.;
		struc.o3str[i] = 0.;
		struc.heatstr[i] = 0.;
		struc.coolstr[i] = 0.;
		struc.pdenstr[i] = 0.;
	}

	holod.fcool = 0.;
	holod.dfcool = 0.;
	holod.feheat = 0.;
	pmp2s.p1909 = 0.;
	pmp2s.p2326 = 0.;
	pmp2s.p1895 = 0.;
	pmp2s.p1666 = 0.;
	pmp2s.p1401 = 0.;

	/* this is the default allowed relative error in the electron density */
	ElecDen.EdenError = 0.01f;
	/* extra electron density, set with eden command */
	ElecDen.EdenExtra = 0.;
	/* forced electron density, set with set eden command */
	ElecDen.EdenSet = 0.;

	fluct.flong = 0.;
	fluct.flcPhase = 0.;
	HCTOnOff.HCTOn = 0.;

	/* lowest level is 2s,
	 * in versions 90 this was equal to 10, but was 1 in 91
	 * >>chng 99 jan 16, changed back to 10*/
	/*HTopOff.nHTopOff = 1;*/
	HTopOff.nHTopOff = 10;

	/* 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( HTopOff.chHTopType, "scal" );*/
	strcpy( HTopOff.chHTopType, " add" );

	PhotRate.PhotScaleOn = 1.;
	vrrfit.lgVrrFit = TRUE;
	PhFitOn.lgPhFit = FALSE;


	/* age of the cloud, to check for time-steady */
	timesc.CloudAgeSet = -1.f;
	/* some timescale for CO and H2 */
	timesc.AgeH2MoleDest = 0.;
	timesc.AgeCOMoleDest = 0.;
	timesc.BigH2MoleForm = 0.;
	timesc.BigCOMoleForm = 0.;

	peimbt.tsqden = 1e7;
	codish.CODissHeat = 0.;
	PrtMaser.lgPrtMaser = FALSE;
	NumDeriv.lgNumDeriv = FALSE;
	fluct.lgDenFluc = TRUE;

	colzro();

	for( i=0; i < LIMPUN; i++ )
	{
		pnunit.lgPunLstIter[i] = FALSE;
	}

	/* normally is unity, set to zero when hydrogen collisions turned off,
	 * collisional ionization*/
	hydro.HColIonOn = 1.;
	/* collisional excitation */
	hydro.HColExcOn = 1.;
	/* 2s - 2p exchange collisions */
	hydro.lgC2psOn = TRUE;
	/* free free heating, cooling, net */
	hydro.FreeFreeCool = 0.;
	hydro.FreeFreeHeat = 0.;
	hydro.HFFNet = 0.;
	/* option to print emissivity instead of intensity/luminosity */
	hydro.lgHydEmiss = FALSE;

	double_.DoubleTau = 1.;
	SumOutInc.SumOutCon = 0.;
	SumOutInc.SumIncCon = 0.;

	/* which matrix inversion routine to use */
	strcpy( TypMatrx.chMatrix, "linpack" );

	PrnLine.lgPrtCont = FALSE;
	PrnLine.lgPrnDiff = FALSE;
	PrnLine.lgPrnPump = FALSE;
	PrnLine.lgPrnInwd = FALSE;
	PrnLine.lgPrnColl = FALSE;
	PrnLine.lgPrnHeat = FALSE;

	h2opac.H2Opacity = 0.;

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

	Secondaries.heatef = 1.;
	Secondaries.efionz = 0.;
	Secondaries.exctef = 0.;
	Secondaries.csupra = 0.;
	Secondaries.x12tot = 0.;
	Secondaries.HyExc = 0.;
	Secondaries.HyIon = 0.;
	Secondaries.seccmp = 0.;
	Secondaries.scmpla = 0.;
	Secondaries.seche = 0.;
	Secondaries.helax = 0.;
	Secondaries.sec2total = 0.;

	Secondaries.Hx12[0][LIMELM] = 0.;
	Secondaries.Hx12[1][LIMELM] = 0.;
	for( i=0; i< LIMELM; ++i )
	{
		Secondaries.Hx12[0][i] = 0.;
		Secondaries.Hx12[1][i] = 0.;
		for( j=0; j<LIMELM; ++j )
		{
			Secondaries.SecondaryExcite[i][j] = 0.;
			Secondaries.SecondaryIonize[i][j] = 0.;
		}
	}

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

	/* lte populations of some hydrogen molecules */
	HMolLTE.phmlte = 0.;
	HMolLTE.ph2lte = 0.;
	HMolLTE.phplte = 0.;

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

	Max2Pht.xMax2Pht = 0.;
	Fe2Neg.nFe2Neg = 0;

	/* turn off printing of heating agents */
	pheat.lgPrintHeating = FALSE;

	/* turn off flag for insane goings */
	insanity.lgInsane = FALSE;

	/* punch every one continuum point */
	PunchSkip.ncPunchSkip = 1;

	/* flag for Lya masing */
	LaMase.xLaMase = 0.;

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

	/* default condition is burgess supressed, nusbaumer and storye not */
	supdie.lgSupDie[0] = TRUE;
	supdie.lgSupDie[1] = 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 */
	ipSolarCom.npSolar = LIMELM - 1;
	for( i=0; i < ipSolarCom.npSolar; i++ )
	{
		ipSolarCom.ipSolar[i] = i + 2;
	}

	/* option to turn off an element */
	for( i=0; i < LIMELM; i++ )
	{
		difher.CaseBCheck[i] = -FLT_MAX;
		elmton.lgElmtOn[i] = TRUE;
		/* set of scale factors for changing abundances with elements command */
		scalem.ScaleElement[i] = 1.;
		/* option to have abundances from table */
		AbunTablLG.lgAbunTabl[i] = FALSE;
	}

	AbunTablLG.lgAbTaON = FALSE;

	/* zero out diffuse recombinaiton emission arrays */
	for( i=0; i < LIMELM; i++ )
	{
		for( j=0; j < (LIMELM - 1); j++ )
		{
			DiffHeav.DiffusHeavy[j][i] = 0.;
			DiffHeav.xLyaHeavy[j][i] = 0.;
		}
	}

	for( i=0; i < LIMELM; i++ )
	{
		for( j=0; j < LIMELM; j++ )
		{
			for( nd=0; nd < 7; nd++ )
			{
				/* must be zero since ion routines use these when
				 * not yet defined */
				PhotRate.PhotoRate[0][nd][j][i] = 0.;
				PhotRate.PhotoRate[1][nd][j][i] = 0.;
				PhotRate.PhotoRate[2][nd][j][i] = 0.;
			}
		}
	}

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

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

	fluors.lgFluor = TRUE;
	dheton.lgDHetOn = TRUE;
	dheton.lgDColOn = TRUE;
	physok.lgPhysOK = TRUE;

	hmidepCom.hmidep = 1.;
	hmidepCom.h2dep = 1.;
	hmidepCom.h2pdep = 1.;
	hmidepCom.h3pdep = 1.;
	ioRecom = NULL;
	lgioRecom = FALSE;


	/* 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;
	wind.dVeldRad = 0.;

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

	/* main arrays to save ionization fractions*/
	for( i=0; i < LIMELM; i++ )
	{
		for( j=0; j < (LIMELM + 1); j++ )
		{
			xIonFracs[j][i] = 0.;
		}
	}

	ca2mly.Ca2RmLya = 0.;
	pca2ex.popca2ex = 0.;

	gionrc.grecon = 1.;
	cmfrac.CarMolFrac = -1.f;
	lotsco.lgLotsCO = FALSE;
	h2ozer.lgH2Ozer = FALSE;
	fntset.lgFntSet = FALSE;
	grmetl.GrainMetal = 1.;
	comneg.lgComNeg = FALSE;
	coatom.catmic = 0.;
	coatom.oatmic = 0.;

	/* zero out molecular species */
	for( i=0; i < NHEVML; i++ )
	{
		hevmolec.hevmol[i] = 0.;
		hevmolec.hevcol[i] = 0.;
	}

	jmin.rjnmin = FLT_MAX;
	jmin.ajmmin = FLT_MAX;
	drneg.lgDrNeg = FALSE;

	/* two photon absorption deep in cloud */
	hhe2phtots.H12phtOTS = 0.;
	hhe2phtots.He12phtOTS = 0.;
	hhe2phtots.He22phtOTS = 0.;

	/* check whether this is a 32-bit machine
	 * use te in common, which forc will translate as a float.
	 * then set te after test */
	phycon.te = 1e-37f;
	phycon.te /= 1e40;
	/* had to move above up here to prevent ms c from using double without
	 * doing store*/

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

	strcpy( dtrans.chDffTrns, "OU2" );
	turb.TurbVel = 0.;
	cextra.lgCExtraOn = FALSE;
	cextra.CoolExtra = 0.;

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

	pressure.lgRadPresON = FALSE;
	pressure.RadPres = 0.;
	pressure.lgPradCap = FALSE;
	pressure.lgPradDen = FALSE;
	pressure.lgRadPrsON = 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;
	strcpy( pressure.chCPres, "CDEN" );

	for( i=0; i < 20; i++ )
	{
		pressure.plines[i] = 0.;
	}

	pressure.RadBetaMax = 0.;
	pressure.PresMax = 0.;

	/* check whether this is a 32-bit machine
	 * use te in common, which forc will translate as a float.
	 * then set te after test */
	if( phycon.te == 0. )
	{
		bit32.lgBit32 = TRUE;
	}
	else
	{
		bit32.lgBit32 = FALSE;
	}

	/* now set physical conditions array */
	phycon.te = 1e4;
	phycon.eden = -99.f;
	phycon.hden = -99.f;

	globul.glbdst = 0.;
	globul.glbrad = 0.;
	negoi.nNegOI = 0;

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

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

	for( i=0; i < NCELL; i++ )
	{
		opac.OldOpacSave[i] = 0.;
		/* chng>>> 99 apr 11, only opac had been set to 0, eovrlp crashed
		 * with div by zero when called by tauout before opacities had been
		 * defined.  Now set al three to very small number,
		 * CODE CRASHED - changed back to zero, apparetly some part of
		 * code uses opacity of zero to check on continuum bounds */
		opac.opac[i] = 0.;
		/* NB read above before changing set to zero in above */
		opac.scatop[i] = 0.;
		opac.albedo[i] = 0.;
		opac.TauTotal[0][i] = tauminCom.taumin;
		opac.tauabs[0][i] = tauminCom.taumin;
		opac.tausct[0][i] = tauminCom.taumin;
		opac.tmn[i] = 1.;
		opac.FreeFreeOpacity[i] = 0.;

		rfield.DiffLocal[i] = 0.;
		rfield.ConOutNoInter[i] = 0.;
		rfield.ConRefNoInter[i] = 0.;
		rfield.ContBoltz[i] = 0.;

		rfield.reflin[i] = 0.;
		rfield.reflux[i] = 0.;
		rfield.refcon[i] = 0.;
		facexp.ExpZone[i] = 1.;
		e2tau.e2TauAbs[i] = 1.;
		exptau.ExpmTau[i] = 1.;
		gaunts.gff[i] = 1.1f;
		gaunts.gffhe2[i] = 1.1f;
		GrainEmis.GrainEmission[i] = 0.;
	}

	/* 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 */
	opac.lgUseFileOpac = TRUE;

	rfield.nflux = NCELL;

	/* these are flags to let gff sub know to recompute ContBoltz and gss */
	tgff.tgffused = -1.f;
	tgff.tgffsued2 = -1.f;

	dndt.dDensityDT = 0.;
	HPhoto.cintot = 0.;
	resion.otsfe2 = 0.;
	tffCom.ntff = 1;
	tffCom.tff = 0.;
	h2max.BiggestH2 = -1.f;
	negdrg.lgNegGrnDrg = FALSE;
	nustbl.nUnstable = 0;
	nustbl.lgUnstable = FALSE;

	/* various criteria for stopping model */
	StopCalc.tauend = 0.;
	StopCalc.taunu = 0.;
	StopCalc.iptnu = NCELL;
	/* highest allowed temperature */
	StopCalc.T2High = 1e12f;
	/* highest and lowest temperatues to ever allow */
	StopCalc.TeHighest = 1e10f;
	StopCalc.TeLowest = 2.8f;
	/* stopping temperature */
	StopCalc.tend = 4000.;
	/* stop electron fraction, usually not used */
	StopCalc.StopElecFrac = -FLT_MAX;
	/* ending column densities */
	StopCalc.HColStop = 1e30f;
	StopCalc.colpls = 1e30f;
	StopCalc.colnut = 1e30f;
	StopCalc.StopElecDensity = -1e30f;

	/* limits to stop line command */
	StopCalc.nstpl = 0;
	for( i=0; i < MXSTPL; i++ )
	{
		StopCalc.stpint[i] = FLT_MAX;
		StopCalc.ipStopLin1[i] = 1;
		StopCalc.ipStopLin2[i] = 1;
		StopCalc.LineStopWl[i] = 0;
	}

	neutrn.frcneu = 0.;
	neutrn.effneu = 1.;
	neutrn.totneu = 0.;
	neutrn.lgNeutrnHeatOn = FALSE;

	pnunit.npunch = 0;
	hemase.lgHeMase = FALSE;
	uspher.lgUSphON = FALSE;
	stimaxCom.stimax[0] = 0.;
	stimaxCom.stimax[1] = 0.;
	hmi.htwo = 0.;
	hmi.h2plus = 0.;
	hmi.h3plus = 0.;
	hmi.hminus = 0.;
	dcala.oldcla = 0.;
	dcala.Ca3d = 0.;
	dcala.Ca4p = 0.;
	dcala.dstCala = 0.;
	hmihetCom.hmihet = 0.;
	hmihetCom.hmitot = 0.;

	for( i=0; i < NCOLD; i++ )
	{
		coldenCom.colden[i] = 0.;
	}

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

	for( i=0; i < LIMSPC; i++ )
	{
		/* range(i,1) = 1.
		 * >>chng 96 dec 18, from inf mass ryd to H mass ryg */
		IncidSpec.range[0][i] = (float)HIONPOT;
		IncidSpec.range[1][i] = bounds.egamry;
	}

	tcool.totcol = 0.;
	tcool.GrossHeat = 0.;
	tcool.heatl = 0.;
	tcool.coolheat = 0.;

	heat.HeatNet = 0.;
	heat.htot = 1.;

	tcool.ctot = 1.;
	heat.power = 0.;

	radius.rinner = 0.e0;
	radius.Radius = 0.e0;
	radius.drad = 0.e0;
	radius.depth = 1e-30;
	radius.r1r0sq = 1.;
	/* this is changed with the roberto command, to go from out to in */
	radius.dRadSign = 1.;
	/* RDFALT is log of default starting radius (cm) */
	radius.rdfalt = 25.;

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

	radius.dReff = 1.;
	radius.dVeff = 1.;
	radius.dRNeff = 1.;
	radius.lgdR2Small = FALSE;

	ZoneCnt.nzone = 0;
	ZoneCnt.nprint = 1000;
	ZoneCnt.lgZoneSet = FALSE;
	ZoneCnt.lgZoneTrp = FALSE;
	ZoneCnt.lgEndDflt = TRUE;

	/* this is default number of zones
	 * >>chng 96 june 5, from 400 to 500 for thickes corners4 grid */
	ZoneCnt.nEndDflt = 600;

	for( i=0; i < ITRDIM; i++ )
	{
		ZoneCnt.nend[i] = ZoneCnt.nEndDflt;
		radius.router[i] = 1e30;
	}

	hgamsc.BoundCompton = 0.;
	heots.esc584 = 0.;
	heots.esc626 = 0.;
	heots.esc911 = 0.;
	heots.esc910 = 0.;
	heots.esc610 = 0.;
	heots.esc660 = 0.;

	/* save initial condition for talk in case PRINT LAST used */
	intalk.lgInTalk = called.lgTalk;

	/* flags for whether continuum is defined over all energies */
	bndsok.lgMMok = TRUE;
	bndsok.lgHPhtOK = TRUE;
	bndsok.lgXRayOK = TRUE;
	bndsok.lgGamrOK = TRUE;
	hydro.lgHiPop2 = FALSE;
	hydro.pop2mx = 0.;
	pmpcah.escah = 1.40e8;
	ffsum.EdenFFSum = 0.;
	tehigh.thist = 0.;
	tehigh.tlowst = 1e20f;
	nsbigCom.nsbig = 0;
	dalerr.ilt = 0;
	dalerr.iltln = 0;
	dalerr.ilthn = 0;
	dalerr.ihthn = 0;
	dalerr.ifail = 0;

	habing.lgHabing = FALSE;

	drminu.lgDrMinUsed = FALSE;

	resion.otsfe2 = 0.;

	mgexc.popmg2 = 0.;
	mgexc.rmg2l = 0.;
	o3exc.poiii2 = 0.;
	o3exc.poiii3 = 0.;
	o3exc.poiexc = 0.;

	thlo.HydTempLimit = 1000.;

	the1loCom.the1lo = 1000.;
	the1loCom.nhe1lo = 0;

	he3lines.pht2s = 0.;
	he3lines.p2s = 0.;
	he3lines.p2p = 0.;
	he3lines.p3s = 0.;
	he3lines.p3p = 0.;
	he3lines.p3d = 0.;
	he3lines.p10830 = 0.;
	he3lines.p5876 = 0.;
	he3lines.p7065 = 0.;
	he3lines.p3889 = 0.;
	he3lines.he3clx = 0.;
	o3exc.d5007r = 0.;
	o3exc.d5007t = 0.;
	o3exc.d4363 = 0.;
	o3exc.d6300 = 0.;

	/***************************************************
	 * charge transfer ionization and recombination 
	 ***************************************************/
	for( i=0; i < 2; i++ )
	{
		for( j=0; j < 2; ++j )
		{
			CharExc.CTHion[i][j] = 0.;
			CharExc.CTHrec[i][j] = 0.;
		}
	}
	for( i=0; i<5; ++i )
	{
		CharExc.CTHionAgnt[i] = 0.;
	}
	CharExc.chrexCool = 0.;
	CharExc.chrener = 0.;

	/* HCharHeatMax, HCharCoolMax are largest fractions of heating in cur zone
	 * or cooling due to ct */
	CharExc.HCharHeatMax = 0.;
	CharExc.HCharCoolMax = 0.;

	/* flag saying that charge transfer heating should be included,
	 * turned off with no CTHeat commmand */
	CharExc.HCharHeatOn = 1.;
	for( i=0; i< LIMELM; ++i )
	{
		for( j=0; j<LIMELM; ++j )
		{
			CharExc.HCharExcIon[i][j] = 0.;
			CharExc.HCharExcRec[i][j] = 0.;
		}
	}

	taucon.telec = tauminCom.taumin;
	taucon.thmin = tauminCom.taumin;

	pmpcah.tpcah[0] = tauminCom.taumin;
	pmpcah.tpcah[1] = 1e20f;

	he.hei = 0.;
	he.heii = 0.;
	he.heiii = 0.;
	he.hei3 = 0.;
	he.hei1 = 0.;

	touton.lgTauOutOn = FALSE;

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

	deplon.lgDepln = FALSE;

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

	dtMase.dTauMase = 0.;
	dtMase.lgMaserCapHit = FALSE;

	Testit.lgTestIt = FALSE;
	Testit.lgTestOn = FALSE;
	CHMax.CoolHeatMax = 0.;
	CHMax.ipCoolHeatMax = 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 */
	HCTMin.HCTAlex = 1.92e-9f;

	/* some grain variables */
	for( nd=0; nd < NDUST; nd++ )
	{
		GrainVar.lgQHeat[nd] = FALSE;
		GrainVar.lgDustVary[nd] = FALSE;
		GrainVar.DustDftVel[nd] = FLT_MAX;
		GrainVar.avdft[nd] = 0.;
		GrainVar.TeGrainMax[nd] = FLT_MAX;
		GrainVar.lgDustOn1[nd] = FALSE;
		GrainVar.dheat[nd] = FLT_MAX;
		GrainVar.dstfactor[nd] = 1.;
		/* this must be zero for the first solutions to be able to converge */
		GrainVar.tedust[nd] = 0.;
		GrainVar.lgQHeat[nd] = FALSE;
		GrainVar.nqp1[nd] = 0;
		GrainVar.qdep[nd] = 1.0;
		GrainVar.qhnrat[nd] = 1.e7;
	}
	GrainVar.lgDustOn = FALSE;

	GrainVar.nqout1 = 0;
	GrainVar.iqp = NULL;
	GrainVar.qtmax = 0.;
	GrainVar.qhnrat[0] = 1.8732e10f;
	GrainVar.qhnrat[1] = 1.6387e10f;
	GrainVar.qhnrat[2] = 2.7518e11f;
	GrainVar.qhnrat[3] = 2.4073e11f;
	GrainVar.qhnrat[4] = 3.7065e09f;
	GrainVar.qhnrat[5] = 3.7069e11f;
	GrainVar.qhnrat[6] = 1.6387e10f;
	GrainVar.qhnrat[7] = 1.6387e10f;
	GrainVar.qhnrat[8] = 3.6312e07f;
	GrainVar.qhnrat[9] = 6.4191e06f;
	he3pht.he3photo = 0.;
	he3pht.he3lya = 0.;

	for( i=0; i < (NHE1LVL + 1); i++ )
	{
		he1bnCOM.he1bn[i] = 1.;
	}

	/* array of lifetimes for helium lines - sum of a's down/4pi*912 */
	he1dmpCOM.he1dmp[0] = 0.;
	he1dmpCOM.he1dmp[1] = 455.1f;
	he1dmpCOM.he1dmp[2] = 72.4f;
	he1dmpCOM.he1dmp[3] = 21.9f;
	he1dmpCOM.he1dmp[4] = 8.42f;
	he1dmpCOM.he1dmp[5] = 3.77f;
	he1dmpCOM.he1dmp[6] = 1.82f;
	he1dmpCOM.he1dmp[7] = 0.939f;
	he1dmpCOM.he1dmp[8] = 0.438f;

	for( i=0; i < NHE3LVL; i++ )
	{
		he3nCom.he3n[i] = 0.;
		he3bnCom.he3bn[i] = 1.;
	}

	fe2cool.Fe2L16Tot = 0.;
	fe2cool.fe21406 = 0.;
	fe2cool.fe21507 = 0.;
	fe2cool.fe21508 = 0.;
	fe2cool.fe21609 = 0.;
	fe2cool.fe21308 = 0.;
	fe2cool.fe21207 = 0.;
	fe2cool.fe21106 = 0.;
	fe2cool.fe21006 = 0.;
	fe2cool.fe21204 = 0.;
	fe2cool.fe21103 = 0.;
	fe2cool.fe21104 = 0.;
	fe2cool.fe21001 = 0.;
	fe2cool.fe21002 = 0.;
	fe2cool.fe20201 = 0.;
	fe2cool.fe20302 = 0.;
	fe2cool.fe20706 = 0.;
	fe2cool.fe20807 = 0.;
	fe2cool.fe20908 = 0.;
	fe2cool.fe21007 = 0.;
	fe2cool.fe21107 = 0.;
	fe2cool.fe21108 = 0.;
	fe2cool.fe21110 = 0.;
	fe2cool.fe21208 = 0.;
	fe2cool.fe21209 = 0.;
	fe2cool.fe21211 = 0.;
	fe4cool.Fe4CoolTot = 0.;
	fe4cool.fe40401 = 0.;
	fe4cool.fe42836 = 0.;
	fe4cool.fe42829 = 0.;
	fe4cool.fe42567 = 0.;
	fe4cool.fe41207 = 0.;
	fe4cool.fe41206 = 0.;
	fe4cool.fe41106 = 0.;
	fe4cool.fe41007 = 0.;
	fe4cool.fe41008 = 0.;
	fe4cool.fe40906 = 0.;

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