/* 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 */
/*ion_solver solve the bi-diagonal matrix for ionization balance */
#include "cddefines.h"
#include "yield.h"
#include "prt.h"
#include "continuum.h"
#include "iso.h"
#include "dynamics.h"
#include "grainvar.h"
#include "hmi.h"
#include "mole.h"
#include "thermal.h"
#include "lapack.h"
#include "conv.h"
#include "secondaries.h"
#include "phycon.h"
#include "atmdat.h"
#include "heavy.h"
#include "elementnames.h"
#include "dense.h"
#include "radius.h"
#include "ionbal.h"

void ion_solver(
	/* this is element on the c scale, H is 0 */
	long int nelem, 
	/* option to print this element when called */
	int lgPrintIt)
{
	long int ion, 
	  limit, 
	  IonProduced, 
	  nej, 
	  ns,
	  jmax=-1;
	double dennew, rateone;
	int lgNegPop;
	static int lgMustMalloc=TRUE;
	static double *sink, *source;
	int32 nerror;
	static int32 *ipiv;
	long ion_low, ion_range, i, j, ion_to , ion_from;
	static double sum_dense;
	/* only used for debug printout */
	static double *auger;

	double abund_total, renormnew;
	int lgHomogeneous = TRUE;
	static double *xmat , *xmatsave;

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

	/* this is on the c scale, so H is 0 */
	ASSERT( nelem >= 0);
	ASSERT( dense.IonLow[nelem] >= 0 );
	ASSERT( dense.IonHigh[nelem] >= 0 );

	/* H is special because its abundance spills into three routines -
	 * the ion/atom solver (this routine), the H-mole solvers (hmole), and
	 * the heavy molecule solver.  xmolecules only includes the heavy mole
	 * part for H only.  So the difference between gas_phase and xmolecules
	 * includes the H2 part of the molecular network.  This branch does
	 * this special H case, then the general case (heavy elements are
	 * not part of the H2 molecular network) */

	/* >>chng 01 dec 07, define abund_total, total atom and ion abundance here, removing molecules */
	if( nelem == ipHYDROGEN )
	{
		/* Hydrogen is a special case since hmole does not include the
		 * atom/molecules - hmole collapses its network into H = H0 + H+
		 * and then forces the solution determined here for partitioning
		 * beweeen these two */
		abund_total = dense.xIonDense[nelem][0] + dense.xIonDense[nelem][1];
	}
	else
	{
		abund_total = SDIV( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
	}

	/* >>chng 04 nov 22,
	 * if gas nearly all atomic/ionic do not let source - sink terms from molecular network force
	 * system of balance equations to become inhomogeneous
	 * what constitutes a source or sink IS DIFFERENT for hydrogen and the rest 
	 * the H soln must couple with hmole - and its defn of source and sink.  For instance, oxygen charge
	 * transfer goes into source and sink terms for hydrogen.  So we never hose source and sink for H.
	 * for the heavy elements, which couple onto comole, mole.source and sink represent terms that
	 * remove atoms and ions from the ionization ladder.  their presence makes the system of
	 * equations inhomogeneous.  we don't want to do this when comole has a trivial effect on
	 * the ionization balance since the matrix becomes unstable */
	/* >>chng 04 dec 06, limit from 0.01 to 1e-10 as per NA suggestion */
	if( nelem>ipHYDROGEN && dense.xMolecules[nelem]/SDIV(dense.gas_phase[nelem]) < 1e-10 )
	{
		for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
		{
			mole.source[nelem][ion] = 0.;
			mole.sink[nelem][ion] = 0.;
		}
	}

	/* protect against case where all gas phase abundances are in molecules, use the
	 * atomic and first ion density from the molecule solver 
	 * >>chng 04 aug 15, NA change from 10 0000 to 10 pre-coef on
	 * FLT_EPSILON for stability in pdr */
	if( fabs( dense.gas_phase[nelem] -  dense.xMolecules[nelem])/SDIV(dense.gas_phase[nelem] ) <
		10.*FLT_EPSILON )
	{
		double renorm;
		/* >>chng 04 jul 31, add logic to conserve nuclei in fully molecular limit;
		 * in first calls, when searching for soln, we may be VERY far off, and sum of first ion
		 * and atom density may be far larger than difference between total gas and molecular densities,
		 * since they reflect the previous evaluation of the soln.  Do renorm to cover this case */
		/* first form sum of all atoms and ions */
		float sum = 0.;
		for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
			sum += dense.xIonDense[nelem][ion];
		/* now renorm to this sum - this should be unity, and is not if we have
		 * now conserved particles, due to molecular fraction changing */
		renorm = dense.gas_phase[nelem] / SDIV(sum + dense.xMolecules[nelem] );

		abund_total = renorm * sum;
	}

	/* negative residual density */
	if( abund_total < 0. )
	{
		/* >>chng 05 dec 16, do not abort when net atomic/ ionic abundance is negative due to
		 * chem net having too much of a species - this happens naturally when ices become
		 * well formed, but the code will converge away from it.  Rather, we set the ionizaiton
		 * convergence flag to "not convergee" and soldier on 
		 * if negative populations do not go away, we will eventually terminate due to
		 * convergence failures */
		/* print error if not search */
		if(!conv.lgSearch )
			fprintf(ioQQQ,
				"PROBLEM neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
				nelem,
				fabs(abund_total) / SDIV(dense.xMolecules[nelem]),
				conv.nTotalIoniz );
		/* fixup is to use half the positive abundance, assuming chem is trying to
		 * get too much of this species */
		abund_total = -abund_total/2.;

		/* say that ionization is not converged - do not abort - but if cannot converge
		 * away from negative soln, this will become a convergence failure abort */
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "neg ion" );
	}

	/* return if IonHigh is zero, since no ionization at all */
	if( dense.IonHigh[nelem] == 0 )
	{
		/* set the atom to the total gas phase abundance */
		dense.xIonDense[nelem][0] = (float)abund_total ;
#		ifdef DEBUG_FUN
		fputs( " <->ion_solver()\n", debug_fp );
#		endif
		return;
	}

	/* >>chng 01 may 09, add option to force ionization distribution with element name ioniz */
	if( dense.lgSetIoniz[nelem] )
	{
		for( ion=0; ion<nelem+2; ++ion )
		{
			dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*(float)abund_total;
		}
#		ifdef DEBUG_FUN
		fputs( " <->ion_solver()\n", debug_fp );
#		endif
		return;
	}

	/* impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1
	 * HIonFrac(nelem) is stripped to hydrogen */
	/* >>chng 01 oct 30, to assert */
	ASSERT( (dense.IonHigh[nelem] < nelem + 1) || iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. );

	/* zero out the ionization and recombination rates that we will modify here,
	 * but not the iso-electronic stages which are done elsewhere,
	 * the nelem stage of ionization is he-like,
	 * the nelem+1 stage of ionization is h-like */

	/* loop over stages of ionization that we solve for here, 
	 * up through and including one less than nelem-NISO,
	 * never actually do highest NISO stages of ionization since they
	 * come from the ionization ratio from the next lower stage */
	limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);

	/* the full range of ionization - this is number of ionization stages */
	ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
	ASSERT( ion_range <= nelem+2 );

	ion_low = dense.IonLow[nelem];

	/* PARALLEL_MODE true if there can be several instances of this routine running
	 * in the same memory pool - we must have fresh instances of the arrays for each */
	if( PARALLEL_MODE || lgMustMalloc )
	{
		long int ncell=LIMELM+1;
		lgMustMalloc = FALSE;
		if( PARALLEL_MODE )
			ncell = ion_range;

		/* this will be "new" matrix, with non-adjacent coupling included */
		if( (xmat=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ))) == NULL )
			BadMalloc();
		if( (xmatsave=(double*)MALLOC( (sizeof(double)*(unsigned)(ncell*ncell) ))) == NULL )
			BadMalloc();
		if( (sink=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ))) == NULL )
			BadMalloc();
		if( (source=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ))) == NULL )
			BadMalloc();
		if( (ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)ncell ))) == NULL )
			BadMalloc();
		/* auger is used only for debug printout - it is special because with multi-electron
		 * Auger ejection, very high stages of ionization can be produced, well beyond
		 * the highest stage considered here.  so we malloc to the limit */
		if( (auger=(double*)MALLOC( (sizeof(double)*(unsigned)(LIMELM+1) ))) == NULL )
			BadMalloc();
	}

	for( i=0; i<ion_range;i++ )
	{
		source[i] = 0.;
	}

	/* this will be used to address the 2d arrays */
#	ifdef MAT
#		undef MAT
#	endif
#	define MAT(M_,I_,J_)	(*((M_)+(I_)*(ion_range)+(J_)))

	/* zero-out loop comes before main loop since there are off-diagonal
	 * elements in the main ionization loop, due to multi-electron processes,
	 * TotIonizRate and TotRecom were already set in h-like and he-like solvers 
	 * other recombination rates were already set by routines responsible for them */
	for( ion=0; ion <= limit; ion++ )
	{
		ionbal.RateIonizTot[nelem][ion] = 0.;
	}

	/* auger is used only for debug printout - it is special because with multi-electron
	 * Auger ejection, very high stages of ionization can be produced, well beyond
	 * the highest stage considered here.  so we malloc to the limit */
	for( i=0; i< LIMELM+1; ++i )
	{
		auger[i] = 0.;
	}

	/* zero out xmat */
	for( i=0; i< ion_range; ++i )
	{
		for( j=0; j< ion_range; ++j )
		{
			MAT( xmat, i, j ) = 0.;
		}
	}

	{
		/*@-redef@*/
		/* this sets up a fake ionization balance problem, with a trival solution,
		 * for debugging the ionization solver */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nelem==ipCARBON )
		{
			broken();
			dense.IonLow[nelem] = 0;
			dense.IonHigh[nelem] = 3;
			abund_total = 1.;
			ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
			/* make up ionization and recombination rates */
			for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
			{
				double fac=1;
				if(ion)
					fac = 1e-10;
				ionbal.RateRecomTot[nelem][ion] = 100.;
				for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
				{
					/* direct photoionization of this shell */
					ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
				}
			}
		}
	}

	/* Now put in all recombination and ionization terms from CO_mole() that come from 
	   molecular reactions. this traces molecular process that change ionization stages with
	   this ladder - but do not remove from the ladder */
	for( ion_to=dense.IonLow[nelem]; ion_to <= limit; ion_to++ )
	{
		for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
		{
			/* do not do ion onto itself */
			if( ion_to != ion_from )
			{
				/* this is the rate coef for charge transfer from ion to ion_to */
				/*rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];*/
				rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to];
				ASSERT( rateone >= 0. );
				MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
				MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
			}
		}
	}
			
	/* now get actual arrays of ionization and recombination processes,
	 * but only for the ions that are done as two-level systems */
	/* in two-stage system, atom + first ion, limit is zero but must
	 * include gv.GrainChTrRate[nelem][1][0] */
	/* grain charge transfer */
	if( gv.lgDustOn && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn )
	{
		long int low;
		/* do not double count this process for atoms that are in the co network - we use
		 * a net recombination coef derived from the co solution, this includes grain ct */
		/* >>chng 05 dec 23, add mole.lgElem_in_chemistry */
		if( mole.lgElem_in_chemistry[nelem] )
		/*if( nelem==ipHYDROGEN ||nelem==ipCARBON ||nelem== ipOXYGEN ||nelem==ipSILICON ||nelem==ipSULPHUR 
			 ||nelem==ipNITROGEN ||nelem==ipCHLORINE )*/
		{
			 low = MAX2(1, dense.IonLow[nelem] );
		}
		else
			low = dense.IonLow[nelem];

		for( ion_to=low; ion_to <= limit; ion_to++ )
		{
			for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
			{
				/* do not do ion onto itself */
				if( ion_to != ion_from )
				{
					/* this is the rate coef for charge transfer from ion to ion_to */
					rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];
					MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
					MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
				}
			}
		}
	}

	for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
	{
		/* thermal & secondary collisional ionization */
		rateone = ionbal.CollIonRate_Ground[nelem][ion][0] +
			secondaries.csupra[nelem][ion] +
			/* inner shell ionization by UTA lines */
			ionbal.xInnerShellIonize[nelem][ion];;
		ionbal.RateIonizTot[nelem][ion] += rateone;

		/* UTA ionization */
		if( ion+1-ion_low < ion_range )
		{
			/* depopulation processes enter with negative sign */
			MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
			MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
		}

		/* total recombination rate */
		if( ion-1-ion_low >= 0 )
		{
			/* loss of this ion due to recom to next lower ion stage */
			MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
			MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
		}

		/* loop over all atomic sub-shells to include photoionization */
		for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
		{
			/* direct photoionization of this shell */
			ionbal.RateIonizTot[nelem][ion] += ionbal.PhotoRate_Shell[nelem][ion][ns][0];

			/* this is the primary ionization rate - add to diagonal element,
			 * test on ion stage is so that we don't include ionization from the very highest
			 * ionization stage to even higher - since those even higher stages are not considered
			 * this would appear as a sink - but populations of this highest level is ensured to
			 * be nearly trivial and neglecting it production of even higher ionization OK */
			/* >>chng 04 nov 29 RJRW, include following in this branch so only
			 * evaluated when below ions done with iso-sequence */
			if( ion+1-ion_low < ion_range )
			{
				/* this will be redistributed into charge states in following loop */ 
				MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.PhotoRate_Shell[nelem][ion][ns][0];
				
				/* yield.n_elec_eject[nelem][ion][ns] is total number of electrons that can
				 * possibly be freed 
				 * loop over nej, the number of electrons ejected including the primary,
				 * nej = 1 is primary, nej > 1 includes primary plus Auger 
				 * yield.frac_elec_eject is prob of nej electrons */
				for( nej=1; nej <= yield.n_elec_eject[nelem][ion][ns]; nej++ )
				{
					/* this is the ion that is produced by this ejection,
					 * limited by highest possible stage of ionization -
					 * do not want to ignore ionization that go beyond this */
					IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
					rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
						yield.frac_elec_eject[nelem][ion][ns][nej-1] ;
					/* >>chng 04 sep 06, above had included factor of nej to get rate, but
					 * actually want events into particular ion */
					/* number of electrons ejected
					 *(double)nej ; */
					/* it goes into this charge state - recall upper cap due to ion stage trimming 
					 * note that compensating loss term on diagonal was done before this
					 * loop, since frac_elec_eject adds to unity */
					MAT( xmat, ion-ion_low, IonProduced-ion_low ) += rateone;
					
					/* only used for possible printout - multiple electron Auger rate  -
					 * do not count one-electron as Auger */
					/*lint -e771 not init */
					if( nej>1 )
						auger[IonProduced-1] += rateone;
					/*lint +e771 not init */
				}
			}
		}

		/* this is charge transfer ionization of this species by hydrogen and helium */
		rateone = 
			atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 
			atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];
		ionbal.RateIonizTot[nelem][ion] += rateone;
		if( ion+1-ion_low < ion_range )
		{
			MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
			MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
		}
	}
	/* after this loop, ionbal.RateIonizTot and ionbal.RateRecomTot have been defined for the
	 * stages of ionization that are done with simple */
	/* begin loop at first species that is treated with full model atom */
	j = MAX2(0,limit+1);
	/* possible that lowest stage of ionization is higher than this */
	j = MAX2( ion_low , j );
	for( ion=j; ion<=dense.IonHigh[nelem] ; ion++ )
	{
		ASSERT( ion>=0 && ion<nelem+2 );
		/* use total ionization/recombination rates for species done with ISO solver */
		if( ion+1-ion_low < ion_range )
		{
			/* depopulation processes enter with negative sign */
			MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.RateIonizTot[nelem][ion];
			MAT( xmat, ion-ion_low, ion+1-ion_low ) += ionbal.RateIonizTot[nelem][ion];
		}

		if( ion-1-ion_low >= 0 )
		{
			/* loss of this ion due to recom to next lower ion stage */
			MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
			MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
		}
	}

	for( i=0; i<ion_range;i++ )
	{
		/* this will be sum of source and sink terms, will be used to decide if
		 * matrix is singular */
		double totsrc = 0.;
		ion = i+ion_low;

		/* these are the external source and sink terms */
		/* source first */
		/*src[i] += mole.source[nelem][ion];*/
		/* need negative sign to get positive pops */
		source[i] -= mole.source[nelem][ion];

		totsrc += mole.source[nelem][ion];
		/* sink next */
		/*MAT( amat, i, i ) += mole.sink[nelem][ion];*/
		MAT( xmat, i, i ) -= mole.sink[nelem][ion];
		/* matrix is not homogeneous if source is non-zero */
		if( totsrc != 0. )
			lgHomogeneous = FALSE;

	}

	/* chng 03 jan 13 rjrw, add in dynamics if required here,
	 * last test - only do advection if we have not overrun the radius scale */
	if( iteration >= 2 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth)
	{
		for( i=0; i<ion_range;i++ )
		{			 
			ion = i+ion_low;
			/*MAT( amat, i, i ) += dynamics.Rate;*/
			MAT( xmat, i, i ) -= dynamics.Rate;
			/*src[i] += dynamics.Source[nelem][ion];*/
			source[i] -= dynamics.Source[nelem][ion];
			/* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(amat,i,i),dynamics.Rate, dynamics.Source[nelem][ion]);*/
		}
		lgHomogeneous = FALSE;
	}
	for( i=0; i< ion_range; ++i )
	{
		for( j=0; j< ion_range; ++j )
		{
			MAT( xmatsave, i, j ) = MAT( xmat, i, j );
		}
	}

	/* this is true if no source terms */
	if( lgHomogeneous )
	{
		double ionprod=1., recomb /*, scale = 0.*/;
		/* Simple estimate of most abundant ion */
		jmax = 0;
		for( i=0; i<ion_range-1;i++)
		{ 
			ion = i+ion_low;
			ionprod *= ionbal.RateIonizTot[nelem][ion];
			recomb = ionbal.RateRecomTot[nelem][ion];
			/* trips if ion rate zero, so ll the gas will be more neutral than this */
			if( ionprod == 0)
				break;
			/* rec rate is zero */
			if (recomb <= 0.) 
				break;

			ionprod /= recomb;
			if (ionprod > 1.) 
			{
				/* this is peak ionization stage */
				jmax = i;
				ionprod = 1.;
			}
		}
		
		for( i=0; i<ion_range;i++ )
		{
			MAT(xmat,i,jmax) = 1.;
		}
		source[jmax] = abund_total;
	}

	if ( FALSE && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 ) 
	{
		fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
		fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot[nelem][0], ionbal.RateIonizTot[nelem][1]);
		fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
		fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
	}

	/* get new soln */
	nerror = 0;
	/* Use general matrix solver */
	getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror);
	if( nerror != 0 )
	{
		fprintf( ioQQQ, 
			" ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li ion_range=%li, limit=%li, xmat follows\n",
			nelem , ion_range,limit  );
		for( i=0; i<ion_range; ++i )
		{
			for( j=0;j<ion_range;j++ )
			{
				fprintf(ioQQQ,"%e\t",MAT(xmat,j,i));
			}
			fprintf(ioQQQ,"\n");
		}
		fprintf(ioQQQ,"source follows\n");
		for( i=0; i<ion_range;i++ )
		{
			fprintf(ioQQQ,"%e\t",source[i]);
		}
		fprintf(ioQQQ,"\n");
		puts( "[Stop in ion_solver]" );
		cdEXIT(EXIT_FAILURE);
	}
	getrs_wrapper('N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror);
	if( nerror != 0 )
	{
		fprintf( ioQQQ, " ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
			nelem , ion_range );
		puts( "[Stop in ion_solver]" );
		cdEXIT(EXIT_FAILURE);
	}

	{
		/*@-redef@*/
		/* this is to debug following failed assert */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone >380 ) && nelem == ipHYDROGEN )
		{
			fprintf(ioQQQ,"debuggg\t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n", 
				fnzone,
				phycon.te,
				dense.eden,
				ionbal.RateIonizTot[nelem][0] , 
				ionbal.RateRecomTot[nelem][0]);
			fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[ipHYDROGEN][0], mole.source[ipHYDROGEN][1]);
			fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[ipHYDROGEN][0], mole.sink[ipHYDROGEN][1]);
			fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0],
				ionbal.RateIonizTot[nelem][0]/ionbal.RateRecomTot[nelem][0]);
		}
	}


	/* 
	 * >> chng 03 jan 15 rjrw:- terms are now included for
	 * molecular sources and sinks of H and H+.
	 *
	 * When the network is not in equilibrium, this will lead to a
	 * change in the derived abundance of H and H+ when after the
	 * matrix solution -- the difference between `renorm' and 1. is a
	 * measure of the quality of the solution (it will be 1. if the
	 * rate of transfer into Ho/H+ balances the rate of transfer
	 * out, for the consistent relative abundances).
	 *
	 * We therefore renormalize to keep the total H abundance
	 * correct -- only the molecular network is allowed to change
	 * this.
	 *
	 * To do this, only the ion abundances are corrected, as the
	 * molecular abundances may depend on several different
	 * conserved species.
	 *
	 */
	if( lgHomogeneous )
	{
		dense.xIonDense[nelem][ion_low] = (float)abund_total;
		for ( i=1;i < ion_range; i++ )
		{
			dense.xIonDense[nelem][i+ion_low] = 0.;
		}
	}

	renormnew = 1.;
	if (iteration >= 2 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth && 
			nelem == ipHYDROGEN && hmi.lgNoH2Mole)
	{
		/* The normalization out of the matrix solution is correct and
		 * should be retained if: dynamics is on and the total
		 * abundance of HI & H+ isn't being controlled by the
		 * molecular network */
		renormnew = 1.;
	}
	else
	{
		dennew = 0.;
		sum_dense = 0.;
		
		/* find total population to renorm - also here check that negative pops did not occur */
		for( i=0;i < ion_range; i++ )
		{
			ion = i+ion_low;
			sum_dense += dense.xIonDense[nelem][ion];
			dennew += source[i];
		} 
		
		if( dennew > 0.)
		{
			renormnew = sum_dense / dennew;
			/*TODO	2	renorm should == 1 when the molecules and
			 * ionization are in equilibrium.  Should monitor
			 * this figure of merit in calling routine.
			 * */
		}
		else
		{
			renormnew = 1.;
		}
	}
	/* check not negative, should be +ve, can be zero if species has become totally molecular.
	 * this happens for hydrogen if no cosmic rays, or cr ion set very low */
	if( renormnew < 0)
	{
		fprintf(ioQQQ,"PROBLEM impossible value of renorm \n");
	}
	ASSERT( renormnew>=0 );

	/* save resulting abundances into main ionization density array, 
	 * while checking whether any negative level populations occured */
	lgNegPop = FALSE;
	for( i=0; i < ion_range; i++ )
	{
		ion = i+ion_low;

		/*fprintf(ioQQQ," %li %li %.3e %.3e\n",nelem,ion,src[ion-ion_low+1],src[ion-ion_low]);
			pop_ion_ov_neut[ion] = src[ion-ion_low+1]/src[ion-ion_low]; */
		if( source[i] < 0. )
		{
			/* >>chng 04 dec 04, put test on neg abund here, don't print unles value is very -ve */
			if( source[i]<-1e-10 )
				fprintf(ioQQQ,
				" PROBLEM negative ion population in ion_solver, nelem=%li ion=%li val=%.3e Search?%c\n",
				nelem , ion , source[i] , TorF(conv.lgSearch) );
			source[i] = 0.;
			/* if this is one of the iso seq model atoms then must also zero out Lya and Pop2Ion */
			if( ion == nelem+1-NISO )
			{
				long int ipISO = nelem - ion;
				ASSERT( ipISO>=0 && ipISO<NISO );
				iso.Pop2Ion[ipISO][nelem][0] = 0.;
			}
		}

		/* use new solution */
		dense.xIonDense[nelem][ion] = (float)(source[i]*renormnew);

	}

	/* Zero levels with abundances < 1e-25 which which will suffer numerical noise */
	while( dense.IonHigh[nelem] > dense.IonLow[nelem] && 
		dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total )
	{
		ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
		/* zero out abundance and heating due to stage of ionization we are about to zero out */
		dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
		thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
		/* decrement counter */
		--dense.IonHigh[nelem];
	}

	/* sanity check, either offset stages of low and high ionization,
	 * or no ionization at all */
	ASSERT( (dense.IonLow[nelem] < dense.IonHigh[nelem]) ||
		(dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) );

	/* this should not happen */
	if( lgNegPop )
	{
		fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n", 
		  elementnames.chElementNameShort[nelem], nzone );

		fprintf( ioQQQ, " Populations were" );
		for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
		{
			fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, " destroy vector =" );
		for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot[nelem][ion-1] );
		}
		fprintf( ioQQQ, "\n" );

		/* print some extra stuff if destroy was negative */
		if( lgNegPop )
		{
			fprintf( ioQQQ, " CTHeavy  vector =" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%9.1e", atmdat.HeCharExcIonOf[nelem][ion] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " HCharExcIonOf vtr=" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%9.1e", atmdat.HCharExcIonOf[nelem][ion] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " CollidRate  vtr=" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
			}
			fprintf( ioQQQ, "\n" );

			/* photo rates per subshell */
			fprintf( ioQQQ, " photo rates per subshell, ion\n" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%3ld", ion );
				for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
				{
					fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
				}
				fprintf( ioQQQ, "\n" );
			}
		}

		/* now check out creation vector */
		fprintf( ioQQQ, " create  vector =" );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		ContNegative();
		ShowMe();
		puts( "[Stop in bidiag]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* option to print ionization and recombination arrays
	 * prt flag set with print array print arrays command */
	if( prt.lgPrtArry || lgPrintIt )
	{
		/* say who we are, what we are doing .... */
		fprintf( ioQQQ, 
			"\n %s ion_solver DEBUG ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n", 
			elementnames.chElementSym[nelem],
			elementnames.chElementName[nelem],
			fnzone,
			phycon.te , 
			dense.eden,
			dense.gas_phase[nelem],
			abund_total ,
			dense.xMolecules[nelem] );
		/* total ionization rate, all processes */
		fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.RateIonizTot[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* sinks from the chemistry network */
		fprintf(ioQQQ," %s mole sink   ",
			elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf(ioQQQ," %9.2e", mole.sink[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* sum of all creation processes */
		fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.RateRecomTot[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* sources from the chemistry network */
		fprintf(ioQQQ," %s mole source ",
			elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf(ioQQQ," %9.2e", mole.source[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* collisional ionization */
		fprintf( ioQQQ, " %s Coll ioniz  " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.CollIonRate_Ground[nelem][ion][0] );
		}
		fprintf( ioQQQ, "\n" );

		/* UTA ionization */
		fprintf( ioQQQ, " %s UTA ioniz   " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.xInnerShellIonize[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* photo ionization */
		fprintf( ioQQQ, " %s Phot ioniz  " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", 
				ionbal.PhotoRate_Shell[nelem][ion][Heavy.nsShells[nelem][ion]-1][0] );
		}
		fprintf( ioQQQ, "\n" );

		/* auger ionization */
		fprintf( ioQQQ, " %s Auger ioniz " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", 
				auger[ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* secondary ionization */
		fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", 
				secondaries.csupra[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* grain ionization - not total rate but should be dominant process */
		fprintf( ioQQQ, " %s ion on grn  "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion][ion+1] );
		}
		fprintf( ioQQQ, "\n" );

		/* charge exchange ionization */
		fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			double sum = atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 
				atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];

			if( nelem==ipHELIUM && ion==0 )
			{
				sum += atmdat.HeCharExcIonTotal;
			}
			else if( nelem==ipHYDROGEN && ion==0 )
			{
				sum += atmdat.HCharExcIonTotal;
			}
			fprintf( ioQQQ, " %9.2e", sum );
		}
		fprintf( ioQQQ, "\n" );

		/* radiative recombination */
		fprintf( ioQQQ, " %s radiati rec "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.RR_rate_coef_used[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* old DR recombination */
		fprintf( ioQQQ, " %s dr old  rec "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_old_rate_coef[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* Badnell DR recombination */
		fprintf( ioQQQ, " %s drBadnel rec"  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* grain recombination - not total but from next higher ion, should
		 * be dominant */
		fprintf( ioQQQ, " %s form on grn "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", gv.GrainChTrRate[nelem][ion+1][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* charge exchange recombination */
		fprintf( ioQQQ, " %s chr trn rec "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			double sum = 
				atmdat.HCharExcRecTo[nelem][ion]*
				iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*
				dense.xIonDense[ipHYDROGEN][1] +

				atmdat.HeCharExcRecTo[nelem][ion]*
				iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe1s1S]*
				dense.xIonDense[ipHELIUM][1];

			if( nelem==ipHELIUM && ion==0 )
			{
				sum += atmdat.HeCharExcRecTotal;
			}
			else if( nelem==ipHYDROGEN && ion==0 )
			{
				sum += atmdat.HCharExcRecTotal;
			}
			fprintf( ioQQQ, " %9.2e", sum );
		}
		fprintf( ioQQQ, "\n" );

		/* the "new" abundances the resulted from the previous ratio */
		fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", dense.xIonDense[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );
	}
	
	if( PARALLEL_MODE )
	{
		free(ipiv);
		free(source);
		free(sink);
		free(xmat);
		free(xmatsave);
		free(auger);
	}

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

