/* 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 */
/*IterStart set and save values of many variables at start of iteration */
/*IterRestart restart iteration */
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "iso.h"
#include "grainvar.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "dynamics.h"
#include "prt.h"
#include "hyperfine.h"
#include "dense.h"
#include "magnetic.h"
#include "continuum.h"
#include "geometry.h"
#include "h2.h"
#include "he.h"
#include "grains.h"
#include "atomfeii.h"
#include "pressure.h"
#include "stopcalc.h"
#include "converge.h"
#include "mean.h"
#include "ca.h"
#include "thermal.h"
#include "atoms.h"
#include "wind.h"
#include "jeans_mass.h"
#include "opacity.h"
#include "timesc.h"
#include "trace.h"
#include "colden.h"
#include "secondaries.h"
#include "hmi.h"
#include "radius.h"
#include "phycon.h"
#include "called.h"
#include "co.h"
#include "intalk.h"
#include "rfield.h"
#include "ionbal.h"
#include "atmdat.h"
#include "lines.h"
#include "aver.h"
#include "molcol.h"
#include "input.h"
#include "negdrg.h"
#include "rt.h"
#include "iterations.h"
/* these are the static saved variables, set in IterStart and reset in IterRestart */

static double h2plus_heat_save, HeatH2Dish_used_save, HeatH2Dexc_used_save,HeatH2Dish_TH85_save,
	HeatH2Dexc_TH85_save, HeatH2Dish_BigH2_save, HeatH2Dexc_BigH2_save, deriv_HeatH2Dexc_BigH2_save,
	hmihet_save, hmitot_save, H2_Solomon_dissoc_rate_used_save, H2_Solomon_dissoc_rate_used_H2g_save,
	H2_Solomon_dissoc_rate_used_H2s_save, H2_H2g_to_H2s_rate_used_save, H2_photodissoc_used_save,
	H2_photodissoc_BigH2_save,H2_H2g_to_H2s_rate_BigH2_save,H2_Solomon_dissoc_rate_BigH2_H2s_save,
	H2_Solomon_dissoc_rate_BigH2_H2g_save;

/* args are atomic number, ionization stage*/
static float xIonFsave[LIMELM+3][LIMELM+1];
static double HeatSave[LIMELM][LIMELM];
static float supsav[LIMELM][LIMELM],
dndtSave;
/* save values of dr and drNext */
static double drSave , drNextSave;

/* also places to save them between iterations */
static long int IonLowSave[LIMELM], 
	IonHighSave[LIMELM];

/* these are used to reset the level popls of the h and he model atoms */
/*static float hnsav[LIMELM][LMHLVL+1], */
static int lgHNSAV = FALSE;

static float  ***HOpacRatSav;

static float hmsav[N_H_MOLEC];
static double ***hnsav,
  edsav;

static float xMolFracsSave[LIMELM],
	gas_phase_save[LIMELM];

void IterStart(void)
{
	int lgHit;
	long int i, 
	  ii, 
	  ion, 
	  ipISO,
	  n ,
	  nelem;
	double fhe1nx, 
	  ratio;

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

	/* allocate two matrices if first time in this coreload */
	/*TODO	2	this does not belong here - move to where main data created */
	if( !lgHNSAV )
	{
		/* set flag so we never do this again */
		lgHNSAV = TRUE;

		if(  (HOpacRatSav = (float ***)MALLOC(sizeof(float **)*NISO ))==NULL  )
			BadMalloc();

		if(  (hnsav = (double ***)MALLOC(sizeof(double **)*NISO ))==NULL  )
			BadMalloc();

		for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
		{

			if(  (HOpacRatSav[ipISO] = (float **)MALLOC(sizeof(float *)*LIMELM ))==NULL  )
				BadMalloc();

			if(  (hnsav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ))==NULL  )
				BadMalloc();

			/* now do the second dimension */
			for( nelem=0; nelem<LIMELM; ++nelem )
			{
				if( (HOpacRatSav[ipISO][nelem] = (float *)MALLOC(sizeof(float)*(unsigned)(iso.numLevels[ipISO][nelem]+1) ) )==NULL )
					BadMalloc();

				if( (hnsav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipISO][nelem]+1) ) )==NULL )
					BadMalloc();
			}
		}
	}

	/* IterStart is called at the start of EVERY iteration */
	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " IterStart called.\n" );
	}

	/* check if this is the last iteration */
	if( iteration > iterations.itermx )
	{
		iterations.lgLastIt = TRUE;
	}
	else
	{
		iterations.lgLastIt = FALSE;
	}

	/* reset some vars dealing with H2
	H2_Reset(); */

	/* flag set true in RT_line_one_tauinc if maser ever capped */
	rt.lgMaserCapHit = FALSE;

	/* zero out charge transfer heating and cooling fractions */
	atmdat.HCharHeatMax = 0.;
	atmdat.HCharCoolMax = 0.;

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

	/* save exansion rate here */
	dndtSave = dynamics.dDensityDT;

	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.;
	thermal.CoolHeatMax = 0.;
	hydro.HCollIonMax = 0.;

	atmdat.HIonFracMax = 0.;

	/* will save highest n=2p/1s ratio for hydrogen atom*/
	hydro.pop2mx = 0.;
	hydro.lgHiPop2 = FALSE;

	hydro.nLyaHot = 0;
	hydro.TLyaMax = 0.;
	hydro.xLaMase = 0.;

	/* evaluate the gas and radiation pressure */
	/* this sets values of pressure.PresTotlCurr, also calls tfidle */
	pressure.PresInteg = 0.;
	pressure.pinzon = 0.;

	PresTotCurrent();
	pressure.PresTotlInit = pressure.PresTotlCurr;


	dynamics.HeatMax = 0.;
	dynamics.CoolMax = 0.;

	/* reset these since we now have a valid solution */
	pressure.pbeta = 0.;
	pressure.RadBetaMax = 0.;
	pressure.lgPradCap = FALSE;
	pressure.lgPradDen = FALSE;
	/* this flag will say we hit the sonic point */
	pressure.lgSonicPoint = FALSE;
	pressure.lgStrongDLimbo = FALSE;

	/* define two timescales for equilibrium, Compton and thermal */
	timesc.tcmptn = 1.e0/(rfield.qtot*6.65e-25*dense.eden);
	timesc.ttherm = 1.5*dense.pden*BOLTZMANN*phycon.te/thermal.ctot;

	/* this will be total mass in computed structure */
	dense.xMassTotal = 0.;

	for( nelem=0; nelem < LIMELM; nelem++ )
	{

		if( dense.lgElmtOn[nelem] )
		{

			/* now zero out the ionic fractions */
			for( ion=0; ion < (nelem + 2); ion++ )
			{
				xIonFsave[nelem][ion] = dense.xIonDense[nelem][ion];
			}

			for( ion=0; ion < LIMELM; ion++ )
			{
				HeatSave[nelem][ion] = thermal.heating[ion][nelem];
			}
		}
	}

	/* >>chng 02 jan 28, add ipISO loop */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* >>chng 03 apr 11, from 0 to ipISO */
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				for( n=0; n < iso.numLevels[ipISO][nelem]; n++ )
				{
					HOpacRatSav[ipISO][nelem][n] = iso.ConOpacRatio[ipISO][nelem][n];
					hnsav[ipISO][nelem][n] = iso.Pop2Ion[ipISO][nelem][n];
				}
			}
			iso.TwoNu_induc_dn_max[ipISO][nelem] = 0.;
		}
	}

	rfield.ipEnergyBremsThin = 0;
	rfield.EnergyBremsThin = 0.;

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

	/* save ionization bal for second row elements */
	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		/* save molecular densities */
		xMolFracsSave[nelem] = dense.xMolecules[nelem];
		gas_phase_save[nelem] = dense.gas_phase[nelem];
		IonLowSave[nelem] = dense.IonLow[nelem];
		IonHighSave[nelem] = dense.IonHigh[nelem];
	}

	edsav = dense.eden;

	/* save molecular fractions */
	for ( i=0; i<N_H_MOLEC; i++) 
	{
		hmsav[i] = hmi.Hmolec[i];
	}
	h2plus_heat_save = hmi.h2plus_heat;
	HeatH2Dish_used_save = hmi.HeatH2Dish_used;
	HeatH2Dexc_used_save = hmi.HeatH2Dexc_used;
	HeatH2Dish_TH85_save = hmi.HeatH2Dish_TH85;
	HeatH2Dexc_TH85_save = hmi.HeatH2Dexc_TH85;
	HeatH2Dish_BigH2_save = hmi.HeatH2Dish_BigH2;
	HeatH2Dexc_BigH2_save = hmi.HeatH2Dexc_BigH2;
	deriv_HeatH2Dexc_BigH2_save = hmi.deriv_HeatH2Dexc_BigH2;
	hmihet_save = hmi.hmihet;
	hmitot_save = hmi.hmitot;
	H2_Solomon_dissoc_rate_used_save = hmi.H2_Solomon_dissoc_rate_used;
	H2_Solomon_dissoc_rate_used_H2g_save = hmi.H2_Solomon_dissoc_rate_used_H2g;
	H2_Solomon_dissoc_rate_used_H2s_save = hmi.H2_Solomon_dissoc_rate_used_H2s;
	H2_H2g_to_H2s_rate_used_save = hmi.H2_H2g_to_H2s_rate_used;
	H2_photodissoc_used_save = hmi.H2_photodissoc_used;

	H2_photodissoc_BigH2_save = hmi.H2_photodissoc_BigH2;
	H2_H2g_to_H2s_rate_BigH2_save = hmi.H2_H2g_to_H2s_rate_BigH2;
	H2_Solomon_dissoc_rate_BigH2_H2s_save = hmi.H2_Solomon_dissoc_rate_BigH2_H2s;
	H2_Solomon_dissoc_rate_BigH2_H2g_save = hmi.H2_Solomon_dissoc_rate_BigH2_H2g;

	/* save zone thicknes, and next zone thickness */
	drSave = radius.drad;
	drNextSave = radius.drNext;

	geometry.nprint = iterations.IterPrnt[iteration-1];

	jeans_mass.TotMassColl = 0.;
	jeans_mass.tmas = 0.;
	jeans_mass.wmas = 0.;

	jeans_mass.rjnmin = FLT_MAX;
	jeans_mass.ajmmin = FLT_MAX;

	thermal.nUnstable = 0;
	thermal.lgUnstable = FALSE;

	rfield.extin_mag_B_point = 0.;
	rfield.extin_mag_V_point = 0.;
	rfield.extin_mag_B_extended = 0.;
	rfield.extin_mag_V_extended = 0.;

	/* plus 1 is to zero out point where unit integration occurs */
	for( i=0; i < rfield.nupper+1; i++ )
	{
		/* diffuse fields continuum */
		/* >>chng 03 nov 08, recover SummedDif */
		rfield.SummedDifSave[i] = rfield.SummedDif[i];
		rfield.ConEmitReflec[i] = 0.;
		rfield.ConEmitOut[i] = 0.;
		rfield.ConInterOut[i] = 0.;
		/* save OTS continuum for next iteration */
		rfield.otssav[i][0] = rfield.otscon[i];
		rfield.otssav[i][1] = rfield.otslin[i];
		rfield.outlin[i] = 0.;
		rfield.outlin_noplot[i] = 0.;
		rfield.ConOTS_local_OTS_rate[i] = 0.;
		rfield.ConOTS_local_photons[i] = 0.;

		/* save initial absorption, scattering opacities for next iteration */
		opac.opacity_abs_savzon1[i] = opac.opacity_abs[i];
		opac.opacity_sct_savzon1[i] = opac.opacity_sct[i];

		/* will accumulate optical depths through cloud */
		opac.TauAbsFace[i] = opac.taumin;
		opac.TauScatFace[i] = opac.taumin;
		/* >>chng 99 dec 04, having exactly 1 zone thickness for first zone caused discontinuity
		 * for heating in very high T models in map.in.  zone 1 and 2 were 20% different,
		 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
		 * these were not here at all - changed to dr/2 */
		/* attenuation of flux by optical depths IN THIS ZONE 
		 * AngleIllum is 1/COS(theta), is usually 1, reset with illuminate command,
		 * option for illumination of slab at an angle */
		/* >>chng 04 oct 09, from drad to radius.drad_x_fillfac - include fill fac, PvH */
		opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.AngleIllum);

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

		/* e2(tau) in inward direction, up to illuminated face, this is used to get the
		 * recombination escape probability */
		opac.e2TauAbs[i] = (float)e2(opac.TauAbsFace[i],opac.ExpmTau[i]);
	}

	/* this zeroes out arrays to hold mean ionization fractions
	 * later enered by mean
	 * read out by setlim */
	MeanZero();

	/* zero out the column densities */
	for( i=0; i < NCOLD; i++ )
	{
		colden.colden[i] = 0.;
	}
	colden.He123S = 0.;
	for( i=0; i < NUM_COMOLE_CALC; i++ )
	{
		co.hevcol[i] = 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.;
	}

	/* these are some line of sight emission measures */
	colden.dlnenp = 0.;
	colden.dlnenHep = 0.;
	colden.dlnenHepp = 0.;

	/* now zero heavy element molecules */
	molcol("ZERO",ioQQQ);

	/* this will be sum of all free free heating over model */
	thermal.FreeFreeTotHeat = 0.;

	thermal.thist = 0.;
	thermal.tlowst = 1e20f;

	wind.AccelAver = 0.;
	wind.acldr = 0.;
	ionbal.ilt = 0;
	ionbal.iltln = 0;
	ionbal.ilthn = 0;
	ionbal.ihthn = 0;
	ionbal.ifail = 0;

	secondaries.SecHIonMax = 0.;
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		for( ion=0; ion<nelem+1; ++ion )
		{
			supsav[nelem][ion] = secondaries.csupra[nelem][ion];
		}
	}
	secondaries.x12sav = secondaries.x12tot;
	secondaries.hetsav = secondaries.heatef;
	secondaries.savefi = secondaries.efionz;
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
	{
		for( ion=0; ion<nelem+1 ; ++ion )
		{
			ionbal.CompRecoilHeatRateSave[nelem][ion] = ionbal.CompRecoilHeatRate[nelem][ion];
			ionbal.CompRecoilIonRateSave[nelem][ion] = ionbal.CompRecoilIonRate[nelem][ion];
		}
	}

	/* these will keep track of the number of convergence failures that occur */
	conv.nTeFail = 0;
	conv.nTotalFailures = 0;
	conv.nPreFail = 0;
	conv.failmx = 0.;
	conv.nIonFail = 0;
	conv.nPopFail = 0;
	conv.nNeFail = 0;
	conv.nGrainFail = 0;
	/* abort flag */
	conv.lgAbort = FALSE;

	/* zero out averaging routine */
	aver("zero",0.,0.,"          ");

	GrainStartIter();

	rfield.comtot = 0.;

	co.codfrc = 0.;
	co.codtot = 0.;

	co.COCoolBigFrac = 0.;
	co.lgCOCoolCaped = FALSE;

	hmi.HeatH2DexcMax = 0.;
	hmi.CoolH2DexcMax = 0.;
	hmi.h2dfrc = 0.;
	hmi.h2line_cool_frac = 0.;
	hmi.h2dtot = 0.;
	thermal.HeatLineMax = 0.;
	thermal.GBarMax = 0.;
	hyperfine.cooling_max = 0.;
	hydro.cintot = 0.;
	geometry.lgZoneTrp = FALSE;

	gv.TotalDustHeat = 0.;
	gv.dphmax = 0.;
	hmi.h2pmax = 0.;
	gv.dclmax = 0.;
	gv.GrnElecDonateMax = 0.;
	gv.GrnElecHoldMax = 0.;

	/************************************************************************
	 *
	 * allocate space for lines arrays 
	 *
	 ************************************************************************/

	/* set up line array */
	/* first count number of lines */
	LineSave.ipass = -1;
	lines();
	/* in a grid this may not be the first time here, return old memory
	 * and grab new appropriate for this size, since number of lines to
	 * be stored can change
	 * as now coded this memory is freed and reallocated on every iteration.
	 * this allows some details of line emission to change for last iteration,
	 * but may not really be necessary */
	if( LineDSv!=NULL )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, "IterStart has freed LindDSv \n" );
		}
		free(LineDSv );
		LineDSv = NULL;
	}

	if( LineSv!=NULL )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, "IterStart has freed LindSv \n" );
		}
		free(LineSv );
		LineSv = NULL;
	}

	/* only allocate space if this sum is positive - will be zero if no grains present */
	if( LineSave.ndsum > 0 )
	{
		if( (LineDSv=(LinDstSv *)MALLOC((unsigned)LineSave.ndsum*sizeof( LinDstSv ) ) ) == NULL )
		{
			fprintf( ioQQQ, "IterStart could not allocate space for LineSave.ndsum array, needed %li\n",
				LineSave.ndsum );
			puts( "[Stop in IterStart]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* this is the large main line array */
	if( (LineSv=(LinSv*)MALLOC((unsigned)LineSave.nsum*sizeof( LinSv ) ) ) == NULL )
	{
		fprintf( ioQQQ, "IterStart could not allocate space for LineSave.nsum array, needed %li\n",
			LineSave.nsum );
		puts( "[Stop in IterStart]" );
		cdEXIT(EXIT_FAILURE);
	}

	for( i=0; i < LineSave.nsum; i++ )
	{
		/* this is array of line energies, only defined for some lines,
		 * this is sentinel saying whether line should be added into total continuum */
		LineSv[i].xLineEnergy = 0.;
	}

	LineSave.nsum = 0;

	/* second call to set up line labels */
	LineSave.ipass = 0;
	lines();
	/* in the future calls to lines will result in integrations */
	LineSave.ipass = 1;

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "%7ld lines printed in main line array\n", 
		  LineSave.nsum );
	}

	/* now zero out some variables set by last call to LINES */
	hydro.cintot = 0.;
	thermal.totcol = 0.;
	rfield.comtot = 0.;
	thermal.FreeFreeTotHeat = 0.;
	thermal.power = 0.;
	thermal.HeatLineMax = 0.;
	thermal.GBarMax = 0.;
	hyperfine.cooling_max = 0.;

	gv.TotalDustHeat = 0.;
	gv.dphmax = 0.;
	hmi.h2pmax = 0.;
	gv.dclmax = 0.;
	gv.GrnElecDonateMax = 0.;
	gv.GrnElecHoldMax = 0.;

	co.codfrc = 0.;
	co.codtot = 0.;

	hmi.h2dfrc = 0.;
	hmi.h2line_cool_frac = 0.;
	hmi.h2dtot = 0.;
	timesc.sound = 0.;

	if( LineSave.lgNormFound == FALSE && LineSave.lgNormSet == TRUE )
	{
		LineSave.ipNormWavL = cdLine( LineSave.chNormLab , LineSave.WavLNorm , &fhe1nx , &ratio );
		if( LineSave.ipNormWavL < 0 )
		{
			/* did not find the line if return is negative */
			fprintf( ioQQQ, "There was a problem resetting the normalisation line.\n");
			fprintf( ioQQQ, 
			"IterStart could not find a line with a wavelength of %f and a label of %s\n", 
			  LineSave.WavLNorm, LineSave.chNormLab );
			fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
			fprintf( ioQQQ, "Sorry.\n");
			puts( "[Stop in IterStart]" );
			cdEXIT(EXIT_FAILURE);
		}

		LineSave.ipEmerNormWavL = cdDLine( LineSave.chNormLab , LineSave.WavLNorm , &fhe1nx , &ratio );
		/* possible that valid line specied but it did not happen to be one of the
		 * predicted emergent lines - use hbeta */
		if( LineSave.ipEmerNormWavL < 0 )
		{
			LineSave.ipEmerNormWavL = 0;
		}
	}
#	if 0
	/* find which line to normalise array to (default is H BETA) */
	if( LineSave.lgNormFound == FALSE && LineSave.lgNormSet == TRUE )
	{
		LineSave.lgNormFound = TRUE;

		lgHit = FALSE;

		i = 1;

		/* this must be reset or it will crash */
		j = -LONG_MAX;

		while( !lgHit && i < LineSave.nsum )
		{
			/* if( fabs(LineSv[i].wavelength - LineSave.WavLNorm)/SDIV(LineSave.WavLNorm) < 
				LineSave.errorwave ) */
			if( fabs(LineSv[i].wavelength - LineSave.WavLNorm) < LineSave.errorwave )
			{
				/* no label was set, so any will do */
				if( strcmp( LineSave.chNormLab,"    ") == 0 )
				{
					j = i;
					lgHit = TRUE;
				}
				else
				{
					/* check whether labels match  */
					cap4(chCAPS , (char*)LineSv[i].chALab);

					if( strcmp( LineSave.chNormLab,chCAPS) == 0 )
					{
						/* when chNormLab set, only count as hit when label matches too */
						j = i;
						lgHit = TRUE;
					}
				}
			}
			++i;
		}

		if( !lgHit )
		{
			/* did not find the line */
			fprintf( ioQQQ, "There was a problem resetting the normalisation line.\n");
			fprintf( ioQQQ, 
			"IterStart could not find a line with a wavelength of %f and a label of %s\n", 
			  LineSave.WavLNorm, LineSave.chNormLab );
			fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
			fprintf( ioQQQ, "Sorry.\n");
			puts( "[Stop in IterStart]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* following confirms that index within bounds, mostly for lint */
		ASSERT( j>0 );
		ASSERT( j<LineSave.nsum);

		LineSave.ipNormWavL = j;

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				" Reormalize line spectrum to line number %li wavelength is %f scale=%5.2f label is %s\n", 
			  LineSave.ipNormWavL, 
			  LineSv[j-1].wavelength, 
			  LineSave.ScaleNormLine, 
			  LineSv[j-1].chALab  );
		}
	}
#	endif

	/* if STPINT not 1E37 then stop integration when line LineStopWl
	 * reaches intensity STPINT relative to H BETA */
	/* nstpl is number of stop line commands, 0 if none entered */
	for( ii=0; ii < StopCalc.nstpl; ii++ )
	{
		if( iteration <= 1 )
		{
			i = 0;
			lgHit = FALSE;
			/* >>chng 96 dec 12, got rid of if, change to do while */
			/* search through lines to find wavelength that matches */
			while( !lgHit && i < LineSave.nsum )
			{
				/* >>chng 01 aug 15, add test for label */
				if( fabs(LineSv[i].wavelength - StopCalc.StopLineWl[ii])/
					SDIV(StopCalc.StopLineWl[ii])<StopCalc.ErrorLineWl[ii] &&
					strcmp(StopCalc.chStopLabel[ii] , LineSv[i].chALab) == 0 )
				{
					/* set to array index of proper line */
					StopCalc.ipStopLin1[ii] = i;
					lgHit = TRUE;
				}
				++i;
			}
			/* back off since incremented too many times above */
			--i;

			if( !lgHit )
			{
				fprintf( ioQQQ, " IterStart 1 couldnt find line%5ld\n", 
				  StopCalc.ipStopLin1[ii] );
				puts( "[Stop in IterStart]" );
				cdEXIT(EXIT_FAILURE);
			}

			if( trace.lgTrace )
			{
				fprintf( ioQQQ, 
					" stop line 1 is number %5ld wavelength is %f label is %4.4s\n", 
				  StopCalc.ipStopLin1[ii], 
				  LineSv[StopCalc.ipStopLin1[ii]].wavelength, 
				  LineSv[StopCalc.ipStopLin1[ii]].chALab );
			}

			/* line index was set to zero in zero.c, if non-zero here then we want
			 * the seconc line - if zero the second line not entered */
			if( StopCalc.StopLineWl2[ii] != 0 )
			{
				i = 0;
				lgHit = FALSE;
				while( !lgHit && i < LineSave.nsum )
				{
					/* >>chng 96 dec 12, more to do while logic
					 * do i=1,nsum */
					if( fabs(LineSv[i].wavelength - StopCalc.StopLineWl2[ii])/SDIV(StopCalc.StopLineWl2[ii])<StopCalc.ErrorLineWl2[ii] &&
						strcmp(StopCalc.chStopLabel2[ii] , LineSv[i].chALab) == 0 )
					{
						StopCalc.ipStopLin2[ii] = i;
						lgHit = TRUE;
					}
					++i;
				}
				/* back off since incremented too many times above */
				--i;

				if( !lgHit )
				{
					fprintf( ioQQQ, " IterStart 2 couldnt find line with wavelength %ld\n", 
					  StopCalc.ipStopLin2[ii] );
					puts( "[Stop in IterStart]" );
					cdEXIT(EXIT_FAILURE);
				}

				if( trace.lgTrace )
				{
					fprintf( ioQQQ, 
						" stop line 2 is number %5ld wavelength is %f label is %4.4s\n", 
					  StopCalc.ipStopLin2[ii], 
					  LineSv[StopCalc.ipStopLin2[ii]].wavelength, 
					  LineSv[StopCalc.ipStopLin2[ii]].chALab  );
				}
			}
		}
	}

	/* option to only print last iteration */
	if( prt.lgPrtLastIt )
	{
		if( iteration == 1 )
		{
			called.lgTalk = FALSE;
		}

		/* initial condition of TALK may be off if AMOEBA used
		 * sec part is for print last command 
		 * lgTalkForcedOff is set true when cdTalk is called with false
		 * to turn off printing */
		if( iterations.lgLastIt && !prt.lgPrtStart && !called.lgTalkForcedOff )
		{
			called.lgTalk = intalk.lgInTalk;
		}
	}

	if( trace.lgTrace && trace.lgOptcBug )
	{
		if( opac.lgTauOutOn )
		{
			fprintf( ioQQQ, " Outer optical depths defined.\n" );
		}
		else
		{
			fprintf( ioQQQ, " Outer optical depths NOT defined.\n" );
		}

		/* heavy element lines */
		for( i=0; i < nLevel1; i++ )
		{
			fprintf( ioQQQ, "%6f:%8.1e%8.1e%8.1e", 
				TauLines[i].WLAng, 
			  TauLines[i].TauIn, 
			  TauLines[i].TauTot, 
			  TauLines[i].Pesc );
		}

		fprintf( ioQQQ, "\n" );
	}

	if( opac.lgCaseB )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " IterStart does not change mid-zone optical depth since CASE B\n" );
		}
	}

	/* check if induced recombination can be ignored */
	hydro.FracInd = 0.;
	hydro.fbul = 0.;

	/* remember some things about effects of induced rec on H onlye
	 * don't do ground state since SPHERE turns it off */
	for( i=ipH2p; i < (iso.numLevels[ipH_LIKE][ipHYDROGEN] - 1); i++ )
	{
		ratio = iso.RecomInducRate[ipH_LIKE][ipHYDROGEN][i]*iso.PopLTE[ipH_LIKE][ipHYDROGEN][i]/(iso.RecomInducRate[ipH_LIKE][ipHYDROGEN][i]*iso.PopLTE[ipH_LIKE][ipHYDROGEN][i] + 
			iso.RadRecomb[ipH_LIKE][ipHYDROGEN][i][ipRecRad]*
			iso.RadRecomb[ipH_LIKE][ipHYDROGEN][i][ipRecNetEsc]);
		if( ratio > hydro.FracInd )
		{
			hydro.FracInd = (float)ratio;
			hydro.ndclev = i;
		}

		ratio = EmisLines[ipH_LIKE][ipHYDROGEN][i+1][i].pump/
			(EmisLines[ipH_LIKE][ipHYDROGEN][i+1][i].pump + 
			EmisLines[ipH_LIKE][ipHYDROGEN][i+1][i].Aul);

		if( ratio > hydro.fbul )
		{
			hydro.fbul = (float)ratio;
			hydro.nbul = i;
		}
	}

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

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

/*===============================================================*/
/*reset reset many variables at the start of a new iteration */
void IterRestart(void)
{
	long int i, 
	  ion, 
	  ipISO ,
	  n, 
	  nelem;
	double SumOTS ;

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

	dynamics.dDensityDT = dndtSave;

	/* this is case where temperature floor has been set, if it was hit
	 * then we did a constant temperature calculation, and must go back to 
	 * a thermal solution */
	if( StopCalc.TeFloor > 0. )
	{
		thermal.lgTSetOn = FALSE;
		thermal.ConstTemp = 0.;
	}

	conv.lgAbort = FALSE;

	/* reset some parameters needed for magnetic field */
	Magnetic_reinit();

	hydro.xLaMase = 0.;
	opac.stimax[0] = 0.;
	opac.stimax[1] = 0.;

	H2_Reset();

	ca.oldcla = 0.;

	/* reset molecular fractions */
	for ( i=0; i<N_H_MOLEC; i++) 
	{
		hmi.Hmolec[i] = hmsav[i];
	}

	/* zero out the column densities */
	for( i=0; i < NCOLD; i++ )
	{
		colden.colden[i] = 0.;
	}
	colden.He123S = 0.;

	for( i=0; i < NUM_COMOLE_CALC; i++ )
	{
		/* column densities */
		co.hevcol[i] = 0.;
		/* largest fraction of atoms in molecules */
		co.xMoleFracMax[i] = 0.;
	}

	/* compute timestep if advection is enabled */
	if( dynamics.lgAdvection )
		DynaEndIter();

	hmi.BiggestH2 = -1.f;
	hmi.h2plus_heat = h2plus_heat_save;
	hmi.HeatH2Dish_used = (float)HeatH2Dish_used_save;
	hmi.HeatH2Dexc_used = (float)HeatH2Dexc_used_save;
	hmi.HeatH2Dish_TH85 = (float)HeatH2Dish_TH85_save;
	hmi.HeatH2Dexc_TH85 = (float)HeatH2Dexc_TH85_save;
	hmi.HeatH2Dish_BigH2 = (float)HeatH2Dish_BigH2_save;
	hmi.HeatH2Dexc_BigH2 = (float)HeatH2Dexc_BigH2_save;
	hmi.deriv_HeatH2Dexc_BigH2 = (float)deriv_HeatH2Dexc_BigH2_save;
	hmi.hmihet = hmihet_save;
	hmi.hmitot = hmitot_save;
	hmi.H2_Solomon_dissoc_rate_used = H2_Solomon_dissoc_rate_used_save;
	hmi.H2_Solomon_dissoc_rate_used_H2g = H2_Solomon_dissoc_rate_used_H2g_save;
	hmi.H2_Solomon_dissoc_rate_used_H2s = H2_Solomon_dissoc_rate_used_H2s_save;
	hmi.H2_H2g_to_H2s_rate_used = H2_H2g_to_H2s_rate_used_save;
	hmi.H2_photodissoc_used = H2_photodissoc_used_save;

	hmi.H2_photodissoc_BigH2 = H2_photodissoc_BigH2_save;
	hmi.H2_H2g_to_H2s_rate_BigH2 = H2_H2g_to_H2s_rate_BigH2_save;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2s = H2_Solomon_dissoc_rate_BigH2_H2s_save;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2g = H2_Solomon_dissoc_rate_BigH2_H2g_save;

	rfield.ipEnergyBremsThin = 0;
	rfield.EnergyBremsThin = 0.;
	rfield.lgUSphON = FALSE;

	radius.glbdst = 0.;
	/* zone thickness, and next zone thickness */
	radius.drad = drSave;
	radius.drNext = drNextSave;

	/* was min dr ever used? */
	radius.lgDrMinUsed = FALSE;

	negdrg.lgNegGrnDrg = FALSE;

	/* simple estimate of H-beta intensity */
	continuum.cn4861 = continuum.sv4861;
	continuum.cn1216 = continuum.sv1216;

	/* total luminosity */
	continuum.totlsv = continuum.TotalLumin;

	/* set debugger on now if NZONE desired is 0 */
	if( (trace.nznbug == 0 && iteration >= trace.npsbug) && trace.lgTrOvrd )
	{
		if( trace.lgTrConvg==0 )
		{
			trace.lgTrace = TRUE;
		}
		else
			/* trace convergence entered = but with negative flag = make positive,
			 * abs and not mult by -1 since may trigger more than one time */
			trace.lgTrConvg = abs( trace.lgTrConvg );

		fprintf( ioQQQ, " IterRestart called.\n" );
	}
	else
	{
		trace.lgTrace = FALSE;
	}

	/* zero secondary suptrathermals variables for first ionization attem */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		for( ion=0; ion<nelem+1; ++ion )
		{
			secondaries.csupra[nelem][ion] = supsav[nelem][ion];
		}
	}
	secondaries.x12tot = secondaries.x12sav;
	secondaries.heatef = secondaries.hetsav;
	secondaries.efionz = secondaries.savefi;
	for( nelem=0; nelem<LIMELM; ++nelem)
	{
		for( ion=0; ion<nelem+1 ; ++ion )
		{
			ionbal.CompRecoilHeatRate[nelem][ion] = ionbal.CompRecoilHeatRateSave[nelem][ion];
			ionbal.CompRecoilIonRate[nelem][ion] = ionbal.CompRecoilIonRateSave[nelem][ion];
		}
	}

	wind.lgVelPos = TRUE;
	wind.AccelMax = 0.;
	wind.AccelAver = 0.;
	wind.acldr = 0.;
	wind.windv = wind.windv0;

	thermal.nUnstable = 0;
	thermal.lgUnstable = FALSE;

	pressure.pbeta = 0.;
	pressure.RadBetaMax = 0.;
	pressure.lgPradCap = FALSE;
	pressure.lgPradDen = FALSE;
	/* this flag will say we hit the sonic point */
	pressure.lgSonicPoint = FALSE;
	pressure.lgStrongDLimbo = FALSE;
	pressure.PresInteg = 0.;
	pressure.pinzon = 0.;

	dense.eden = edsav;
	dense.EdenHCorr = edsav;
	dense.EdenTrue = edsav;

	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		/* reset molecular densities */
		dense.gas_phase[nelem] = gas_phase_save[nelem];
		dense.xMolecules[nelem] = xMolFracsSave[nelem];
		dense.IonLow[nelem] = IonLowSave[nelem];
		dense.IonHigh[nelem] = IonHighSave[nelem];
	}
	/* fix number of stages of ionization */


	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				for( n=ipH1s; n < iso.numLevels[ipISO][nelem]; n++ )
				{
					iso.ConOpacRatio[ipISO][nelem][n] = HOpacRatSav[ipISO][nelem][n];
					iso.Pop2Ion[ipISO][nelem][n] = hnsav[ipISO][nelem][n];
				}
			}
		}
	}

	if( trace.lgTrace && trace.lgNeBug )
	{
		fprintf( ioQQQ, "     EDEN set to%12.4e by IterRestart.\n", 
		  dense.eden );
	}

#	if 0
	for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			/* reset total gas phase abundance to the original set */
			dense.gas_phase[nelem] = abund.solar[nelem]*dense.gas_phase[ipHYDROGEN];
		}
		else
		{
			/* >>chng 04 apr 20, set to zero if element is off */
			dense.gas_phase[nelem] = 0.;
		}
	}
#	endif
	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		for( ion=0; ion < (nelem + 2); ion++ )
		{
			dense.xIonDense[nelem][ion] = xIonFsave[nelem][ion];
		}
		for( i=0; i < LIMELM; i++ )
		{
			thermal.heating[nelem][i] = HeatSave[nelem][i];
		}
	}

	GrainRestartIter();

	/* continuum was saved in FluxSave */
	for( i=0; i < rfield.nupper; i++ )
	{
		rfield.flux[i] = rfield.FluxSave[i];
		/* >>chng 03 nov 08, recover SummedDif */
		rfield.SummedDif[i] = rfield.SummedDifSave[i];
		/* >>chng 03 nov 08, update SummedCon here since flux changed */
		rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];

		rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i];
		rfield.otscon[i] = rfield.otssav[i][0];
		rfield.otslin[i] = rfield.otssav[i][1];
		/* >>chng 01 apr 10, add these two resets */
		/* >>chng 03 may 20, had set to otscon and otslin, which is whatever
		 * was left from previous iteration - change to zero */
		rfield.otsconNew[i] = 0.;
		rfield.otslinNew[i] = 0.;
		rfield.ConInterOut[i] = 0.;
		rfield.OccNumbBremsCont[i] = 0.;
		rfield.OccNumbDiffCont[i] = 0.;
		rfield.outlin[i] = 0.;
		rfield.outlin_noplot[i] = 0.;
		rfield.ConOTS_local_OTS_rate[i] = 0.;
		rfield.ConOTS_local_photons[i] = 0.;

		opac.opacity_abs[i] = opac.opacity_abs_savzon1[i];
		opac.OldOpacSave[i] = opac.opacity_abs_savzon1[i];
		opac.opacity_sct[i] = opac.opacity_sct_savzon1[i];
		opac.albedo[i] = 
			opac.opacity_sct[i]/MAX2(1e-30,opac.opacity_sct[i] + opac.opacity_abs[i]);
		opac.tmn[i] = 1.;
		/* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
		 * for heating in very high T models in map.in.  zone 1 and 2 were 20% different,
		 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
		opac.ExpmTau[i] = 1.;
		opac.e2TauAbs[i] = 1.;*/
		/* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
		 * for heating in very high T models in map.in.  zone 1 and 2 were 20% different,
		 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
		 * these were not here at all*/
		/* attenuation of flux by optical depths IN THIS ZONE 
		 * AngleIllum is 1/COS(theta), is usually 1, reset with illuminate command,
		 * option for illumination of slab at an angle */
		opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad/2.*geometry.AngleIllum);

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

		/* e2(tau) in inward direction, up to illuminated face */
		opac.e2TauAbs[i] = (float)e2(opac.TauAbsFace[i],opac.ExpmTau[i]);
		rfield.reflin[i] = 0.;
		rfield.ConEmitReflec[i] = 0.;
		rfield.ConEmitOut[i] = 0.;
		rfield.ConRefIncid[i] = 0.;
	}

	/* update continuum */
	RT_OTS_Update(&SumOTS  , 0.);

	thermal.FreeFreeTotHeat = 0.;

	if( called.lgTalk )
	{
		fprintf( ioQQQ, "\f\n          Start Iteration Number%2ld   %75.75s\n\n\n", 
		  iteration, input.chTitle );
	}

	/* reset variable to do with FeII atom, in FeIILevelPops */
	FeIIReset();

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

/* do some work with ending the iteration */
void IterEnd(void)
{
	long i;
	double fac,
		tau;

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

	/* give indication of geometry */
	fac = radius.depth/radius.rinner;
	if( fac < 0.1 )
	{
		geometry.lgGeoPP = TRUE;
	}
	else
	{
		geometry.lgGeoPP = FALSE;
	}

	/* all continue were attenuated in last zone in radius_increment to represent the intensity
	 * in the middle of the next zone - this was too much since we now want
	 * intensity emergent from outer edge of last zone */
	for( i=0; i < rfield.nflux; i++ )
	{
		{
			/*@-redef@*/
			enum{DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC)
			{
				fprintf(ioQQQ,"i=%li opac %.2e \n", i, 
					(double)opac.opacity_abs[i]*radius.drad_x_fillfac/2.*(double)geometry.AngleIllum );
			}
		}
		tau = opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.AngleIllum;
		ASSERT( tau > 0. );
		fac = sexp( tau );

		/* >>chng 02 dec 14, add first test to see whether product of ratio is within
		 * range of a float - ConInterOut can be finite and fac just above a small float,
		 * so ratio exceeds largest size of a float */
		/*if( fac > SMALLFLOAT )*/
		if( (float)(fac/SDIV(rfield.ConInterOut[i]))>SMALLFLOAT && fac > SMALLFLOAT )
		{
			rfield.ConInterOut[i] /= (float)fac;
			rfield.outlin[i] /= (float)fac;
			rfield.outlin_noplot[i] /= (float)fac;
		}
	}

	/* remember thickness of previous iteration */
	radius.router[iteration-1] = radius.depth;

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