/* 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 */
/*CO_step fills in matrix for heavy elements molecular routines */
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "ionbal.h"
#include "thermal.h"
#include "phycon.h"
#include "hmi.h"
#include "conv.h"
#include "trace.h"
#include "coolheavy.h"
#include "timesc.h"
#include "lapack.h"
#include "mole.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 CO_solve(
	/* set true if we found neg pops */
	int *lgNegPop, 
	/* set true if we tried to compute the pops, but some were zero */
	int *lgZerPop )
{

	int32 ipiv[NUM_COMOLE_CALC], merror;
	long int i, j, k,
		nelem , ion , ion2;
	double
		matrix_value, 
		co_denominator;
	float cartot_mol, oxytot_mol , cartot_ion, oxytot_ion,  siltot_ion,
		 nitrotot_ion,  sulptot_ion, chlotot_ion;
		float abundan;
	
		/* call CO_step(), which will calculate the matrix elements */
		CO_step();

	cartot_mol = dense.xMolecules[ipCARBON] + co.hevmol[ipATC] + co.hevmol[ipCP];
	oxytot_mol = dense.xMolecules[ipOXYGEN] + co.hevmol[ipATO] + co.hevmol[ipOP];

	ASSERT( cartot_mol >= 0. && oxytot_mol >= 0.);


	*lgNegPop = FALSE; 
	*lgZerPop = FALSE;

		/* zero out molecular charge transfer rates */
	for(nelem=ipLITHIUM; nelem < LIMELM; ++nelem)
	{
		for(ion=0; ion < nelem+2; ++ion)
		{
			/*zero out the arrays */				
			mole.sink[nelem][ion]     = 0.;
			mole.source[nelem][ion]   = 0.;

			for(ion2=0; ion2< nelem+2; ++ion2)
			{
				mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
			}
		}
	}


	/* >>chng Oct. 21, 2004 -- Terms that contribute to the ionization balance of C, O, S, Si, and N
	will now be inserted directly into the ion solver.  Before, we just corrected the recombination rate
	by scaling saying that IONIZATION_RATE = RECOMBINATION_RATE * [n(X+)/n(X0)].  This code follows the logic
	of hmole_step, written by Robin Williams. */

	
	/* sink terms should include only terms that either form or destroy an atomic or singly ionized
		element in the network, but not both.  This means that terms that destroy X at the expense of 
		forming X+ can't be included.  The rate of destruction of X or X+ is equal to the formation of 
		the inverse.  We therefore add this rate to sink, which gets rid of this dependence */

	/* The following code is all generalized.  After all the molecules, the total number being NUM_HEAVY_MOLEC,
	   there are 2*NUM_ELEMENTS remaining.  The first set, equal to NUM_ELEMENTS, are the ionized elemental
	   species.  The second set, also equal to NUM_ELEMENTS, are the atomic species.  They are in order such that
		   
	   array number for X = array number for (X+) + NUM_ELEMENTS

	   In the future, if one wants to add another element to the network, such as Na the macros must be changed.
	   Also, the array number for Na and Na+ must obey the equation above.  If this is done, then the sources,
	   sinks, and molecular recombination terms are automatically added!!! */

	/* co.nelem_hevmol[i][j] is just the dominant element for a given molecule.  For an element this
	   string is just the element itself!! */

	for(i=NUM_HEAVY_MOLEC;i<NUM_HEAVY_MOLEC+NUM_ELEMENTS; ++i)
	{ 
		mole.sink[co.nelem_hevmol[i]][0] -= (mole.c[i+NUM_ELEMENTS][i+NUM_ELEMENTS] + 
								mole.c[i+NUM_ELEMENTS][i])*dense.lgElmtOn[co.nelem_hevmol[i]];

		mole.sink[co.nelem_hevmol[i]][1] -= (mole.c[i][i] + mole.c[i][i+NUM_ELEMENTS] )*dense.lgElmtOn[co.nelem_hevmol[i]];
	}

	/* source terms must be multiplied by the density of the corresponding matrix element */

	for(j=NUM_HEAVY_MOLEC;j < NUM_HEAVY_MOLEC+NUM_ELEMENTS; ++j)
	{
		for(i=0; i < NUM_COMOLE_CALC; ++i)
		{
			mole.source[co.nelem_hevmol[j]][1] += co.hevmol[i]*mole.c[i][j]*dense.lgElmtOn[co.nelem_hevmol[j]];
			mole.source[co.nelem_hevmol[j]][0] += co.hevmol[i]*mole.c[i][j+NUM_ELEMENTS]*dense.lgElmtOn[co.nelem_hevmol[j+NUM_ELEMENTS]];
		}
	}

	/* subtract out diagonal terms, they are already in sink.  Also subtract recombination terms,
	   as these are done by mole.xMoleChTrRate */

	for(i=NUM_HEAVY_MOLEC;i < NUM_HEAVY_MOLEC+NUM_ELEMENTS; ++i)
	{ 
		mole.source[co.nelem_hevmol[i]][1]   -= ( mole.c[i][i]*co.hevmol[i] + 
			mole.c[i+NUM_ELEMENTS][i]*co.hevmol[i+NUM_ELEMENTS])*dense.lgElmtOn[co.nelem_hevmol[i]] ;

		mole.source[co.nelem_hevmol[i]][0]   -= ( mole.c[i+NUM_ELEMENTS][i+NUM_ELEMENTS]*co.hevmol[i+NUM_ELEMENTS] + 
			mole.c[i][i+NUM_ELEMENTS]*co.hevmol[i])*dense.lgElmtOn[co.nelem_hevmol[i+NUM_ELEMENTS]] ;

		/* The source terms, as they are right now, include negative contributions.  This is because the 
		linearization "trick" creates source terms.   Take for example the reaction:
		C+        H2O      =>        HCO+   H          
		This reaction destroys C+, but in the act of linearizing, there will be the following term:
		mole.c[ipH2O][ipCP]
		This is only a numerical "trick" and not a real matrix term.  All terms like this
		have to be removed to get the "true" source terms.  Fortunately, the sum of all these terms
		is just mole.b[i].  To remove these terms, we have to add mole.b[i] if the overall contribution from these terms
		was negative, subtract if their contribution was positive.  This is done by subtracting mole.b[i] 
		from the current value of source. */

		mole.source[co.nelem_hevmol[i]][1] = mole.source[co.nelem_hevmol[i]][1] - mole.b[i];
		mole.source[co.nelem_hevmol[i]][0] = mole.source[co.nelem_hevmol[i]][0] - mole.b[i+NUM_ELEMENTS];

	}

	/* negative source terms are actually destruction mechanisms, so should go into the sink vector.
	   source and sinks have different units, so if source is negative divide by the density of
	   the corresponding element to get same units.  For example, if:

	   mole.source[ipCARBON][0]

	   is negative, we divide by dense.xIonDense[ipCARBON][0] to get rate in same units as mole.sink */

	for(i=2; i < LIMELM; ++i)
	{
		for(j=0; j < 2; ++j)
		{
			/*if source term is negative, make it a sink and then set to zero*/ 
			if(mole.source[i][j] < 0)
			{
				mole.sink[i][j] -= (mole.source[i][j]/SDIV(dense.xIonDense[i][j]));
				mole.source[i][j] = 0;
			}
		}
	}
	/* These are rates of recombination for atomic and singly ionized species that are due to molecular processes
	   This will be added to ion_solver to get correct ionization balance */

	for(i=NUM_HEAVY_MOLEC;i < NUM_HEAVY_MOLEC+NUM_ELEMENTS; ++i)
	{ 
		mole.xMoleChTrRate[co.nelem_hevmol[i]][1][0] = (float)(mole.c[i][i+NUM_ELEMENTS] - 
				ionbal.RateRecomTot[co.nelem_hevmol[i]][0])*dense.lgElmtOn[co.nelem_hevmol[i]];

		mole.xMoleChTrRate[co.nelem_hevmol[i]][0][1] =  (float)(mole.c[i+NUM_ELEMENTS][i] - 
				ionbal.RateIonizTot[co.nelem_hevmol[i+NUM_ELEMENTS]][0])*dense.lgElmtOn[co.nelem_hevmol[i+NUM_ELEMENTS]];
	}
	
	/* If rate for going from 1-0 or 0-1 is negative then it should be added to the inverse process */
	for(i=2; i < LIMELM; ++i)
	{
		for(j=0; j < 2; ++j)
		{
			for(k=0; k< 2; ++k)
			{
			/*only possible charge transfers are 0-1 or 1-0 */		
				if(j != k)
				{
					if( mole.xMoleChTrRate[i][j][k] < 0. )
					{
						mole.xMoleChTrRate[i][k][j] -= mole.xMoleChTrRate[i][j][k];
						mole.xMoleChTrRate[i][j][k] = 0.;
					}
				}
			}
		}
	}

	cartot_ion = dense.gas_phase[ipCARBON];
	oxytot_ion = dense.gas_phase[ipOXYGEN];
	siltot_ion = dense.gas_phase[ipSILICON];
	nitrotot_ion = dense.gas_phase[ipNITROGEN];
	sulptot_ion = dense.gas_phase[ipSULPHUR];
	chlotot_ion = dense.gas_phase[ipCHLORINE];


	for ( i=2; i < ipCARBON+2; i++ )
	{
		cartot_ion -= dense.xIonDense[ipCARBON][i];
	}
	for ( i=2; i < ipOXYGEN+2; i++ )
	{
		oxytot_ion -= dense.xIonDense[ipOXYGEN][i];
	}

	for ( i=2; i < ipSILICON+2; i++ )
	{
		siltot_ion -= dense.xIonDense[ipSILICON][i];
	}
	
	for ( i=2; i < ipNITROGEN+2; i++ )
	{
		nitrotot_ion -= dense.xIonDense[ipNITROGEN][i];
	}

	for ( i=2; i < ipSULPHUR+2; i++ )
	{
		sulptot_ion -= dense.xIonDense[ipSULPHUR][i];
	}

	for ( i=2; i < ipCHLORINE+2; i++ )
	{
		chlotot_ion -= dense.xIonDense[ipCHLORINE][i];
	}
	

	/* at this point the cartot_mol, sum of molecules and atom/first ion,
	 * should be equal to the gas_phase minus double and higher ions */
	/*fprintf(ioQQQ," dbuggggas\t %.2f\t%f\t%f\n",fnzone,
		cartot_mol/cartot_ion,
		oxytot_mol/oxytot_ion);*/

	/* these will be used in population equation in case of homogeneous equation */

	mole.b[ipATC]  = 0;
	mole.b[ipATO]  = 0;
	mole.b[ipATSI] = 0;
	mole.b[ipATN]  = 0;
	mole.b[ipATS]  = 0;
	mole.b[ipATCl]  = 0;

	mole.b[ipATC]  = cartot_ion;
	mole.b[ipATO]  = oxytot_ion;
	mole.b[ipATSI] = siltot_ion;
	mole.b[ipATN]  = nitrotot_ion;
	mole.b[ipATS]  = sulptot_ion;
	mole.b[ipATCl] = chlotot_ion;


	/* <<chng 03 Nov 11--Nick Abel,  Set up the non-zero matrix elements that go into the conservation equations
	   for atomic C, O, and Si, this is now set up by looping over the atomic species in
	   co.h and setting the number of atoms of C, O, and Si equal to the variable co.nCarb, 
	   co.nOxyg, and co.nSili respectively.  For every element (excluding Hydrogen) not in
	   the network an if statement will be needed */

	/* loop over last NUM_ELEMENTS elements in co vector - these are atoms of the NUM_ELEMENTS */
	for ( i=NUM_COMOLE_CALC - NUM_ELEMENTS; i < NUM_COMOLE_CALC; i++ )
	{
		for ( j=0; j < NUM_COMOLE_CALC; j++ )

		{
			if ( i == (NUM_COMOLE_CALC - 6) )	
				
			{	
				matrix_value = co.nCarb[j];
				mole.c[j][i] = 0;
				mole.c[j][i] = matrix_value;
			}

			else if ( i == (NUM_COMOLE_CALC - 5 ))	
				
			{
				matrix_value = co.nOxyg[j];
				mole.c[j][i] = 0;
				mole.c[j][i] = matrix_value;
			}
			
			else if ( i == (NUM_COMOLE_CALC - 4) )
				
			{	
				matrix_value = co.nSili[j];
				mole.c[j][i] = 0;
				mole.c[j][i] = matrix_value;
			}

			else if ( i == (NUM_COMOLE_CALC - 3) )
				
			{	
				matrix_value = co.nNitr[j];
				mole.c[j][i] = 0;
				mole.c[j][i] = matrix_value;
			}
		
			else if ( i == (NUM_COMOLE_CALC - 2) )
				
			{	
				matrix_value = co.nSulp[j];
				mole.c[j][i] = 0;
				mole.c[j][i] = matrix_value;
			}

			else if ( i == (NUM_COMOLE_CALC - 1) )
				
			{	
				matrix_value = co.nChlo[j];
				mole.c[j][i] = 0;
				mole.c[j][i] = matrix_value;
			}


		}

	}

	/*--------------------------------------------------------------------
		* */
	/*printf( " Here are all matrix elements\n ");

	for (i=0; i < NUM_COMOLE_CALC; i++)
		
	{
		

			for (j=0; j < NUM_COMOLE_CALC; j++)
			{
				printf( "%4.4s", co.chLab[i] );
				printf( "%4.4s\t", co.chLab[j] );
				printf( "%.8e,\n", mole.c[j][i] );
			}
	}
	printf( "b's are:\n" );
	for (i=0; i < NUM_COMOLE_CALC; i++)
	{
		printf( "%.8e,\n", mole.b[i] );
	}*/


	if ( trace.lgTrace && trace.lgTr_CO_Mole )
	{
		fprintf( ioQQQ, " COMOLE matrix\n           " );

#		if 0
		for ( i=0; i < MIN2(NUM_COMOLE_CALC,8); i++ )
		{
			fprintf( ioQQQ, "%4.4s", co.chLab[i] );
		}
		fprintf( ioQQQ, "     \n" );

		for ( j=0; j < NUM_COMOLE_CALC; j++ )
		{
			fprintf( ioQQQ, " %4.4s", co.chLab[j] );
			fprintf( ioQQQ, " " );
			for ( i=0; i < MIN2(NUM_COMOLE_CALC,8); i++ )
			{
				fprintf( ioQQQ, "%9.1e", mole.c[i][j] );
			}
			fprintf( ioQQQ, "\n" );
		}
#		endif

		fprintf( ioQQQ, " COMOLE matrix\n           " );

		for ( i=0; i < NUM_COMOLE_CALC; i++ )
		{
			fprintf( ioQQQ, "%4.4s", co.chLab[i] );
		}
		fprintf( ioQQQ, "     \n" );

		for ( j=0; j < NUM_COMOLE_CALC; j++ )
		{
			fprintf( ioQQQ, " %4.4s", co.chLab[j] );
			fprintf( ioQQQ, " " );
			for ( i=0; i < NUM_COMOLE_CALC; i++ )
			{
				/*fprintf( ioQQQ, "%9.1e", mole.c[i][j] );*/
				fprintf( ioQQQ, ", %.15e", mole.c[i][j] );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, " COMOLE b\n           " );
		for ( j=0; j < NUM_COMOLE_CALC; j++ )
		{
			fprintf( ioQQQ, " %4.4s", co.chLab[j] );
			fprintf( ioQQQ, ", %.15e\n",mole.b[j] );
		}

		fprintf( ioQQQ, 
			" COMOLE H2 den:%10.3e, H2,3+=%10.2e%10.2e Carb sum=%10.3e Oxy sum=%10.3e\n", 
			hmi.H2_total, 
			hmi.Hmolec[ipMH2p], 
			hmi.Hmolec[ipMH3p], 
			mole.c[NUM_HEAVY_MOLEC][ipATC], 
			mole.c[NUM_HEAVY_MOLEC][ipATO] );

	}

	/* remember longest CO timescale */

	if ( -mole.c[ipCO][ipCO] > SMALLFLOAT )
	{
		/* this is rate CO is destroyed, equal to formation rate in equilibrium */
		timesc.AgeCOMoleDest = -1./mole.c[ipCO][ipCO];
		/* moved to radinc */
		/*timesc.BigCOMoleForm = (float)MAX2(timesc.BigCOMoleForm,timesc.AgeCOMoleForm);*/
	}
	else
	{
		timesc.AgeCOMoleDest = 0.;
	}

	/* this one may be more robust */
	for ( j=0; j < NUM_COMOLE_CALC; j++ )
	{
		for ( i=0; i < NUM_COMOLE_CALC; i++ )
		{
			mole.amat[i][j] = mole.c[i][j];
		}
	}

	merror = 0;

	/*lint -e740 unusual pointer case */
	getrf_wrapper(NUM_COMOLE_CALC,NUM_COMOLE_CALC,(double*)mole.amat,NUM_COMOLE_CALC, ipiv,&merror);
	/*DGETRF(NUM_COMOLE_CALC,NUM_COMOLE_CALC,(double*)amat,NUM_COMOLE_CALC, ipiv,&merror);*/

	if ( merror != 0 )
	{
		fprintf( ioQQQ, " CO_solve getrf_wrapper error\n" );
		puts( "[Stop in CO_solve]" );
		cdEXIT(EXIT_FAILURE);
	}

	getrs_wrapper('N',NUM_COMOLE_CALC,1,(double*)mole.amat,NUM_COMOLE_CALC,ipiv,mole.b,NUM_COMOLE_CALC,&merror);
	/*lint +e740 unusual pointer case */

	if ( merror != 0 )
	{
		fprintf( ioQQQ, " CO_solve: dgetrs finds singular or ill-conditioned matrix\n" );
		puts( "[Stop in CO_solve]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* check for negative populations, which happens when 100% co */
	*lgNegPop = FALSE;
	*lgZerPop = FALSE;
	co_denominator = 0;
#	if 0
	if(fnzone > 1.02)
	{
			for ( i=0; i < NUM_COMOLE_CALC; i++)
			{			

				for ( j=61; j < 62; j++  )
				{
					
					printf("ZONE\t%.5e\tSPECIES_1\t%s\tSPECIES_2\t%s\tDEST_RATE\t%.3e\tbvec\t%.3e\n",
						fnzone, co.chLab[i],co.chLab[j],mole.c[i][j],mole.b[i]);

				}
			}
	}
#	endif
	for ( i=0; i < NUM_COMOLE_CALC; i++ )
	{			

		if ( mole.b[i] < 0. )
		{

			/* >>chng 03 sep 11, the following*/
			/* comparison between the solution vector from 32 bit machines vs
			 * maple on a pc shows that very small negative numbers are produced by
			 * the linear algebra package used by the code, but is positive with
			 * the math package.  allow very small negative numbers to occur,
			 * and silently reset to a postive number. 
			 *
			 * >>chng 04 feb 25
			 * Here we want to check if the negative abundance is a significant
			 * fraction of any of the species in the network.  The code below
			 * checks to see which elements the molecule in question is made of,
			 * then finds which one has the least abundance.  The one with
			 * the least abundance will then be used as the divisor in checking 
			 * if the negative abundance is too large to be ignored.*/
			
			/*  >>chng 04 apr 21, when an element in the chemical network is 
			 *  turned off, the molecular abundance of a species containing that
			 *  element was not zero but rather a small number of order 1e-20. When 
			 *  these abundances went negative, the code thought the abundances were
			 *  important because it checks the ratio of the species to its gas phase
			 *  abundance.  When the gas phase abundance is zero, the denominator is set 
			 *  to SMALLFLOAT, which is ~1e-35.  This led to a ratio of molecule to gas
			 *	phase of ~1e15, clearly unphysical.  Here the value of co_denominator 
			 *  will be set to a high value if the gas phase abundance is zero.*/

			if( dense.lgElmtOn[co.nelem_hevmol[i]] )
			{
				co_denominator = dense.gas_phase[co.nelem_hevmol[i]];
			}
			else
			{
				/* >>chng 04 apr 20, set to zero if element is off */
				co_denominator = 1e10;
			}

			
			/* we must have a positive value by now, unlesss element is turned off */
			ASSERT( co_denominator > SMALLFLOAT || !dense.lgElmtOn[co.nelem_hevmol[i]] );

			/* >>chng 04 feb 28, threshold from -1e-10 to -1e-6, the 
			 * level of roundoff in a float */
			/*if( mole.b[i] / MAX2 (co_denominator, SMALLFLOAT) > -1e-10) */
			/*if( mole.b[i] / SDIV(co_denominator) > -1e-6 ) */
			/* >>chng 04 oct 31, change limit from -1e-6 to -1e-5, produced only
			 * a few comments in map.in */
			if( mole.b[i] / SDIV(co_denominator) > -1e-5 ) 
			{
				/* this case, it is only slightly negative relative to
				 * the parent species - multiply by -1 to make positive
				 * and press on */
				mole.b[i] *= -1.;
			}
			else
			{
				/*>>chng 04 mar 06 press on */
				*lgNegPop = TRUE;
				/* 05 jul 16, there had been a bug such that CO was always evaluated as long as
				 * the temperature was below 20,000K.  Once fixed, one sim turned on CO mid way
				 * into the H+ zone and had neg pops during first trys - change check on whether it
				 * is to be commented on only if this is not first zone with CO soln */
				/*if( nzone>0 )*/
				if( nzone>co.co_nzone )
				{
					static long int nzFail=-2;
					long int nzFail2 = (long)fnzone*100;
					/* during a map we may be in search phase but not in first zone */
					if( !conv.lgSearch )
						fprintf(ioQQQ," PROBLEM ");
					fprintf(ioQQQ,
						" CO_solve neg pop for species %li %s, value is %.2e rel value is %.2e zone %.2f Te %.4e Search?%c\n",
						i , 
						co.chLab[i] , 
						/* abs value */
						mole.b[i], 
						/* relative value */
						mole.b[i] / SDIV(co_denominator) ,
						fnzone,
						phycon.te,TorF( conv.lgSearch ) );
					/* if well beyond search phase and still having problems, announce them 
					 * only announce one failure per sweep through solvers by ConvFail */
					if( nzFail2 !=nzFail )
					{
						nzFail = nzFail2;
						ConvFail( "pops" , "CO");
					}
				}
				mole.b[i] *= -1;
				/*>>chng 04 jan 24 set to zero instead of *-1,
				 * h2_cr.in co density went to infinite, with some negative
				 * var getting to - inf 
				 *>>chng 04 nov 03, remove this - not needed h2_cr works without
				mole.b[i] = SMALLFLOAT;*/
			} 
		}
		else if ( mole.b[i] == 0. )
		{
			/* this is not used for anything in calling routine and 
			 * could be cleaned up */
			*lgZerPop = TRUE;
			/* >>chng 04 feb 28, zero pop not really a problem in very deep neutral
			 * gas - this happens */
#			if 0
			if( /*nzone>0 ||*/ CODEBUG>1 )
			{
				fprintf(ioQQQ," PROBLEM ");
				fprintf(ioQQQ,
					" CO_solve zero pop for species %li %s, value is %.2e zone %li\n",
					i , co.chLab[i] , mole.b[i], nzone);
			}
#			endif
			mole.b[i] = SMALLFLOAT;
		}
		co.hevmol[i] = (float)mole.b[i];
	}
	/* >>chng 04 mar 06 pass negative pop as a problem, but not a fatal one,
	 * also do not call this during 0th zone when search for conditions
	 * is underway */
	/* 05 jul 16, there had been a bug such that CO was always evaluated as long as
	 * the temperature was below 20,000K.  Once fixed, on sim turned on CO mid way
	 * into the H+ zone and had neg pops during first trys - change check on whether it
	 * is to be commented on only if this is not first zone with CO soln */
	/*if( *lgNegPop && (nzone>0 &&!conv.lgSearch)  )*/
	if( *lgNegPop && (nzone>co.co_nzone &&!conv.lgSearch)  )
	{
		conv.lgConvPops = FALSE;
		fprintf(ioQQQ," CO network negative population occurred, Te=%.4e, calling ConvFail. ",
			phycon.te);
		fprintf( ioQQQ, " CO/C=%9.1e", co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]) );
		fprintf( ioQQQ, "\n" );
		/*ConvFail("pops", "CO");*/
		*lgNegPop = FALSE;
	}
	/* if negative pops were present need to renorm b to get
	 * proper sum rule */
	if( 0 && *lgNegPop )
	{

		for(j=0; j<2; ++j )
		{
			double sumcar = 0.;
			double sumoxy = 0.;
			double renorm;
			for ( i=0; i < NUM_COMOLE_CALC; i++ )
			{
				/* this case, different molecules */
				sumcar += co.hevmol[i]*co.nCarb[i];
				sumoxy += co.hevmol[i]*co.nOxyg[i];
			}
			renorm = (cartot_ion + oxytot_ion) / SDIV( sumcar+sumoxy);
			if(j)
				fprintf(ioQQQ,"\t%f\n", renorm);
			else
				fprintf(ioQQQ,"renormco\t%.2f\t%f", fnzone, renorm);
			for ( i=0; i < NUM_COMOLE_CALC; i++ )
			{
				co.hevmol[i] *= (float)renorm;
			}
		}
	}

	if ( merror != 0 )
	{
		fprintf( ioQQQ, " COMOLE matrix inversion error, MERROR=%5ld zone=%5ld\n", 
				 (long)merror, nzone );
		ShowMe();
		fprintf( ioQQQ, " Product matrix\n           " );

		for ( i=0; i < MIN2(NUM_COMOLE_CALC,8); i++ )
		{
			fprintf( ioQQQ, "%4.4s", co.chLab[i] );
		}
		fprintf( ioQQQ, "     \n" );

		for ( j=0; j < NUM_COMOLE_CALC; j++ )
		{
			fprintf( ioQQQ, " %4.4s", co.chLab[j] );
			fprintf( ioQQQ, " " );

			for ( i=0; i < MIN2(NUM_COMOLE_CALC,8); i++ )
			{
				fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
						 co.hevmol[i] );
			}
			fprintf( ioQQQ, "\n" );
		}

		if ( NUM_COMOLE_CALC >= 9 )
		{
			fprintf( ioQQQ, " COMOLE matrix\n           " );
			for ( i=8; i < NUM_COMOLE_CALC; i++ )
			{
				fprintf( ioQQQ, "%4.4s", co.chLab[i] );
			}
			fprintf( ioQQQ, "     \n" );

			for ( j=0; j < NUM_COMOLE_CALC; j++ )
			{
				fprintf( ioQQQ, " %4.4s", co.chLab[j] );
				fprintf( ioQQQ, " " );
				for ( i=8; i < NUM_COMOLE_CALC; i++ )
				{
					fprintf( ioQQQ, "%9.1e", 
							 mole.amat[i][j]* co.hevmol[i] );
				}
				fprintf( ioQQQ, "\n" );
			}
		}

		fprintf( ioQQQ, " Mole dens:" );
		for ( j=0; j < NUM_COMOLE_CALC; j++ )
		{
			fprintf( ioQQQ, " %4.4s:%.2e", co.chLab[j]
					 , co.hevmol[j] );
		}
		fprintf( ioQQQ, " \n" );

		puts( "[Stop in CO_solve]" );
		cdEXIT(EXIT_FAILURE);
	}

	if ( trace.lgTrace && trace.lgTr_CO_Mole )
	{
		fprintf( ioQQQ, " Product matrix\n           " );

		for ( i=0; i < MIN2(NUM_COMOLE_CALC,8); i++ )
		{
			fprintf( ioQQQ, "%4.4s", co.chLab[i] );
		}
		fprintf( ioQQQ, "     \n" );

		for ( j=0; j < NUM_COMOLE_CALC; j++ )
		{
			fprintf( ioQQQ, " %4.4s", co.chLab[j] );
			fprintf( ioQQQ, " " );
			for ( i=0; i < MIN2(NUM_COMOLE_CALC,8); i++ )
			{
				fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
						 co.hevmol[i] );
			}
			fprintf( ioQQQ, "\n" );

		}

		if ( NUM_COMOLE_CALC >= 9 )
		{
			fprintf( ioQQQ, " COMOLE matrix\n           " );
			for ( i=8; i < NUM_COMOLE_CALC; i++ )
			{
				fprintf( ioQQQ, "%4.4s", co.chLab[i] );
			}
			fprintf( ioQQQ, "     \n" );

			for ( j=0; j < NUM_COMOLE_CALC; j++ )
			{
				fprintf( ioQQQ, " %4.4s", co.chLab[j] );
				fprintf( ioQQQ, " " );

				for ( i=8; i < NUM_COMOLE_CALC; i++ )
				{
					fprintf( ioQQQ, "%9.1e", mole.amat[i][j]* co.hevmol[i] );
				}
				fprintf( ioQQQ, "\n" );
			}
		}

		fprintf( ioQQQ, " Mole dens:" );
		for ( j=0; j < NUM_COMOLE_CALC; j++ )
		{
			fprintf( ioQQQ, " %4.4s:%.2e", co.chLab[j]
					 , co.hevmol[j] );
		}
		fprintf( ioQQQ, " \n" );
	}

	/* heating due to CO photodissociation */
	co.CODissHeat = (float)(co.CO_photo_dissoc_rate*co.hevmol[ipCO]*1e-12);

	thermal.heating[0][9] = co.CODissHeat;

	/* the real multi-level model molecule */
	abundan = co.hevmol[ipCO];
	/* IonStg and nelem were set to 2, 0 in makelevlines */
	ASSERT( C12O16Rotate[0].IonStg < LIMELM+2 );
	ASSERT( C12O16Rotate[0].nelem-1 < LIMELM+2 );
	dense.xIonDense[ C12O16Rotate[0].nelem-1][C12O16Rotate[0].IonStg-1] = abundan;

	CO_PopsEmisCool(&C12O16Rotate, nCORotate,abundan , "12CO",
		&CoolHeavy.C12O16Rot,&CoolHeavy.dC12O16Rot );
	/*if( nzone>400 )
		fprintf(ioQQQ,"DEBUG co cool\t%.2f\t%.4e\t%li\t%.4e\t%.4e\n",
			fnzone,CoolHeavy.C12O16Rot,nCORotate,abundan,phycon.te );*/

	abundan = co.hevmol[ipCO]/co.RatioC12O16_2_C13O16;
	/* IonStg and nelem were set to 3, 0 in makelevlines */
	ASSERT( C13O16Rotate[0].IonStg < LIMELM+2 );
	ASSERT( C13O16Rotate[0].nelem-1 < LIMELM+2 );
	dense.xIonDense[ C13O16Rotate[0].nelem-1][C13O16Rotate[0].IonStg-1] = abundan;

	CO_PopsEmisCool(&C13O16Rotate, nCORotate,abundan ,"13CO",
		&CoolHeavy.C13O16Rot,&CoolHeavy.dC13O16Rot );

	/* now set total density of each element locked in gas phase */
	for ( i=0;i<LIMELM; ++i )
	{
		dense.xMolecules[i] = 0.;
	}

	/* total number of H per unit vol in molecules,
	 * of course not including H0/H+ */
	for (i=0;i<N_H_MOLEC;i++) 
	{
		dense.xMolecules[ipHYDROGEN] += hmi.Hmolec[i]*hmi.nProton[i];
	}
	dense.xMolecules[ipHYDROGEN] -= (hmi.Hmolec[ipMH] + hmi.Hmolec[ipMHp]);

	/* >>chng 02 sep 05, dense.xMolecules[[ipHYDROGEN] is the part of H 
	 * that is not computed in hmole
	 * add in gas phase abundances locked in molecules
	 * H is special since treated in bidiag with all hmole molecules
	 * xMolecules are the densities of these species that are done in co,
	 * this sum is only up to NUM_HEAVY_MOLEC and so does not include the atoms/ions */
	for(i=0; i<NUM_HEAVY_MOLEC; ++i)
	{
		/* these were all initialized to zero above */
		dense.xMolecules[ipHYDROGEN] +=co.nHydr[i]*co.hevmol[i];
		dense.xMolecules[ipCARBON]   +=co.nCarb[i]*co.hevmol[i];
		dense.xMolecules[ipOXYGEN]   +=co.nOxyg[i]*co.hevmol[i];
		dense.xMolecules[ipSILICON]  +=co.nSili[i]*co.hevmol[i];
		dense.xMolecules[ipSULPHUR]  +=co.nSulp[i]*co.hevmol[i];
		dense.xMolecules[ipNITROGEN] +=co.nNitr[i]*co.hevmol[i];
		dense.xMolecules[ipCHLORINE] +=co.nChlo[i]*co.hevmol[i];
	}

	/* This takes the average of an element's atomic and singly ionized density from ion_solver 
	   and comole and sets the solution to the updated co.hevmol[i] value.  The macro NUM_HEAVY_MOLEC
	   is the number of heavy element molecules in the network.  In addition there are 2*NUM_ELEMENTS
	   in the network, half atomic elements and half ionized elements.  Finally, the total number of species
	   in the network, including elements, equals NUM_COMOLE_CALC.  The following will first take the average of 
	   the ionized species (between NUM_HEAVY_MOLEC and NUM_HEAVY_MOLEC+NUM_ELEMENTS) and the other statement will
	   take the average of the atomic species (between NUM_HEAVY_MOLEC+NUM_ELEMENTS and NUM_COMOLE_CALC).  
	   This is generalized, so that if one wants to insert, for example, a sodium chemistry, once the macro's
	   are updated this if statement immediately works */

	for(i=NUM_HEAVY_MOLEC;i<NUM_HEAVY_MOLEC+NUM_ELEMENTS; ++i)
	{ 
		/*printf("MOL\t%3e\tION\t%e\tSPECIES\t%s\n", co.hevmol[i],dense.xIonDense[co.nelem_hevmol[i]][1],co.chLab[i]);*/
		co.hevmol[i] = ((dense.xIonDense[co.nelem_hevmol[i]][1]+ co.hevmol[i])/2)*dense.lgElmtOn[co.nelem_hevmol[i]];
	}

	for(i=NUM_HEAVY_MOLEC+NUM_ELEMENTS;i<NUM_COMOLE_CALC; ++i)
	{ 
		/*printf("MOL\t%3e\tION\t%e\tSPECIES\t%s\n", co.hevmol[i],dense.xIonDense[co.nelem_hevmol[i]][0],co.chLab[i]);*/
		co.hevmol[i] = ((dense.xIonDense[co.nelem_hevmol[i]][0]+ co.hevmol[i])/2)*dense.lgElmtOn[co.nelem_hevmol[i]];
	}


	/* check whether ion and chem solvers agree yet */
	if( dense.lgElmtOn[ipCARBON] &&
		fabs(dense.xIonDense[ipCARBON][0]- co.hevmol[ipATC])/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "CO C0 con");
		conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][0];
		conv.BadConvIoniz[1] = co.hevmol[ipATC];
	}
	else if( dense.lgElmtOn[ipCARBON] &&
		fabs(dense.xIonDense[ipCARBON][1]- co.hevmol[ipCP])/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "CO C1 con");
		conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][1];
		conv.BadConvIoniz[1] = co.hevmol[ipCP];
	}
	else if( dense.lgElmtOn[ipOXYGEN] &&
		fabs(dense.xIonDense[ipOXYGEN][0]- co.hevmol[ipATO])/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "CO O0 con");
		conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][0];
		conv.BadConvIoniz[1] = co.hevmol[ipATO];
	}
	else if( dense.lgElmtOn[ipOXYGEN] &&
		fabs(dense.xIonDense[ipOXYGEN][1]- co.hevmol[ipOP])/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "CO O1 con");
		conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][1];
		conv.BadConvIoniz[1] = co.hevmol[ipOP];
	}

	/* now update ionization distribution of the elements we just did */
	dense.xIonDense[ipCARBON][0] = co.hevmol[ipATC]*dense.lgElmtOn[ipCARBON] ; 
	dense.xIonDense[ipCARBON][1] = co.hevmol[ipCP]*dense.lgElmtOn[ipCARBON]; 
	dense.IonLow[ipCARBON] = 0;
	dense.IonHigh[ipCARBON] = MAX2( dense.IonHigh[ipCARBON] , 1 );

	dense.xIonDense[ipOXYGEN][0] = co.hevmol[ipATO]*dense.lgElmtOn[ipOXYGEN] ; 
	dense.xIonDense[ipOXYGEN][1] = co.hevmol[ipOP]*dense.lgElmtOn[ipOXYGEN]; 
	dense.IonLow[ipOXYGEN] = 0;
	dense.IonHigh[ipOXYGEN] = MAX2( dense.IonHigh[ipOXYGEN] , 1 );

	dense.xIonDense[ipSILICON][0] = co.hevmol[ipATSI]*dense.lgElmtOn[ipSILICON] ; 
	dense.xIonDense[ipSILICON][1] = co.hevmol[ipSIP]*dense.lgElmtOn[ipSILICON]; 
	dense.IonLow[ipSILICON] = 0;
	dense.IonHigh[ipSILICON] = MAX2( dense.IonHigh[ipSILICON] , 1 );

	dense.xIonDense[ipSULPHUR][0] = co.hevmol[ipATS]*dense.lgElmtOn[ipSULPHUR] ; 
	dense.xIonDense[ipSULPHUR][1] = co.hevmol[ipSP]*dense.lgElmtOn[ipSULPHUR]; 
	dense.IonLow[ipSULPHUR] = 0;
	dense.IonHigh[ipSULPHUR] = MAX2( dense.IonHigh[ipSULPHUR] , 1 );

	dense.xIonDense[ipNITROGEN][0] = co.hevmol[ipATN]*dense.lgElmtOn[ipNITROGEN] ; 
	dense.xIonDense[ipNITROGEN][1] = co.hevmol[ipNP]*dense.lgElmtOn[ipNITROGEN]; 
	dense.IonLow[ipNITROGEN] = 0;
	dense.IonHigh[ipNITROGEN] = MAX2( dense.IonHigh[ipNITROGEN] , 1 );
	
	dense.xIonDense[ipCHLORINE][0] = co.hevmol[ipATCl]*dense.lgElmtOn[ipCHLORINE] ; 
	dense.xIonDense[ipCHLORINE][1] = co.hevmol[ipClP]*dense.lgElmtOn[ipCHLORINE]; 
	dense.IonLow[ipCHLORINE] = 0;
	dense.IonHigh[ipCHLORINE] = MAX2( dense.IonHigh[ipCHLORINE] , 1 );


	/* if populations not conserved then not converged */
#	define EPS_MOLE	0.1
	if( dense.lgElmtOn[ipHYDROGEN] &&
		dense.xMolecules[ipHYDROGEN] > dense.gas_phase[ipHYDROGEN]*(1.+EPS_MOLE) )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "COcon%2i",ipHYDROGEN );
		conv.BadConvIoniz[0] = dense.xMolecules[ipHYDROGEN];
		conv.BadConvIoniz[1] = dense.gas_phase[ipHYDROGEN];
	}
	else if( dense.lgElmtOn[ipCARBON] &&
		dense.xMolecules[ipCARBON] > dense.gas_phase[ipCARBON]*(1.+EPS_MOLE) )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "COcon%2i",ipCARBON );
		conv.BadConvIoniz[0] = dense.xMolecules[ipCARBON];
		conv.BadConvIoniz[1] = dense.gas_phase[ipCARBON];
	}
	else if( dense.lgElmtOn[ipOXYGEN] &&
		dense.xMolecules[ipOXYGEN] > dense.gas_phase[ipOXYGEN]*(1.+EPS_MOLE) )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "COcon%2i",ipOXYGEN );
		conv.BadConvIoniz[0] = dense.xMolecules[ipOXYGEN];
		conv.BadConvIoniz[1] = dense.gas_phase[ipOXYGEN];
	}
	else if( dense.lgElmtOn[ipSILICON] &&
		dense.xMolecules[ipSILICON] > dense.gas_phase[ipSILICON]*(1.+EPS_MOLE) )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "COcon%2i",ipSILICON );
		conv.BadConvIoniz[0] = dense.xMolecules[ipSILICON];
		conv.BadConvIoniz[1] = dense.gas_phase[ipSILICON];
	}
	else if( dense.lgElmtOn[ipSULPHUR] &&
		dense.xMolecules[ipSULPHUR] > dense.gas_phase[ipSULPHUR]*(1.+EPS_MOLE) )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "COcon%2i",ipSULPHUR );
		conv.BadConvIoniz[0] = dense.xMolecules[ipSULPHUR];
		conv.BadConvIoniz[1] = dense.gas_phase[ipSULPHUR];
	}
	else if( dense.lgElmtOn[ipNITROGEN] &&
		dense.xMolecules[ipNITROGEN] > dense.gas_phase[ipNITROGEN]*(1.+EPS_MOLE) )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "COcon%2i",ipNITROGEN );
		conv.BadConvIoniz[0] = dense.xMolecules[ipNITROGEN];
		conv.BadConvIoniz[1] = dense.gas_phase[ipNITROGEN];
	}

	else if( dense.lgElmtOn[ipCHLORINE] &&
		dense.xMolecules[ipCHLORINE] > dense.gas_phase[ipCHLORINE]*(1.+EPS_MOLE) )
	{
		conv.lgConvIoniz = FALSE;
		sprintf( conv.chConvIoniz, "COcon%2i",ipCHLORINE );
		conv.BadConvIoniz[0] = dense.xMolecules[ipCHLORINE];
		conv.BadConvIoniz[1] = dense.gas_phase[ipCHLORINE];
	}
#	undef EPS_MOLE

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


/*lint +e550 */


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

