/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*CO_step fills in matrix for heavy elements molecular routines */
#include "cddefines.h"
#include "called.h"
#include "dense.h"
#include "deuterium.h"
#include "ionbal.h"
#include "thermal.h"
#include "phycon.h"
#include "hmi.h"
#include "dynamics.h"
#include "conv.h"
#include "trace.h"
#include "timesc.h"
#include "mole.h"
#include "mole_priv.h"
#include "grainvar.h"
#include "h2.h"
#include "newton_step.h"
/* Nick Abel between July and October of 2003 assisted Dr. Ferland in
 * improving the heavy element molecular network in Cloudy. Before
 * this routine would predict negative abundances if the fraction of
 * carbon in the form of molecules came close to 100%. A reorganizing
 * of the reaction network detected several bugs.  Treatment of
 * "coupled reactions", in which both densities in the reaction rate
 * were being predicted by Cloudy, were also added.  Due to these
 * improvements, Cloudy can now perform calculations where 100% of the
 * carbon is in the form of CO without predicting negative abundances
 *
 * Additional changes were made in November of 2003 so that our reaction 
 * network would include all reactions from the TH85 paper.  This involved 
 * adding silicon to the chemical network.  Also the reaction rates were
 * labeled to make identification with the reaction easier and the matrix 
 * elements of atomic C, O, and Si are now done in a loop, which makes 
 * the addition of future chemical species (like N or S) easy.
 * */
void check_co_ion_converge(void);

STATIC void funjac(GroupMap &MoleMap, const valarray<double> &b2vec,
						 double * const ervals, double * const amat, const bool lgJac, bool *lgConserve);

STATIC void mole_h_fixup(void);

STATIC void grouped_elems(const double bvec[], double mole_elems[]);

#define SMALLABUND 1e-24

double mole_solve()
{
	int nBad, nPrevBad, i;
	realnum
		eqerror, error;
	valarray<double> oldmols(mole_global.num_compacted),
		newmols(mole_global.num_compacted);
	GroupMap MoleMap( atom_list.size() );

	DEBUG_ENTRY( "mole_solve()" );

	if (hmi.H2_frac_abund_set>0.)
	{
		mole_h_fixup();
		fixit(); 
		fprintf(stderr,"Need to treat hmi.H2_frac_abund_set in revised network\n");
		fprintf(stderr,"%g\n",hmi.H2_frac_abund_set);
		exit(-1);
	}

	ASSERT(lgElemsConserved());

	/* will test against this error for quitting */
	const double BIGERROR = 1e-8;
	/* >>chng 04 jul 20, upper limit had been 6, why?  change to 20
	 * to get closer to soln */
	const int LIM_LOOP = 100;
	/* loop until run to limit, or no fix ups needed and error is small */
	/* will be used to keep track of neg solns */
	nBad = nPrevBad = 0;
	eqerror = -1.;

	valarray<double> escale(mole_global.num_compacted);
	double rlimit=-1., rmax=0.0;
		
	// Run enough times through to converge nonlinear system.
	for(i=0; i < LIM_LOOP;i++) 
	{
		{
			/* option to print deduced abundances */
			/*@-redef@*/
			enum {DEBUG_LOC=false};
			/*@+redef@*/
			if( DEBUG_LOC && (nzone>140) )
			{
				fprintf(ioQQQ,"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
					fnzone,
					phycon.te,
					dense.eden);
				for( long k=0; k<mole_global.num_calc; k++ )
					fprintf(ioQQQ,"\t%.2e", mole.species[k].den );
				fprintf(ioQQQ,"\n" );
			}
		}
		/* nBad returns number of times negative abundances were fixed
		 * for latest step, should be zero at return for valid soln */

		nPrevBad = nBad;

		MoleMap.setup(get_ptr(oldmols));

		long j;

		// Need to allow enough loops for j that rlimit increases until
		// no species get to -ve densities.  Often, though, this will
		// just be the first time through.
		conv.incrementCounter(MOLE_SOLVE);
		for (j=0; j<30; j++) 
		{
			conv.incrementCounter(MOLE_SOLVE_STEPS);
			bool lgOK = newton_step(MoleMap, oldmols, newmols, &eqerror, &error, mole_global.num_compacted, &rlimit, &rmax, escale, funjac); 
		
			/* check for negative populations */
			if (lgOK)
			{
				nBad = 0; 
				double fracneg = 0.;
				long iworst = -1;
				for( long mol=0; mol < mole_global.num_compacted; mol++ )
				{
					if( newmols[mol] < 0. ) 
					{
						if(-newmols[mol] > fracneg*SDIV(oldmols[mol]))
						{
							fracneg = -newmols[mol]/SDIV(oldmols[mol]);
							iworst = mol;
						}
						if( newmols[mol] < -SMALLABUND) 
						{
							++nBad;
						}
						newmols[mol] = 0.;
					}
				}
				if ( 0 != nBad)
				{
					lgOK = false;
					if (0)
					{
						fprintf(ioQQQ,"-ve density in inner chemistry loop, worst is %ld\n",iworst);
					}
				}
			}
			if ( lgOK )
			{
				//fprintf(ioQQQ,"OK at z %ld l %d j %ld t %g e %g\n",nzone,i,j,1./rlimit,eqerror);
				break;
			}

			ASSERT( rlimit < BIGFLOAT );
			ASSERT( rlimit > 0.0 );
			rlimit = 10.*rlimit;
		}
		
		//fprintf(stdout,"Mole zone %ld -- %7.2f error %15.8g eqerror %15.8g rlimit %15.8g nBad %d ninner %ld loop %d\n",nzone,fnzone,error,eqerror,rlimit,nBad,j+1,i);

		// If the accuracy of the solution obtained was strongly limited
		// by the pseudo-timestep limit, extend it
		if ( error < 0.01*eqerror )
			rlimit = 0.1*rlimit;
		
		MoleMap.updateMolecules( newmols );
			
		/* Stop early if good enough */
		if (eqerror < BIGERROR && nBad == 0 && nPrevBad == 0) 
			break;
	}

	//fprintf(stdout,"Mole: zn %ld -- %7.2f err %13.6g rlimit %13.6g nBad %d loop %d\n",nzone,fnzone,eqerror,rlimit,nBad,i);

	if( (i == LIM_LOOP && eqerror > BIGERROR)  || nBad != 0 )
	{
		conv.setConvIonizFail("Chem conv", eqerror, nBad);
	}
	
	// Effect_delta determines whether changes due to the molecular
	// network rescale values enough to invalidate other solutions
	realnum effect_delta = mole_return_cached_species(MoleMap);

	realnum delta_threshold = 0.2*conv.IonizErrorAllowed;
	if (effect_delta > delta_threshold)
	{
		conv.setConvIonizFail("chem eft chg", effect_delta, 0.0);
	}		

	if (trace.lgTrace)
	{
		fprintf(ioQQQ,"Mole zn %3ld -- %7.2f\n",nzone,fnzone); 
		fprintf(ioQQQ,"Internal error %15.8g nBad %4d loops %3d\n",
				  eqerror,nBad,i);
		fprintf(ioQQQ,"Scaling delta on ions %15.8g; threshold %15.8g\n", 
				  effect_delta, delta_threshold);
	}

	check_co_ion_converge();

	{
		/* option to print deduced abundances */
		/*@-redef@*/
		enum {DEBUG_LOC=false};
		/*@+redef@*/
		if( DEBUG_LOC /*&& (nzone>68)*/ )
		{
			fprintf(ioQQQ,"DEBUG mole out\t%i\t%.2f\t%.5e\t%.5e",
							i,
							fnzone,
							phycon.te,
							dense.eden);
			fprintf(ioQQQ,
							"\terror\t%.4e\tH0\t%.4e\tH+\t%.4e\tsink\t%.4e\t%.4e\tsour\t%.4e\t%.4e\tion\t%.4e\trec\t%.4e", 
							eqerror,
							dense.xIonDense[ipHYDROGEN][0],
							dense.xIonDense[ipHYDROGEN][1],
							mole.sink[ipHYDROGEN][0],
							mole.sink[ipHYDROGEN][1],
							mole.source[ipHYDROGEN][0] , 
							mole.source[ipHYDROGEN][1] ,
							ionbal.RateIonizTot(ipHYDROGEN,0),
							ionbal.RateRecomTot[ipHYDROGEN][0]);
			
			/*for( j=0; j<mole_global.num_calc; j++ )
				fprintf(ioQQQ,"\t%.4e", HMOLEC(j) );*/
			fprintf(ioQQQ,"\n" );
		}
	}
	
	return eqerror;

/*lint +e550 */
/*lint +e778 constant expression evaluates to 0 in operation '-' */
}


void check_co_ion_converge(void)
{
	DEBUG_ENTRY("check_co_ion_converge()");
	/* check whether ion and chem solvers agree yet */
	if( dense.lgElmtOn[ipCARBON] &&
		fabs(dense.xIonDense[ipCARBON][0]- findspecieslocal("C")->den)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
	{
		conv.setConvIonizFail("CO C0 con",
									 dense.xIonDense[ipCARBON][0],
									 findspecieslocal("C")->den);
	}
	else if( dense.lgElmtOn[ipCARBON] &&
		fabs(dense.xIonDense[ipCARBON][1]- findspecieslocal("C+")->den)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
	{
		conv.setConvIonizFail("CO C1 con",
									 dense.xIonDense[ipCARBON][1],
									 findspecieslocal("C+")->den);
	}
	else if( dense.lgElmtOn[ipOXYGEN] &&
		fabs(dense.xIonDense[ipOXYGEN][0]- findspecieslocal("O")->den)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
	{
		conv.setConvIonizFail(
			"CO O0 con",dense.xIonDense[ipOXYGEN][0],
			findspecieslocal("O")->den);
	}
	else if( dense.lgElmtOn[ipOXYGEN] &&
		fabs(dense.xIonDense[ipOXYGEN][1]- findspecieslocal("O+")->den)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
	{
		conv.setConvIonizFail(
			"CO O1 con", dense.xIonDense[ipOXYGEN][1],
			findspecieslocal("O+")->den);
	}
}

STATIC void mole_h_fixup(void)
{
	long int mol;
	
	/* there are two "no molecules" options, the no co, which turns off everything,
	 * and the no n2, which only turns off the h2.  in order to not kill the co
	 * part we still need to compute the hydrogen network here, and then set h2 to
	 * small values */
	/* >> chng 03 jan 15 rjrw -- suddenly switching off molecules confuses the solvers... */
	DEBUG_ENTRY( "mole_h_fixup()" );

	/* >>chng 02 jun 19, add option to force H2 abundance, for testing h2 molecules,
	 * hmi.H2_frac_abund_set is fraction in molecules that is set by set h2 fraction command */
	if( hmi.H2_frac_abund_set>0.)
	{
		for(mol=0;mol<mole_global.num_calc;mol++) 
		{
			mole.species[mol].den = 0.;
		}
		/* >>chng 03 jul 19, from 0 to SMALLFLOAT, to pass asserts in ConvBase,
		 * problem is that ion range has not been reset for hydrogen */
		dense.xIonDense[ipHYDROGEN][0] = dense.xIonDense[ipHYDROGEN][1] = 
			2.f*SMALLFLOAT*dense.gas_phase[ipHYDROGEN];
		/* put it all in the ground state */
		findspecieslocal("H2")->den = (realnum)(dense.gas_phase[ipHYDROGEN] * hmi.H2_frac_abund_set);
		findspecieslocal("H2*")->den = 0.;

		hmi.H2_total = findspecieslocal("H2")->den + findspecieslocal("H2*")->den;
		/* first guess at ortho and para densities */
		h2.ortho_density = 0.75*hmi.H2_total;
		h2.para_density = 0.25*hmi.H2_total;
		{
			hmi.H2_total_f = (realnum)hmi.H2_total;
			h2.ortho_density_f = (realnum)h2.ortho_density;
			h2.para_density_f = (realnum)h2.para_density;
		}

		hmi.hmihet = 0.;
		hmi.h2plus_heat = 0.;
		hmi.H2Opacity = 0.;
		hmi.hmicol = 0.;
		hmi.HeatH2Dish_TH85 = 0.;
		hmi.HeatH2Dexc_TH85 = 0.;
		hmi.deriv_HeatH2Dexc_TH85 = 0.;
		hmi.hmidep = 1.;

		for( size_t nd=0; nd < gv.bin.size(); nd++ )
		{
			gv.bin[nd]->rate_h2_form_grains_used = 0.;
		}

		return;
	}
}

#define ABSLIM  1e-12
#define ERRLIM  1e-12
#	ifdef MAT
#		undef MAT
#	endif
#	define MAT(a,I_,J_)	((a)[(I_)*(mole_global.num_compacted)+(J_)])


STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c);
enum {PRINTSOL = false};

STATIC void funjac(GroupMap &MoleMap, const valarray<double> &b2vec, double * const ervals,
						 double * const amat, const bool lgJac, bool *lgConserve)
{
	static multi_arr<double,2> c(mole_global.num_total,mole_global.num_total);
	bool printsol = PRINTSOL;
	valarray<double> b(mole_global.num_total);

	DEBUG_ENTRY( "funjac()" );

	MoleMap.updateMolecules( b2vec );		
		
	if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection && 
			dynamics.Rate != 0.0) 
	{
		ASSERT(dynamics.Rate > 0.0);
		*lgConserve = false; 
	}
	else
	{
		*lgConserve = true; 
	}		

	/* Generate chemical balance vector (mole.b[]) and Jacobian array 
		 (c[][], first iteration) from reaction list */
	
	mole_eval_dynamic_balance(mole_global.num_total,get_ptr(b),lgJac,c);	
	
	/*------------------------------------------------------------------ */
	if(printsol || (trace.lgTrace && trace.lgTraceMole ))
	{
		/* print the full matrix */
		fprintf( ioQQQ, "                ");
		for( long i=0; i < mole_global.num_calc; i++ )
		{
			fprintf( ioQQQ, "      %-4.4s", mole_global.list[i]->label.c_str() );
		}
		fprintf( ioQQQ, " \n" );
		
		fprintf(ioQQQ,"       MOLE old abundances\t%.2f",fnzone);
		for( long i=0; i<mole_global.num_calc; i++ )
			fprintf(ioQQQ,"\t%.2e", mole.species[i].den );
		fprintf(ioQQQ,"\n" );
		
		for( long i=0; i < mole_global.num_calc; i++ )
		{
			fprintf( ioQQQ, "       MOLE%2ld %-4.4s", i ,mole_global.list[i]->label.c_str());
			for( long j=0; j < mole_global.num_calc; j++ )
			{
				fprintf( ioQQQ, "%10.2e", c[j][i] );
			}
			fprintf( ioQQQ, "%10.2e", b[i] );
			fprintf( ioQQQ, "\n" );
		}
	}

	/* add positive ions and neutral atoms: ratios are set by ion_solver,
	 * we determine abundance of the group as a whole here */
	
	if (lgJac) {
		for (unsigned long j=0; j<atom_list.size(); ++j )
		{
			if (atom_list[j]->ipMl[0] != -1)
			{
				for(long i=0;i<mole_global.num_calc;i++) 
				{
					ASSERT(mole_global.list[i]->index == i);
					double sum = 0.;
					for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++) 
					{
						if (atom_list[j]->ipMl[ion] != -1)
						{
							sum += MoleMap.fion[j][ion]*c[atom_list[j]->ipMl[ion]][i];
						}
					}
					c[atom_list[j]->ipMl[0]][i] = sum;
				}
			}
		}
		for (unsigned long j=0; j<atom_list.size(); ++j )
		{
			if (atom_list[j]->ipMl[0] != -1)
			{
				for(long i=0;i<mole_global.num_calc;i++) 
				{
					double sum = 0.0;

					for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
					{
						if (atom_list[j]->ipMl[ion] != -1)
						{
							sum += c[i][atom_list[j]->ipMl[ion]];
						}
					}
					c[i][atom_list[j]->ipMl[0]] = sum; 
				}					
			}
		}
	}
	
	for (unsigned long j=0; j<atom_list.size(); ++j )
	{
		//if (j == 8)
		//	fprintf(ioQQQ,"Nsum1 %s %ld %d %d\n",
		//			  atom_list[j]->label().c_str(), j, groupspecies[29]->index,atom_list[j]->ipMl[0]);
		if (atom_list[j]->ipMl[0] != -1)
		{
			double sum = 0.0;
			for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
			{	
				if (atom_list[j]->ipMl[ion] != -1)
				{
					//			if (j == 8)
					//	fprintf(ioQQQ,"Nsum %d %15.8g\n",atom_list[j]->ipMl[ion],
					//			  b[atom_list[j]->ipMl[ion]]);
					sum += b[atom_list[j]->ipMl[ion]];
					b[atom_list[j]->ipMl[ion]] = 0.0;
				}
			}
			b[atom_list[j]->ipMl[0]] = sum;
		}
	}
	
	/* Species are now grouped -- only mole_global.num_compacted elements now active */

	if (*lgConserve)
	{
		// Replace atom rows with conservation constraints
		if (lgJac)
		{
			for ( unsigned long j = 0; j < atom_list.size(); ++j )
			{
				long ncons = atom_list[j]->ipMl[0];
				if (ncons != -1)
				{
					ASSERT( mole_global.list[ncons]->isMonatomic() );
					double scale=1.0/MoleMap.molElems[j]; //fabs(c[ncons][ncons]);// 
					for( long i=0;i<mole_global.num_compacted;i++)
					{
						if( groupspecies[i]->nAtom.find(atom_list[j]) != groupspecies[i]->nAtom.end() )
							c[groupspecies[i]->index][ncons] = groupspecies[i]->nAtom[atom_list[j]]*scale;
						else
							c[groupspecies[i]->index][ncons] = 0.;
						
					}
				}
			}
		}

		valarray<double> molnow(atom_list.size());
		grouped_elems(get_ptr(b2vec), get_ptr(molnow));
		for ( unsigned long j = 0; j < atom_list.size(); ++j )
		{
			long ncons = atom_list[j]->ipMl[0];
			if (ncons != -1)
			{
				ASSERT( mole_global.list[ncons]->isMonatomic() );
				double scale = c[ncons][ncons];
				b[ncons] = (molnow[j]-MoleMap.molElems[j])*scale;
				if (false)
					if (b[ncons] != 0.0)
						fprintf(ioQQQ,"Cons %s err %g rel %g\n",
								  atom_list[j]->label().c_str(),b[ncons],
								  (molnow[j]-MoleMap.molElems[j])/SDIV(MoleMap.molElems[j]));
			}
		}
	}
	
	for( long i=0; i < mole_global.num_compacted; i++ )
	{
		ervals[i] = b[groupspecies[i]->index];
	}
	
	if (lgJac)
	{

		for( long i=0; i < mole_global.num_compacted; i++ )
		{
			for( long j=0; j < mole_global.num_compacted; j++ )
			{
				MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
			}			 
		}

		// Replace rows for species with no sources and sinks with
		// an identity
		for(long i=0;i<mole_global.num_compacted;i++)
		{
			double sum1 = 0.;
			for (long j=0;j<mole_global.num_compacted;j++)
			{
				sum1 += fabs(MAT(amat,j,i));
			}
			if (sum1 == 0.0)
			{
				ASSERT(ervals[i] == 0.0);
				MAT(amat,i,i) = 1.0;
				// fprintf(ioQQQ,"Fixing %ld %s\n",i,groupspecies[i]->label.c_str());
			}
		}
	}
}

STATIC void grouped_elems(const double bvec[], double mole_elems[])
{
	map<chem_atom*, long> atom_to_index;
	for (unsigned long j=0; j<atom_list.size(); ++j )
	{
		mole_elems[j] = 0.;
		atom_to_index[atom_list[j].get_ptr()] = j;
	}
	for( long i=0; i < mole_global.num_compacted; i++ )
	{
		if( groupspecies[i]->parentLabel.empty() == false )
			continue;
		
		for (molecule::nAtomsMap::const_iterator el = groupspecies[i]->nAtom.begin(); 
			  el != groupspecies[i]->nAtom.end(); ++el)
		{
			mole_elems[atom_to_index[el->first.get_ptr()]] += bvec[i]*el->second;
		}
	}
}

void GroupMap::setup(double *b0vec)
{
	valarray<double> calcv(mole_global.num_total);
	bool lgSet;

	for( long i=0;i<mole_global.num_total;i++)
	{
		calcv[i] = mole.species[i].den;
	}

	for (unsigned long j=0; j<atom_list.size(); ++j )
	{
		if (atom_list[j]->ipMl[0] != -1)
		{
			double sum = 0.;
			for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
			{
				if (atom_list[j]->ipMl[ion] != -1)
					sum += calcv[atom_list[j]->ipMl[ion]];
			}
			if (sum > SMALLFLOAT)
			{
				double factor = 1./sum;
				for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
				{
					if (atom_list[j]->ipMl[ion] != -1)
						fion[j][ion] = calcv[atom_list[j]->ipMl[ion]]*factor;
					else
						fion[j][ion] = 0.;
				}
			}
			else
			{
				lgSet = false;
				for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
				{
					if (atom_list[j]->ipMl[ion] != -1 && !lgSet)
					{
						fion[j][ion] = 1.0;
						lgSet = true;
					}
					else
					{
						fion[j][ion] = 0.;
					}
				}
			}
			
			lgSet = false;
			for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
			{
				if (atom_list[j]->ipMl[ion] != -1)
				{
					if (!lgSet) 
						calcv[atom_list[j]->ipMl[ion]] = sum;
					else
						calcv[atom_list[j]->ipMl[ion]] = 0.;
					lgSet = true;
				}
			}
		}
	}	

	for( long i=0; i < mole_global.num_compacted; i++ )
	{
		b0vec[i] = calcv[groupspecies[i]->index];
	}
	
	grouped_elems(get_ptr(b0vec),get_ptr(molElems));

	for (unsigned long i = 0; i<atom_list.size(); ++i)
	{
		double densAtom = 0.;
		// deuterium is special case
		if( atom_list[i]->el->Z==1 && atom_list[i]->A==2 )
		{
			ASSERT( deut.lgElmtOn );
			densAtom = deut.gas_phase;
		}
		// skip other isotopes
		else if( !atom_list[i]->lgMeanAbundance() )
			continue;
		else
		{
			int nelem = atom_list[i]->el->Z-1;
			densAtom = dense.gas_phase[nelem];
		}
		bool lgTest = 
			( densAtom < SMALLABUND && molElems[i] < SMALLABUND ) ||
			( fabs(molElems[i]- densAtom) <= conv.GasPhaseAbundErrorAllowed*densAtom );
		if( !lgTest )
			fprintf( ioQQQ, "PROBLEM molElems BAD  %li\t%s\t%.12e\t%.12e\t%.12e\n", 
						i, atom_list[i]->label().c_str(), atom_list[i]->frac, densAtom, molElems[i] );
		
		ASSERT( lgTest );
		molElems[i] = densAtom;
	}

}

void GroupMap::updateMolecules(const valarray<double> & b2 )
{
	DEBUG_ENTRY("updateMolecules()");

	for (long mol=0;mol<mole_global.num_calc;mol++)
	{
		mole.species[mol].den = 0.;
	}
	for (long mol=0;mol<mole_global.num_compacted;mol++)
	{
		mole.species[ groupspecies[mol]->index ].den = b2[mol];	/* put derived abundances back into appropriate molecular species */
	}

	// calculate isotopologue densities
	for (long mol=0;mol<mole_global.num_calc;mol++)
	{
		if( mole_global.list[mol]->parentIndex >= 0 )
		{
			ASSERT( !mole_global.list[mol]->parentLabel.empty() );
			long parentIndex = mole_global.list[mol]->parentIndex;
			mole.species[mol].den = mole.species[parentIndex].den;
			for( nAtoms_i it = mole_global.list[mol]->nAtom.begin(); it != mole_global.list[mol]->nAtom.end(); ++it )
			{
				if( !it->first->lgMeanAbundance() )
				{
					mole.species[mol].den *= pow( it->first->frac, it->second );
				}
			} 
		}
	}
	
	// calculate densities for monatomic species that had been collapsed into a group
	for (unsigned long j=0; j<atom_list.size(); ++j )
	{
		if (atom_list[j]->ipMl[0] != -1)
		{
			double grouptot = mole.species[atom_list[j]->ipMl[0]].den;
			double sum = 0.0;
			for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
			{
				if (atom_list[j]->ipMl[ion] != -1)
				{
					mole.species[atom_list[j]->ipMl[ion]].den = grouptot * fion[j][ion];
					sum += mole.species[atom_list[j]->ipMl[ion]].den;
				}
			}
			ASSERT(fabs(sum-grouptot) <= 1e-10 * fabs(grouptot));
		}
	}
	
	mole.set_isotope_abundances();

	return;
}

STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c)
{
	double source;
	
	DEBUG_ENTRY("mole_eval_dynamic_balance()");

	mole_eval_balance(num_total, b, lgJac, c);

	/* >>chng 06 mar 17, comment out test on old full depth - keep old solution if overrun scale */
	if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
			dynamics.Rate != 0.0 ) 
	{
		/* Don't use conservation form in matrix solution -- dynamics rate terms make c[][] non-singular */		

		for( long i=0;i<mole_global.num_calc;i++)
		{
			if (lgJac)
				c[i][i] -= dynamics.Rate;

			if( !mole_global.list[i]->parentLabel.empty() )
				continue;

			b[i] -= mole.species[i].den*dynamics.Rate;

			if (!mole_global.list[i]->isMonatomic() || mole_global.list[i]->charge < 0 || ! mole_global.list[i]->lgGas_Phase ||
				( mole_global.list[i]->isMonatomic() && !mole_global.list[i]->nAtom.begin()->first->lgMeanAbundance() ) )
			{
				b[i] += dynamics.molecules[i]*dynamics.Rate;
			}
			else if (mole_global.list[i]->charge == 0)
			{
				ASSERT( mole_global.list[i]->isMonatomic() );
				ASSERT( (int)mole_global.list[i]->nAtom.size() == 1 );
				count_ptr<chem_atom> atom = mole_global.list[i]->nAtom.begin()->first;
				long nelem = atom->el->Z-1;
				if( nelem >= LIMELM )						
					continue;
				source = 0.0;
				for ( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
				{
					source += dynamics.Source[nelem][ion] * atom->frac;
					if (unresolved_atom_list[nelem]->ipMl[ion] == -1)
						source -= dense.xIonDense[nelem][ion]*dynamics.Rate * atom->frac;
				}
				b[i] += source;
			}
		}
	}
}