/* 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 */
/*ConvBase main routine to drive ionization solution for all species, find total opacity
 * called by ConvIoniz */
/*lgConverg check whether ionization of element nelem has converged */
#include "cddefines.h"
#include "dynamics.h"
#include "hydrogenic.h"
#include "trace.h"
#include "elementnames.h"
#include "phycon.h"
#include "secondaries.h"
#include "grainvar.h"
#include "highen.h"
#include "dense.h"
#include "hmi.h"
#include "rfield.h"
#include "pressure.h"
#include "taulines.h"
#include "helike.h" 
#include "rt.h"
#include "grains.h"
#include "atmdat.h"
#include "ionbal.h"
#include "opacity.h"
#include "cooling.h"
#include "thermal.h"
#include "mole.h"
#include "iso.h"
#include "conv.h"

/*lgIonizConverg check whether ionization of element nelem has converged, called by ionize,
 * returns TRUE if element is converged, FALSE if not */
static int lgIonizConverg(
	/* atomic number on the C scale, 0 for H, 25 for Fe */
	long int nelem ,
	/* this is allowed error as a fractional change.  Most are 0.15 */
	double delta ,
	/* option to print abundances */
	int lgPrint )
{
	int lgConverg_v;
	long int ion;
	double Abund, 
	  bigchange ,
	  change ,
	  one ;
	double abundold=0. , abundnew=0.;

	/* this is fractions [ion stage][nelem], ion stage = 0 for atom, nelem=0 for H*/
	static float OldFracs[LIMELM+1][LIMELM+1];

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

	if( !dense.lgElmtOn[nelem] )
		return TRUE;

	/* args are atomic number, ionization stage
	 * OldFracs[nelem][0] is abundance of nelem (cm^-3) */

	/* this function returns true if ionization of element
	 * with atomic number nelem has not changed by more than delta,*/

	/* check whether abundances exist yet, only do this for after first zone */
	/*if( nzone > 0 )*/
	/* >>chng 03 sep 02, check on changed ionization after first call through,
	 * to insure converged constant eden / temperature models */
	if( conv.nPres2Ioniz )
	{
		/* >>chng 04 aug 31, this had been static, caused last hits on large changes
		 * in ionization */
		bigchange = 0.;
		change = 0.;
		Abund = dense.gas_phase[nelem];

		/* loop over all ionization stages, loop over all ions, not just active ones,
		 * since this also sets old values, and so will zero out non-existant ones */
		for( ion=0; ion <= (nelem+1); ++ion )
		{
			/*lint -e727 OlsdFracs not initialized */
			if( OldFracs[nelem][ion]/Abund > 1e-4 && 
				dense.xIonDense[nelem][ion]/Abund > 1e-4 )
			/*lint +e727 OlsdFracs not initialized */
			{
				/* check change in old to new value */
				one = fabs(dense.xIonDense[nelem][ion]-OldFracs[nelem][ion])/
					OldFracs[nelem][ion];
				change = MAX2(change, one );
				/* remember abundances for largest change */
				if( change>bigchange )
				{
					bigchange = change;
					abundold = OldFracs[nelem][ion]/Abund;
					abundnew = dense.xIonDense[nelem][ion]/Abund;
				}
			}
			/* now set new value */
			OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
		}

		if( change < delta )
		{
			lgConverg_v = TRUE;
		}
		else
		{
			lgConverg_v = FALSE;
			conv.BadConvIoniz[0] = abundold;
			conv.BadConvIoniz[1] = abundnew;
			ASSERT( abundold>0. && abundnew>0. );
		}
	}
	else
	{
		for( ion=0; ion <= (nelem+1); ++ion )
		{
			OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
		}

		lgConverg_v = TRUE;
	}

	/* option to print abundances */
	if( lgPrint )
	{
		fprintf(ioQQQ," element %li converged? %c ",nelem, TorF(lgConverg_v));
		for( ion=0; ion<(nelem+1); ++ion )
		{
			fprintf(ioQQQ,"\t%.2e", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]);
		}
		fprintf(ioQQQ,"\n");
	}

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

/*ConvBase main routine to drive ionization solution for all species, find total opacity
 * called by ConvIoniz
 * return 0 if ok, 1 if abort */
int ConvBase( 
	/* this is zero the first time ConvBase is called by convIoniz,
	 * counts number of call thereafter */
	 long loopi )
{
	double O_HIonRate_New , O_HIonRate_Old ;
	double HeatOld,
		EdenTrue_old,
		EdenFromMolecOld,
		EdenFromGrainsOld,
		HeatingOld ,
		CoolingOld ;
	float comole_save[NUM_COMOLE_CALC];
	float hmole_save[N_H_MOLEC];
	static double SecondOld;
	static long int nzoneOTS=-1;
#	define LOOP_ION_LIMIT	10
	long int nelem, 
		ion,
		loop_ion,
		i,
		ipISO;
	static double SumOTS=0. , OldSumOTS[2]={0.,0.} ;
	double save_iso_grnd[NISO][LIMELM];

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

	/* this is set to phycon.te in tfidle, is used to insure that all temp
	 * vars are properly updated when conv_ionizeopacitydo is called 
	 * NB must be same type as phycon.te */
	/*lint -e777 test float for equal */
	ASSERT( phycon.te == thermal.te_update );
	/*lint +e777 test float for equal */

	/* this allows zone number to be printed with slight increment as zone converged
	 * conv.nPres2Ioniz is incremented at the bottom of this routine */
	fnzone = (double)nzone + (double)conv.nPres2Ioniz/100.;

	/* reevaluate pressure */
	/* this sets values of pressure.PresTotlCurr, also calls tfidle,
	 * and sets the total energy content of gas, which may be important for acvection */
	PresTotCurrent();

	/* >>chng 04 sep 15, move EdenTrue_old up here, and will redo at bottom
	 * to find change 
	 * find and save current true electron density - if this changes by more than the
	 * tolerance then ionizaiton soln is not converged */
	/* >>chng 04 jul 27, move eden_sum here from after this loop, so that change in EdenTrue
	 * can be monitored */
	/* update EdenTrue, eden itself is actually changed in ConvEdenIoniz */
	/* must not call eden_sum on very first time since for classic PDR total
	 * ionization may still be zero on first call */
	if( conv.nTotalIoniz )
	{
		if( eden_sum() )
		{
			/* non-zero return indicates abort condition */
			++conv.nTotalIoniz;
#			ifdef DEBUG_FUN
			fputs( " <->ConvBase()\n", debug_fp );
#			endif
			return 1;
		}
	}

	/* the following two were evaluated in eden_sum 
	 * will confirm that these are converged */
	EdenTrue_old = dense.EdenTrue;
	EdenFromMolecOld = co.comole_eden;
	EdenFromGrainsOld = gv.TotalEden;
	HeatingOld = thermal.htot;
	CoolingOld = thermal.ctot;

	/* remember current ground state populatoin - will check if converged */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem<LIMELM;++nelem )
		{
			if( dense.lgElmtOn[nelem] )
			{
				/* save the ground state population */
				save_iso_grnd[ipISO][nelem] = iso.Pop2Ion[ipISO][nelem][0];
			}
		}
	}

	for( i=0; i < NUM_COMOLE_CALC; i++ )
	{
		comole_save[i] = co.hevmol[i];
	}

	for( i=0; i < N_H_MOLEC; i++ )
	{
		hmole_save[i] = hmi.Hmolec[i];
	}

	if( loopi==0 )
	{
		/* these will be used to look for oscillating ots rates */
		OldSumOTS[0] = 0.;
		OldSumOTS[1] = 0.;
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			"   ConvBase called. %.2f Te:%.3e  HI:%.3e  HII:%.3e  H2:%.3e  Ne:%.3e  htot:%.3e  CSUP:%.2e Conv?%c\n", 
		  fnzone,
		  phycon.te, 
		  dense.xIonDense[ipHYDROGEN][0], 
		  dense.xIonDense[ipHYDROGEN][1], 
		  hmi.H2_total,
		  dense.eden, 
		  thermal.htot, 
		  secondaries.csupra[ipHYDROGEN][0] ,
		  TorF(conv.lgConvIoniz) );
	}
	/* want this flag to be TRUE when we exit, varioius problems will set falst */
	conv.lgConvIoniz = TRUE;

	/* this routine is in heatsum.c, and zeroes out the heating array */
	HeatZero();

	/* if this is very first call, say not converged, since no meaningful way to
	 * check on changes in quantities.  this counter is false only on very first
	 * call, when equal to zero */
	if( !conv.nTotalIoniz )
	{
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "first call" );
	}

	/* this will be flag to check whether any ionization stages
	 * were trimed down */
	conv.lgIonStageTrimed = FALSE;

	/* must redo photoionizaiton rates for all species on second and later tries */
	/* always reevaluate the rates when . . . */
	/* this flag says not to redo opac and photo rates, and following test will set
	 * true if one of several tests not done*/
	opac.lgRedoStatic = FALSE;
	if( 
		/* opac.lgOpacStatic (usually TRUE), is set FALSE with no static opacity command */
		!opac.lgOpacStatic || 
		/* we are in search mode */
		conv.lgSearch || 
		/* this is the first call to this zone */
		conv.nPres2Ioniz == 0 )
	{
		/* we need to redo ALL photoionization rates */
		opac.lgRedoStatic = TRUE;
	}

	/* reevaluate DR rates if temperature has changed 
	 * >>chng 05 dec 19, add this rate */
	{
		static double TeUsed = -1;
		if( fabs(phycon.te/TeUsed - 1. ) > 0.001 )
		{
			/*>>chng 06 jan 31 save mean rec coef */
			long int nsumrec[LIMELM];
			for( ion=0; ion<LIMELM; ++ion )
			{
				ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 0.;
				nsumrec[ion] = 0;
			}
			TeUsed = phycon.te;
			/* NB - for following loop important to go over all elements and all
			 * ions, not just active ones, since mean DR is used as the guess for
			 * DR rates for unknown species.  */
			for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
			{
				for( ion=0; ion<nelem+1; ++ion )
				{
					long int n_bnd_elec_before_recom ,
						     n_bnd_elec_after_recom;

					n_bnd_elec_before_recom = nelem-ion;
					n_bnd_elec_after_recom = nelem-ion+1;

					/* Badnell dielectronic recombination rate coefficients */
					if( (ionbal.DR_Badnell_rate_coef[nelem][ion] = Badnell_DR_rate_eval(
						/* atomic number on physics scale - He - 2 */
						nelem+1, 
						/* number of core electrons before capture of free electron */
						n_bnd_elec_before_recom )) >= 0. )
					{
						ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = TRUE;
						/* keep track of mean DR rate for this ion */
						++nsumrec[ion];
						ionbal.DR_Badnell_rate_coef_mean_ion[ion] += ionbal.DR_Badnell_rate_coef[nelem][ion];
					}
					else
					{
						ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = FALSE;
						ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.;
					}
					/* Badnell radiative recombination rate coefficients */
					if( (ionbal.RR_Badnell_rate_coef[nelem][ion] = Badnell_RR_rate_eval(
						/* atomic number on physics scale - He - 2 */
						nelem+1, 
						/* number of core electrons before capture of free electron */
						n_bnd_elec_before_recom )) >= 0. )
					{
						ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = TRUE;
					}
					else
					{
						ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = FALSE;
						ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.;
					}
					/* save D. Verner's radiative recombination rate coef 
					 * needed for rec cooling, cm3 s-1 */
					ionbal.RR_Verner_rate_coef[nelem][ion] = 
						atmdat_rad_rec( 
						/* number of bound electrons after recom */
						nelem+1 , 
						/* number of protons on physics scale */
						n_bnd_elec_after_recom , 
						phycon.te );
				}
			}
			/* now get mean of the DR rates - may be used for ions with no DR rates */
			for( ion=0; ion<LIMELM; ++ion )
			{
				if( nsumrec[ion] > 0 )
					ionbal.DR_Badnell_rate_coef_mean_ion[ion] /= nsumrec[ion];
					/*fprintf(ioQQQ,"DEBUG %li %.2e\n", ion , ionbal.DR_Badnell_rate_coef_mean_ion[ion] );*/
			}
			/*exit(1);*/

			/* this set true with PRINT on any of the Badnell set recom commands */
			if( ionbal.lgRecom_Badnell_print )
			{
				fprintf(ioQQQ,"DEBUG Badnell recom RR, then DR, T=%.3e\n", phycon.te );
				for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
				{
					fprintf(ioQQQ,"nelem=%li, RR then DR\n", nelem+1 );
					for( ion=0; ion<nelem+1; ++ion )
					{
						if( ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] )
							fprintf(ioQQQ," %.2e", ionbal.RR_Badnell_rate_coef[nelem][ion] );
					}
					fprintf(ioQQQ,"\n" );
					for( ion=0; ion<nelem+1; ++ion )
					{
						if( ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
							fprintf(ioQQQ," %.2e", ionbal.DR_Badnell_rate_coef[nelem][ion] );
					}
					fprintf(ioQQQ,"\n\n" );
				}
				/* now print mean recom and standard deviation */
				fprintf(ioQQQ,"mean recom \n" );
				for( ion=0; ion<LIMELM; ++ion )
				{
					double stddev;
					stddev = 0.;
					for( nelem=ion; nelem<LIMELM; ++nelem )
					{
						stddev += POW2( 
							ionbal.DR_Badnell_rate_coef[nelem][ion]- 
							ionbal.DR_Badnell_rate_coef_mean_ion[ion] );
					}
					stddev = sqrt( stddev / MAX2( 1 , nsumrec[ion]-1 ) );
					fprintf(ioQQQ," %2li %.2e %.2e %.2e\n", ion , ionbal.DR_Badnell_rate_coef_mean_ion[ion] , stddev ,
						stddev / SDIV(ionbal.DR_Badnell_rate_coef_mean_ion[ion]) );
				}
				cdEXIT( EXIT_SUCCESS );
			}
		}
	}
	/* this adjusts the lowest and highest stages of ionization we will consider,
	 * only safe to call when lgRedoStatic is true since this could lower the 
	 * lowest stage of ionization, which needs all its photo rates */

	/* conv.nTotalIoniz is only 0 (FALSE) on the very first call to ConvBase,
	 * when we do not know what the ionization distribution is, since not yet done */
	if( conv.nTotalIoniz )
	{
		int lgion_trim_called = FALSE;
		/* >>chng 04 jun 04, bring he into ion_trim */
		/* ion_trim only used for Li and heavier, not H or He */
		/* only do this one time per zone since up and down cycle can occur */
		/* >>chng 05 jan 15, increasing temperature above default first conditions, also
		 * no trim during search - this fixed major logical error when sim is
		 * totally mechanically heated to coronal temperatures -
		 * problem id'd by Ronnie Hoogerwerf */
		if( !conv.nPres2Ioniz && !conv.lgSearch /*|| conv.lgSearch*/ )
		{
			lgion_trim_called = TRUE;
			for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
			{
				if( dense.lgElmtOn[nelem] )
				{
					/* ion_trim will set conv.lgIonStageTrimed TRUE is any ion has its
					 * lowest stage of ionization dropped or trimmed */
					ion_trim(nelem);
				}
			}
		}

		/* following block only set of asserts */
#		if !defined(NDEBUG)
		/* check that proper abundances are either positive or zero */
		for( nelem=ipHELIUM; nelem<LIMELM; ++nelem)
		{
			if( dense.lgElmtOn[nelem] )
			{
				for( ion=0; ion<dense.IonLow[nelem]; ++ion )
				{
					ASSERT( dense.xIonDense[nelem][ion] == 0. );
				}
				/*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/
				for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
				{
					/* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to
					* trip over VERY small floats that failed on alphas, but not 386
					* 
					* in case where lower ionization stage was just lowered or
					* trimmed down the abundance
					* was set to SMALLFLOAT so test must be < SMALLFLOAT */
					/* >>chng 02 feb 19, add check for search phase.  During this search
					* models with extreme ionization (all neutral or all ionized) can
					* have extreme but non-zero abundances far from the ionization peak for
					* element with lots of electrons.  These will go away once the model
					* becomes stable */
					/* >>chng 03 dec 01, add check on whether ion trim was called 
					 * conserve.in threw assert when iontrim not called and abund grew small */
					ASSERT( conv.lgSearch || !lgion_trim_called ||
						/* this can happen if all C is in the form of CO */
						(ion==0 && dense.IonHigh[nelem]==0 ) ||
						dense.xIonDense[nelem][ion] >= SMALLFLOAT || 
						dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT );
				}
				for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
				{
					ASSERT( dense.xIonDense[nelem][ion] == 0. );
				}
			}
		}
#		endif
	}

	/* now check if anything trimed down */
	if( conv.lgIonStageTrimed )
	{
		/* something was trimed down, so say that ionization not yet stable */
		/* say that ionization has not converged, secondaries changed by too much */
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "IonTrimmed" );
	}

	/* reevaluate advective terms if turned on */
	if( dynamics.lgAdvection )
		DynaIonize();

	/* >>chng 04 feb 15, add loop over ionization until converged.  Non-convergence
	 * in this case means ionization ladder pop changed, probably because of way
	 * that Auger ejection is treated - loop exits when conditions tested at end are met */
	loop_ion = 0;
	do
	{

		/* >>chng 04 sep 11, moved to before grains and other routines - must no clear this
		 * flag before their cals are done */
		/* set the convergece flag to TRUE, 
		 * if anything changes in ConvBase, it will be set false */
		if( loop_ion )
			conv.lgConvIoniz = TRUE;

		/* >>chng 01 apr 29, move charge transfer evaluation here, had been just before
		 * HeLike, but needs to be here so that same rate coef used for H ion and other recomb */
		/* fill in master charge transfer array, and zero out some arrays that track effects,
		 * O_HIonRateis rate oxygen ionizes hydrogen - will need to converge this */
		ChargTranEval( &O_HIonRate_Old );

		/* update rate coefficients, only temp part - in mole_co_etc.c */
		CO_update_chem_rates();

		/* find grain temperature, charge, and gas heating rate */
		/* >>chng 04 apr 15, moved here from after call to HeLike(), PvH */
		GrainDrive();

		/* do all hydrogenic species, and fully do hydrogen itself, including molecules */
		/*fprintf(ioQQQ,"DEBUG h2\t%.2f\t%.3e\t%.3e", fnzone,hmi.H2_total,hmi.Hmolec[ipMH2g]);*/
		Hydrogenic();
		/*fprintf(ioQQQ,"\t%.3e\n", hmi.Hmolec[ipMH2g]);*/

		/* evaluate Compton heating, bound E compton, cosmic rays */
		highen();

		/* find corrections for three-body rec - collisional ionization */
		atmdat_3body();

		/* deduce dielctronic supression factors */
		atmdat_DielSupres();

		/* >>chng 02 mar 08 move helike to before helium */
		/* do all he-like species */
		HeLike();

		/*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures
		 * which cascaded downstream if we did not abort */
		/* return if one of above solvers has declared abort condition */
		if( lgAbort )
		{

			++conv.nTotalIoniz;
#			ifdef DEBUG_FUN
			fputs( " <->ConvBase()\n", debug_fp );
#			endif

			return 1;
		}

		/* find grain temperature, charge, and gas heating rate */
		/*GrainDrive();*/

		/* inner shell ionization */
		for( nelem=0; nelem< LIMELM; ++nelem )
		{
			for( ion=0; ion<nelem+1; ++ion )
			{
				ionbal.xInnerShellIonize[nelem][ion] = 0.;
			}
		}
		/* inner shell ionization by absorbing lines */
		/* this flag turned off with no uta command */
		if( ionbal.lgInnerShellLine_on )
		{
			for( i=0; i<nUTA; ++i )
			{
				if( UTALines[i].Aul > 0. )
				{
					/* cs is actually the negative of the branching ratio for autoionization, 
					* rateone is inverse lifetime of level against autoionizaiton */
					double rateone = UTALines[i].pump * -UTALines[i].cs;
					ionbal.xInnerShellIonize[UTALines[i].nelem-1][UTALines[i].IonStg-1] += rateone;
				}
			}
		}

		/* helium ionization balance */
		IonHelium();

		/* do the ionization balance */
		/* evaluate current opacities, OpacityAddTotal is called only here during actual calculation */
		/* option to only evaluate opacity one time per zone, 
		 * rfield.lgOpacityReevaluate normally true, set false with no opacity reevaluate command */
		/* >>chng 04 jul 16, move these out of block below, so always done 
		 * these all produce ions that are used in the CO network */
		IonCarbo();
		IonOxyge();
		IonNitro();
		IonSilic();
		IonSulph();
		/* do carbon monoxide after oxygen */
		CO_drive();
		if( !conv.nPres2Ioniz || rfield.lgIonizReevaluate || conv.lgIonStageTrimed || conv.lgSearch )
		{
			IonLithi();
			IonBeryl();
			IonBoron();
			/*IonCarbo();
			IonOxyge();*/
			/* do carbon monoxide after oxygen */
			/*fprintf(ioQQQ,"DEBUG co\t%.2f\t%.3e", fnzone,co.hevmol[ipCO]);*/
			/*CO_drive();*/
			/*fprintf(ioQQQ,"\t%.3e\n", co.hevmol[ipCO]);*/
			IonFluor();
			IonNeon();
			IonSodiu();
			IonMagne();
			IonAlumi();
			IonPhosi();
			IonChlor();
			IonArgon();
			IonPotas();
			IonCalci();
			IonScand();
			IonTitan();
			IonVanad();
			IonChrom();
			IonManga();
			IonIron();
			IonCobal();
			IonNicke();
			IonCoppe();
			IonZinc();
		}

		/* all elements have now had ionization reevaluated - in some cases we may have upset
		 * the ionization that was forced with an "element ionization" command - here we will 
		 * re-impose that set ionization */
		/* >>chng 04 oct 13, add this logic */
		for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
		{
			if( dense.lgSetIoniz[nelem] )
			{
				dense.IonLow[nelem] = 0;
				dense.IonHigh[nelem] = nelem + 1;
				while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
					++dense.IonLow[nelem];
				while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
					--dense.IonHigh[nelem];
			}
		}

		/* redo Oxygen ct ion rate of H to see if it changed */
		ChargTranEval( &O_HIonRate_New );

		/* lgIonizConverg is a function to check whether ionization has converged
		 * check whether ionization changed by more than relative error
		 * given by second number */
		/* >>chng 04 feb 14, loop over all elements rather than just a few */
		for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
		{
			if( !lgIonizConverg(nelem, 0.05  ,FALSE ) )
			{
				conv.lgConvIoniz = FALSE;
				sprintf( conv.chConvIoniz , "%2s ion" , elementnames.chElementSym[nelem] );
			}
		}
		
		if( fabs(EdenTrue_old - dense.EdenTrue)/SDIV(dense.EdenTrue) > conv.EdenErrorAllowed/2. )
		{
			conv.lgConvIoniz = FALSE;
			sprintf( conv.chConvIoniz , "EdTr cng" );
			conv.BadConvIoniz[0] = EdenTrue_old;
			conv.BadConvIoniz[1] = dense.EdenTrue;
		}
		++loop_ion;
	}
	/* loop ig not converged, less than loop limit, and we are reevaluating */
	while( !conv.lgConvIoniz && loop_ion < LOOP_ION_LIMIT && rfield.lgIonizReevaluate);

	/* >>chng 05 oct 29, move CT heating here from heat_sum since this sometimes contributes
	 * cooling rather thann heat and so needs to be sorted out before either heating or cooling 
	 * are derived first find net heating - */
	thermal.char_tran_heat = ChargTranSumHeat();
	/* net is cooling if negative */
	thermal.char_tran_cool = MAX2(0. , -thermal.char_tran_heat );
	thermal.char_tran_heat = MAX2(0. ,  thermal.char_tran_heat );

	/* get total cooling, thermal.ctot = does not occur since passes as pointer.  This can add heat.
	 * it calls coolSum at end to sum up the total cooling */
	CoolEvaluate(&thermal.ctot );

	/* get total heating rate - first save old quantities to check how much it changes */
	HeatOld = thermal.htot;

	/* HeatSum will update ElecFrac, 
	 * secondary electron ionization and excitation efficiencies,
	 * and sum up current secondary rates - remember old values before we enter */
	SecondOld = secondaries.csupra[ipHYDROGEN][0];

	/* update the total heating - it was all set to zero in HeatZero at top of this routine,
	 * occurs before secondaries bit below, since this updates elec fracs */
	HeatSum();

	/* test whether we can just set the rate to the new one, or whether we should
	 * take average and do it again.  secondaries.sec2total was set in hydrogenic, and
	 * is ratio of secondary to total hydrogen destruction rates */
	/* >>chng 02 nov 20, add test on size of csupra - primal had very close to underflow */
	if( (secondaries.csupra[ipHYDROGEN][0]>SMALLFLOAT) && (secondaries.sec2total>0.001) && 
		fabs( 1. - SecondOld/SDIV(secondaries.csupra[ipHYDROGEN][0]) ) > 0.1 && 
		SecondOld > 0. &&  secondaries.csupra[ipHYDROGEN][0] > 0.)
	{
		/* say that ionization has not converged, secondaries changed by too much */
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "SecIonRate" );
		conv.BadConvIoniz[0] = SecondOld;
		conv.BadConvIoniz[1] = secondaries.csupra[ipHYDROGEN][0];
	}

#	if 0
	static float hminus_old=0.;
	/* >>chng 04 apr 15, add this convergence test */
	if( conv.nTotalIoniz )
	{
		if( fabs( hminus_old-hmi.Hmolec[ipMHm])/SDIV(hmi.Hmolec[ipMHm]) >
			conv.EdenErrorAllowed/4. )
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "Big H- chn" );
			conv.BadConvIoniz[0] = hminus_old;
			conv.BadConvIoniz[1] = hmi.Hmolec[ipMHm];
	}
	hminus_old = hmi.Hmolec[ipMHm];
#	endif

	if( HeatOld > 0. && !thermal.lgTSetOn )
	{
		/* check if heating has converged - toler is final match */
		if( fabs(1.-thermal.htot/HeatOld) > conv.HeatCoolRelErrorAllowed*.5 )
		{
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "Big d Heat" );
			conv.BadConvIoniz[0] = HeatOld;
			conv.BadConvIoniz[1] = thermal.htot;
		}
	}

	/* abort flag may have already been set - if so bail */
	if( lgAbort )
	{
#		ifdef DEBUG_FUN
		fputs( " <->ConvBase()\n", debug_fp );
#		endif

		return 1;
	}

	/* evaluate current opacities, OpacityAddTotal is called only here during actual calculation */
	/* option to only evaluate opacity one time per zone, 
	 * rfield.lgOpacityReevaluate normally true, set false with no opacity reevaluate command */
	if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate  || conv.lgIonStageTrimed )
		OpacityAddTotal();

	/* >>chng 02 jun 11, call even first time that this routine is called -
	 * this seemed to help convergence */

	/* do OTS rates for all lines and all continua since
	* we now have ionization balance of all species.  Note that this is not
	* entirely self-consistnet, since dest probs here are not the same as
	* the ones used in the model atoms.  Problems??  if near convergence
	* then should be nearly identical */
	if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || rfield.lgIonizReevaluate ||
		conv.lgIonStageTrimed || conv.lgSearch || nzone!=nzoneOTS )
	{
		RT_OTS();
		nzoneOTS = nzone;

		/* remember old ots rates */
		OldSumOTS[0] = OldSumOTS[1];
		OldSumOTS[1] = SumOTS;
		/*fprintf(ioQQQ," calling RT_OTS zone %.2f SumOTS is %.2e\n",fnzone,SumOTS);*/

		/* now update several components of the continuum, this only happens after
		 * we have gone through entire soln for this zone at least one time. 
		 * there can be wild ots oscillation on first call */
		/* the rel change of 0.2 was chosen by running hizqso - smaller increased
		 * itrzn but larger did not further decrease it. */
		RT_OTS_Update(&SumOTS , 0.20 );
		/*fprintf(ioQQQ,"RT_OTS_Updateee\t%.3f\t%.2e\t%.2e\n", fnzone,SumOTS , OldSumOTS[1] );*/
	}
	else
		SumOTS = 0.;

	/* now check whether the ots rates changed */
	if( SumOTS> 0. )
	{
		/* the ots rate must be converged to the error in the electron density */
		/* how fine should this be converged?? orginally had above, 10%, but take
		 * smaller ratio?? */
		if( fabs(1.-OldSumOTS[1]/SumOTS) > conv.EdenErrorAllowed )
		{
			/* this branch, ionization not converged due to large change in ots rates.
			 * check whether ots rates are oscillating, if loopi > 1 so we have enough info*/
			if( loopi > 1 )
			{
				/* here we have three values, are they changing sign? */
				if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
				{
					/* ots rates are oscillating, if so then use small fraction of new dest rates 
					 * when updating Lyman line dest probs in HydroPesc*/
					conv.lgOscilOTS = TRUE;
				}
			}

			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "OTSRatChng" );
			conv.BadConvIoniz[0] = OldSumOTS[1];
			conv.BadConvIoniz[1] = SumOTS;
		}
		/* produce info on the ots fields if either "trace ots" or 
		 * "trace convergence xxx ots " was entered */
		if( ( trace.lgTrace && trace.lgOTSBug ) ||
			( trace.lgTrConvg && trace.lgOTSBug ) )
		{
			RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
		}
		/*fprintf(ioQQQ,"DEBUG opac\t%.2f\t%.3e\t%.3e\n",fnzone,
			dense.xIonDense[ipNICKEL][0] ,
			dense.xIonDense[ipZINC][0] );*/
		{
			/* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC && (nzone>110)/**/ )
			{
#				if 0
#				include "lines_service.h"
				DumpLine( &EmisLines[ipH_LIKE][ipHYDROGEN][2][0] );
#				endif
				/* last par 'l' for line, 'c' for continua, 'b' for both,
				 * the numbers printed are:
				 * cell i, [i], so 1 less than ipoint
				 * anu[i], 
				 * otslin[i], 
				 * opacity_abs[i],
				 * otslin[i]*opacity_abs[i],
				 * rfield.chLineLabel[i] ,
				 * rfield.line_count[i] */
			}
		}
	}
	{
		/* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone>200)/**/ )
		{
			fprintf(ioQQQ,"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n",
				nzone,
				EmisLines[0][1][15][3].ots,
				TauLines[26].ots,
				opac.opacity_abs[2069]);
		}
	}

	/* option to print OTS continuum with TRACE OTS */
	if( trace.lgTrace && trace.lgOTSBug )
	{
		/* find ots rates here, so we only print fraction of it,
		 * SumOTS is both line and continuum contributing to ots, and is mult by opacity */
		/* number is weakest rate to print */
		RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
	}

	/* now reevaluate only destruction probabilities if call is not first*/
	/* >>chng 02 jun 13, had always been false, now reeval escape probs 
	 * and pumping rates on first call for each zone */
	if( conv.nPres2Ioniz && !conv.lgSearch )
	{
		RT_line_all(FALSE , FALSE );
	}
	else
	{
		RT_line_all(TRUE , FALSE );
	}

	/* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
	/* reevaluate pressure */
	/* this sets values of pressure.PresTotlCurr, also calls tfidle */
	PresTotCurrent();

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			"   ConvBase return. %.2f Te:%.3e  HI:%.3e  HII:%.3e  H2:%.3e  Ne:%.3e  htot:%.3e  CSUP:%.2e Conv?%c\n", 
		  fnzone,
		  phycon.te, 
		  dense.xIonDense[ipHYDROGEN][0], 
		  dense.xIonDense[ipHYDROGEN][1], 
		  hmi.H2_total,
		  dense.eden, 
		  thermal.htot, 
		  secondaries.csupra[ipHYDROGEN][0] ,
		  TorF(conv.lgConvIoniz) );
	}

	/* this checks that all molecules and ions add up to gas phase abundance
	 * but only if we have a converged solution */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		/* check that species add up if element is turned on, and
		 * ionization is not set */
		if( dense.lgElmtOn[nelem] && !dense.lgSetIoniz[nelem])
		{
			/* this sum of over the molecules did not include the atom and first
			 * ion that is calculated in the molecular solver */
			double sum = dense.xMolecules[nelem];
			for( ion=0; ion<nelem+2; ++ion )
			{
				sum += dense.xIonDense[nelem][ion];
			}
			ASSERT( sum > SMALLFLOAT );
			/* check that sum agrees with total gas phase abundance */
			/*TODO	2	this assert is not passed if error made much smaller.  This
			 * error should be related to a check on convergence of the molecular
			 * networks and their sum rules, with a criteria used here and there */
			/*if( fabs(sum-dense.gas_phase[nelem])/sum > 1e-2 )
			{
				fprintf(ioQQQ,"DEBUG non conservation element %li sum %.3e gas phase %.3e rel error %.3f\n",
					nelem, sum , dense.gas_phase[nelem],(sum-dense.gas_phase[nelem])/sum);
			}*/
			/* >>chng 04 jul 23, needed to increase to 2e-2 from 1e-2 to get conserve.in to work,
			 * not clear what caused increase in error */
			if( fabs(sum-dense.gas_phase[nelem])/sum > 2e-2 )
			{
				/*fprintf(ioQQQ,"DEBUG non-conv nelem\t%li\tsum\t%e\ttot gas\t%e\trel err\t%.3e\n",
					nelem,
					sum,
					dense.gas_phase[nelem],
					(sum-dense.gas_phase[nelem])/sum);*/
				conv.lgConvIoniz = FALSE;
				strcpy( conv.chConvIoniz, "CO con4" );
				sprintf( conv.chConvIoniz, "COnelem%2li",nelem );
				conv.BadConvIoniz[0] = sum;
				conv.BadConvIoniz[1] = dense.gas_phase[nelem];
			}
		}
	}

	/* update some counters that keep track of how many times this routine
	 * has been called */

	/* this counts number of times we are called by ConvPresTempEdenIoniz, 
	 * number of calls in this zone so first call is zero
	 * reset to zero each time ConvPresTempEdenIoniz is called */
	++conv.nPres2Ioniz;

	/* this is abort option set with SET PRESIONIZ command,
	 * test on nzone since init can take many iterations
	 * this is seldom used except in special cases */
	if( conv.nPres2Ioniz > conv.limPres2Ioniz && nzone > 0)
	{
		fprintf(ioQQQ,"PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz.\n");
		fprintf(ioQQQ,"Their values are %li and %li \n",conv.nPres2Ioniz , conv.limPres2Ioniz);
		lgAbort = TRUE;

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

		return 1;
	}

	/* various checks on the convergence of the current solution */
	/* >>chng 04 sep 15, add this check if the ionization so converged */
	if( eden_sum() )
	{
		/* non-zero return indicates abort condition */
#		ifdef DEBUG_FUN
		fputs( " <->ConvBase()\n", debug_fp );
#		endif
		return 1;
	}
	/* is electron density converged? */
	/*TODO	0	PvH prefers test against err/10 */
	if( fabs(EdenTrue_old - dense.EdenTrue)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "eden chng" );
		conv.BadConvIoniz[0] = EdenTrue_old;
		conv.BadConvIoniz[1] = dense.EdenTrue;
	}

	/* check on molec elec den */
	if( fabs(EdenFromMolecOld - co.comole_eden)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "edn chnCO" );
		conv.BadConvIoniz[0] = EdenFromMolecOld;
		conv.BadConvIoniz[1] = dense.EdenTrue;
	}

	/* check on molec elec den */
	if( fabs(EdenFromGrainsOld - gv.TotalEden)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "edn grn e" );
		conv.BadConvIoniz[0] = EdenFromGrainsOld;
		conv.BadConvIoniz[1] = gv.TotalEden;
	}

	/* check on heating and cooling if vary temp model */
	if( !thermal.lgTSetOn )
	{
		if( fabs(HeatingOld - thermal.htot)/thermal.htot > conv.HeatCoolRelErrorAllowed/2. )
		{
			conv.lgConvIoniz = FALSE;
			sprintf( conv.chConvIoniz, "heat chg" );
			conv.BadConvIoniz[0] = HeatingOld;
			conv.BadConvIoniz[1] = thermal.htot;
		}

		if( fabs(CoolingOld - thermal.ctot)/thermal.ctot > conv.HeatCoolRelErrorAllowed/2. )
		{
			conv.lgConvIoniz = FALSE;
			sprintf( conv.chConvIoniz, "cool chg" );
			conv.BadConvIoniz[0] = CoolingOld;
			conv.BadConvIoniz[1] = thermal.ctot;
		}
	}

	/* check on sum of grain and molec elec den - often two large num that cancel */
	if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (co.comole_eden-gv.TotalEden) )/dense.eden > 
		conv.EdenErrorAllowed/4. )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "edn grn e" );
		conv.BadConvIoniz[0] = (EdenFromMolecOld-EdenFromGrainsOld);
		conv.BadConvIoniz[1] = (co.comole_eden-gv.TotalEden);
	}

	/* check whether molecular abundances are stable - these were evaluated in CO_drive */
	for( i=0; i < NUM_COMOLE_CALC; ++i )
	{
		if( dense.lgElmtOn[co.nelem_hevmol[i]] && co.hevmol[i]>SMALLFLOAT )
		{
			if( fabs(co.hevmol[i]-comole_save[i])/dense.gas_phase[co.nelem_hevmol[i]]-1. >
				conv.EdenErrorAllowed/2.)
			{
				conv.lgConvIoniz = FALSE;
				sprintf( conv.chConvIoniz, "ch %s",co.chLab[i] );
				conv.BadConvIoniz[0] = comole_save[i]/dense.gas_phase[co.nelem_hevmol[i]];
				conv.BadConvIoniz[1] = co.hevmol[i]/dense.gas_phase[co.nelem_hevmol[i]];
			}
		}
	}

	/* check whether H molecular abundances are stable */
	for( i=0; i < N_H_MOLEC; ++i )
	{
		if( fabs(hmi.Hmolec[i]-hmole_save[i])/dense.gas_phase[ipHYDROGEN]-1. >
			conv.EdenErrorAllowed/2.)
		{
			conv.lgConvIoniz = FALSE;
			sprintf( conv.chConvIoniz, "ch %s",hmi.chLab[i] );
			conv.BadConvIoniz[0] = hmole_save[i]/dense.gas_phase[ipHYDROGEN];
			conv.BadConvIoniz[1] = hmi.Hmolec[i]/dense.gas_phase[ipHYDROGEN];
		}
	}

	/* >>chng 05 mar 26, add this convergence test - important for molecular or advective
	 * sims since iso ion solver must synch up with chemistry */
	/* remember current ground state population - will check if converged */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem<LIMELM;++nelem )
		{
			if( dense.lgElmtOn[nelem] )
			{
				/* only do this check for "significant" levels of ionization */
				/*lint -e644 var possibly not init */
				if( dense.xIonDense[nelem][nelem+1-ipISO]/dense.gas_phase[nelem] > 1e-5 )
				{
					if( fabs(iso.Pop2Ion[ipISO][nelem][0]-save_iso_grnd[ipISO][nelem])/SDIV(iso.Pop2Ion[ipISO][nelem][0])-1. >
						conv.EdenErrorAllowed)
					{
						conv.lgConvIoniz = FALSE;
						sprintf( conv.chConvIoniz,"iso %3li%li",ipISO, nelem );
						conv.BadConvIoniz[0] = save_iso_grnd[ipISO][nelem]/SDIV(iso.Pop2Ion[ipISO][nelem][0]);
						conv.BadConvIoniz[1] = iso.Pop2Ion[ipISO][nelem][0]/SDIV(iso.Pop2Ion[ipISO][nelem][0]);
					}
				}
				/*lint +e644 var possibly not init */
			}
		}
	}
	/*if( nzone >= 292 )
		fprintf(ioQQQ," DEBUG conv_base\t%.3e\t%.3e\t%.5e\t%.5e\t%.5e\t%.5e\n",
			save_iso_grnd[0][0] , iso.Pop2Ion[0][0][0],
			phycon.te , thermal.heating[0][1],
			dense.eden , iso.Pop2Ion[0][0][1]+iso.Pop2Ion[0][0][2]);*/

	/* this counts how many times ConvBase has been called in this entire model,
	 * located here at bottom of routine so that number is false on very first
	 * call, since set to 0 in zero - used to create itr/zn number in printout
	 * often used to tell whether this is the very first attempt at soln in this iteration */
	++conv.nTotalIoniz;

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

	return 0;
}

