/* 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 */
/*lgMolecAver average old and new molecular equilibrium balance from CO_solve */
/*CO_drive - public routine, calls CO_solve to converge molecules */
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "hmi.h"
#include "conv.h"
#include "phycon.h"
#include "trace.h"
#include "thermal.h"
#include "called.h"
#include "map.h"
#include "coolheavy.h"
#include "mole.h"

/* says whether the co network is currently set to zero */
static int lgMoleZeroed=TRUE;

/* the limit to H2/H where we will solve for CO */
static double h2lim;

/* this save first valid solution from previous iteration */
static float hev_reinit[NUM_COMOLE_CALC];

/* this is the relative change in OH, which is used as a convergence
 * criteria in lgMolecAver */
#define	COTOLER_MOLAV	0.10

/* flag to control debug statement, if 0 then none, just loops when 1, more when 2 */
#define	CODEBUG	0
/* this is the maximum number of loops CO_drive will go over while
 * trying to converge the molecular network */
#define	LUPMAX_CODRIV	100

/*lgMolecAver average old and new molecular equilibrium balance from CO_solve,
 * returns TRUE if converged */
static int lgMolecAver(
	/* job == "SETUP", set things up during first call,
	 * job == "AVER" take average of arrays */
	const char chJob[10] ,
	char *chReason )
{
	long int i;
	int lgReturn;
	/* this will be used to rescale old saved abundances in constant pressure model */
	static float hden_save_prev_call;
	double conv_fac;

#	ifdef DEBUG_FUN
	fputs( "<+>lgMolecAver()\n", debug_fp );
#	endif
#	define	DEBUG_MOLECAVER	FALSE

	/* this will become reason not converged */
	strcpy( chReason , "none given" );

	/* when called with SETUP, just about to enter solver loop */
	if( strcmp( "SETUP" , chJob ) == 0 )
	{
		static float hden_save_prev_iter;

		/* >>chng 04 apr 29, use pdr solution, scaled to density of C0, as first guess*/
		/* CO_Init set this to -2 when code initialized, so negative
		 * number shows very first call in this model */
		/* >>chng 05 jul 15, do not init if we have never evaluated CO in this iteration
		 * and we are below limit where it should be evaluated */
		if( iteration!=co.iteration_co && 
			hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim )
		{

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

			lgReturn = TRUE;
			return lgReturn;
		}

		else if( co.iteration_co < 0 || lgMoleZeroed )
		{
			/* these are relative to C0 and are the solution to the first zone of the 
			 * TH85 pdr */
			float pdr_mole_co[NUM_COMOLE_CALC]= {
				1.10E-06f, 6.99E-05f, 2.15E-03f, 2.45E-04f, 1.91E-07f,
				4.05E-05f, 2.16E-06f, 1.60E-11f, 1.27E-10f, 1.17E-06f,
				3.44E-17f, 3.71E-09f, 8.59E-13f, 4.01E-10f, 7.02E-16f,
				2.07E-13f, 1.80E-14f, 5.39E-15f, 3.89E-14f, 2.27E-08f,
				4.59E-19f, 9.12E-28f, 1.03E-28f, 8.16E-28f, 2.46E-17f,
				2.22E-17f, 5.99E-12f, 1.31E-11f, 1.21E-11f, 1.00E-13f,
				3.86E-11f, 1.65E-22f, 1.30E-09f, 2.20E-10f, 6.78E-20f,
				1.71E-16f, 1.25E-29f, 3.03E-23f, 1.15E-33f, 3.38E-11f,
				5.07E-11f, 1.07E-17f, 9.89E-17f, 9.09E-21f, 1.91E-18f,
				1.71E-14f, 4.83E-13f, 6.69E-18f, 4.12E-11f, 2.67E-26f,
				2.41E-23f, 3.12E-11f, 6.81E-13f, 4.00E-15f, 1.20E-07f,
				7.03E-12f, 7.60E-14f, 1.05E-11f, 8.91E-17f, 8.83E-24f,
				1.72E-21f, 3.71E-09f, 8.59E-13f, 4.00E-15f, 4.00E-15f,
				4.00E-15f, 4.00E-15f, 4.00E-15f, 4.00E-15f, 4.00E-15f, 
				4.00E-15f, 4.00E-15f, 4.00E-15f, 4.00E-15f, 4.00E-15f, 
				4.00E-15f, 4.00E-15f, 4.00E-15f, 4.00E-15f, 1.00e-15f, 
				4.00E-15f, 4.00E-15f, 1.00e-15f, 1.00e-15f, 1.00e-15f, 
				6.79E+04f, 1.58E+01f, 1.79E+02f, 2.62E-10f, 1.79E+03f, 
				3.71E-09f, 1.00E+00f, 1.13E+05f, 3.56E-04f, 2.25E+00f, 
				6.90E-03f, 6.90E-03f  };

			/* very first attempt at ever obtaining a solution -
			 * called one time per model since co.iteration_co set negative
			 * when cdInit called */
	
			/* >>chng 05 jun 24, during map phase do not reset molecules to zero
			 * since we probably have a better estimate right now */
			if( !map.lgMapBeingDone || lgMoleZeroed )
			{
				for ( i=0; i < NUM_COMOLE_CALC; i++ )
				{
					co.hevmol[i] = dense.xIonDense[ipCARBON][0]*pdr_mole_co[i];
					co.co_save[i] = co.hevmol[i];
				}
			}

			/* we should have a C0 soln at this point */
			ASSERT( dense.xIonDense[ipCARBON][0]>SMALLFLOAT );

			/* make sure order of molecules has not changed */
			ASSERT( 
					ipCH == 0
				&&	ipCHP	== 1
				&&	ipOH==	 2
				&&	ipOHP==	 3
				&&	ipO2  ==	 4
				&&	ipCO==	 5
				&&	ipCOP==	 6
				&&	ipH2O	== 7
				&&	ipH2OP	 ==8
				&&	ipO2P	 ==9
				&&	ipH3OP	 ==10
				&&	ipCH2P	== 11
				&& ipCH2  ==  12
				&& ipHCOP ==  13
				&& ipCH3P ==  14
				&& ipSIH2P == 15
				&& ipSIH   == 16
				&& ipHOSIP == 17
				&& ipSIO  ==  18
				&& ipSIOP ==  19
				&& ipCH3  ==  20
				&& ipCH4  ==  21
				&& ipCH4P ==  22
				&& ipCH5P ==  23
				&&	ipN2	== 24
				&&	ipN2P==	 25
				&&	ipNO	== 26
				&&	ipNOP==	 27
				&&	ipS2  == 28
				&&	ipS2P	== 29
				&&	ipOCN	== 30 );
			ASSERT( 
				  	ipOCNP==	 31
				&&	ipNH 	== 32
				&&	ipNHP 	 ==33
				&&	ipNH2 ==	 34
				&&	ipNH2P ==	 35
				&& ipNH3   == 36
				&& ipNH3P  == 37
				&& ipNH4P ==  38
				&& ipCN    == 39
				&& ipCNP  ==  40
				&& ipHCN  ==  41
				&& ipHCNP ==  42
				&& ipHNO  ==  43
				&& ipHNOP ==  44
				&& ipHS   ==  45
				&& ipHSP  ==  46
				&& ipCS   ==  47
				&&	ipCSP	== 48
				&&	ipNO2==	 49
				&&	ipNO2P==	 50
				&&	ipNS==	 51
				&&	ipNSP==  	 52
				&&	ipSO	== 53
				&&	ipSOP	== 54
				&&	ipSIN	== 55
				&&	ipSINP ==	 56
				&&	ipN2O==	 57
				&&	ipHCSP==	 58
				&&	ipOCS==	 59
				&&	ipOCSP==	 60
				&&	ipC2==	 61
				&&	ipC2P==	 62
				&& ipCCl  ==  63         
				&& ipClO   == 64      
				&& ipHClP  == 65  
				&& ipHCl   == 66  );
			ASSERT( 
				   ipH2ClP == 67 
				&& ipCClP  == 68
				&& ipH2CClP== 69
				&& ipClOP  == 70
				&& ipHNC == 71
				&& ipHCNHP  == 72
				&& ipC2H  == 73
				&& ipC2HP  == 74
				&& ipC2H2  == 75
				&& ipC2H2P  == 76
				&& ipC3H  == 77
				&& ipC3HP  == 78
				&& ipC2H3P  == 79
				&& ipC3  == 80
				&& ipC3P  == 81
				&& ipCOGR   == 82
				&& ipH2OGR   == 83
				&& ipOHGR   == 84
				&& ipCP    ==  85
				&& ipOP    ==  86
				&& ipSIP   ==  87
				&& ipNP    ==  88
				&& ipSP    ==  89
				&& ipClP   ==  90
				&& ipATC   ==  91
				&& ipATO   ==  92
				&& ipATSI  ==  93
				&& ipATN   ==  94
				&& ipATS   ==  95
				&& ipATCl  ==  96
				&&  NUM_COMOLE_CALC==97 );

			/* first guess of these come from ionization solvers */
			co.hevmol[ipATC] = dense.xIonDense[ipCARBON][0];
			co.hevmol[ipCP] = dense.xIonDense[ipCARBON][1];
			co.hevmol[ipATO] = dense.xIonDense[ipOXYGEN][0];
			co.hevmol[ipOP] = dense.xIonDense[ipOXYGEN][1];
			co.hevmol[ipATSI] = dense.xIonDense[ipSILICON][0];
			co.hevmol[ipSIP] = dense.xIonDense[ipSILICON][1];
			co.hevmol[ipATS] = dense.xIonDense[ipSULPHUR][0];
			co.hevmol[ipSP] = dense.xIonDense[ipSULPHUR][1];
			co.hevmol[ipATN] = dense.xIonDense[ipNITROGEN][0];
			co.hevmol[ipNP] = dense.xIonDense[ipNITROGEN][1];
			co.hevmol[ipATCl] = dense.xIonDense[ipCHLORINE][0];
			co.hevmol[ipClP] = dense.xIonDense[ipCHLORINE][1];
		
			/* set iteration flag */
			co.iteration_co = iteration;
			co.co_nzone = nzone;

			/* for constant pressure models when molecules are reset on second
			 * and higher iterations, total density will be different, so
			 * must take this into account when abundances are reset */
			hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
			hden_save_prev_call = dense.gas_phase[ipHYDROGEN];

			if( DEBUG_MOLECAVER )
				fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li zeroing since very first call. H2/H=%.2e\n", 
				iteration,nzone,hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
		}
		else if( iteration > co.iteration_co )
		{
			float ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_iter;

			/* this is first call on new iteration, reset
			 * to first initial abundances from previous iteration */
			for ( i=0; i < NUM_COMOLE_CALC; i++ )
			{
				co.hevmol[i] = hev_reinit[i]*ratio;
				co.co_save[i] = hev_reinit[i]*ratio;
			}

			co.iteration_co = iteration;
			co.co_nzone = nzone;

			if( DEBUG_MOLECAVER )
				fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li resetting since first call on new iteration. H2/H=%.2e\n", 
				iteration,
				nzone,
				hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
		}
		else if( iteration == co.iteration_co && nzone==co.co_nzone+1 )
		{
			/* this branch, second zone with solution, so save previous
			 * zone's solution to reset things in next iteration */
			for ( i=0; i < NUM_COMOLE_CALC; i++ )
			{
				hev_reinit[i] = co.hevmol[i];
			}

			co.co_nzone = -2;
			hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
			hden_save_prev_call = dense.gas_phase[ipHYDROGEN];

			if( DEBUG_MOLECAVER )
				fprintf(ioQQQ,"DEBUG lgMolecAver iteration %li zone %li second zone on new iteration, saving reset.\n", iteration,nzone);
		}
		
		/* didn't do anything, but say converged */
		lgReturn = TRUE;
	}

	else if( strcmp( "AVER" , chJob ) == 0 )
	{
		/* >>chng 04 jan 22, co.hevmol is current value, we want to save old
		 * value, which is in co.co_save */
		/*float oldoh = co.hevmol[ipOH];*/
		float oldoh = co.co_save[ipOH];
		lgReturn = TRUE;

		/* get new numbers - take average of old and new */
		/* >>chng 03 sep 16, only use damper for molecular species,
		 * not ion/atoms */
		for ( i=0; i < NUM_COMOLE_CALC; i++ )
		{
			/* parameter to mix old and new,
			 * damper is fraction of old solution to take */
			float damper = 0.2f;

			/* >>chng 04 sep 11, correct for den chng in var den models */
			/* the ratio of densities is unity for constant density model, but in const pres or other
			 * models, account for change in entire density scale */
			co.co_save[i] *= dense.gas_phase[ipHYDROGEN] / hden_save_prev_call;

			/* if co.co_save is zero then don't take average of it and 
			 * new solution */
			if( co.co_save[i]< SMALLFLOAT )
			{
				co.co_save[i] = co.hevmol[i];
			}
			else
			{
				/* >>chng 04 jan 24, add logic to check how different old and
				 * new solutions are.  if quite different then take average
				 * in the log */
				double ratio = (double)co.hevmol[i]/(double)co.co_save[i];
				ASSERT( co.co_save[i]>0. && co.hevmol[i]>=0. );
				if( fabs(ratio-1.) < 1.5 ||
					fabs(ratio-1.) > 0.66 )
				{
					/* this branch moderate changes */
					co.hevmol[i] = co.hevmol[i]*(1.f - damper) + 
						co.co_save[i]*damper;
				}
				else
				{
					/* this branch - large changes take average in the log */
					co.hevmol[i] = (float)pow(10. , 
						(log10(co.hevmol[i]) + log10(co.co_save[i]))/2. );
				}
				co.co_save[i] = co.hevmol[i];
			}
		}
		/* remember the current H density */
		hden_save_prev_call = dense.gas_phase[ipHYDROGEN];

		/* >>chng 05 jul 14, add logic to not be as fine a tolerance when
		 * only trace amounts of molecules are present */
		if( oldoh/SDIV(dense.gas_phase[ipOXYGEN]) < 1e-10  )
			conv_fac = 3.;
		else
			conv_fac = 1.;

		if(fabs(oldoh/SDIV(co.hevmol[ipOH])-1.) < COTOLER_MOLAV*conv_fac )
		{
			lgReturn = TRUE;
		}
		else
		{
			/* say we are not converged */
			lgReturn = FALSE;
			sprintf( chReason , "OH change from %.3e to %.3e",
				oldoh , 
				co.hevmol[ipOH]);
		}
	}
	else
		TotalInsanity();

	/* also not converged if co exceeds c */
	if( co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]) > 1.+FLT_EPSILON )
	{
		lgReturn = FALSE;
		sprintf( chReason , "CO/C >1, value is %e",
			co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]) );
	}

	if( (trace.lgTrace && trace.lgTr_CO_Mole) || DEBUG_MOLECAVER )
	{
		fprintf( ioQQQ, " DEBUG lgMolecAver converged? %c" ,TorF(lgReturn) );
		fprintf( ioQQQ, " CO/C=%9.1e", co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]) );
		fprintf( ioQQQ, " OH/O=%9.1e", co.hevmol[ipOH]/SDIV(dense.gas_phase[ipOXYGEN]) );
		fprintf( ioQQQ, "\n" );
	}

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

	return lgReturn;
}


/*CO_drive main driver for heavy molecular equilibrium routines      */
void CO_drive(void)
{
	int lgNegPop, 
	  lgZerPop,
	  lgFirstTry;
	long int i;
	int lgConverged;
	/* count how many times through loop */
	long int loop;
	long int nelem , ion;


	/* this will give reason CO not converged */
	char chReason[100];
	/*static int lgMoleZeroed=TRUE;*/

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

	/* h2lim is smallest ratio of H2/HDEN for which we will
	 * even try to invert heavy element network */
	if( hmi.lgNoH2Mole )
	{
		/* H2 network is turned off, ignore this limit */
		h2lim = 0.;
	}
	else
	{
		/* this limit is important since the co molecular network is first turned
		 * on when this is hit.  It's conservation law will only ever include the
		 * initital O, O+, and C, C+, so must not start before nearly all
		 * C and O is in these forms */
		h2lim = 1e-15;
		/* >>chng 05 jul 15, from 1e-15 to 1e-12, see whether results are stable,
		 * this does include CO in the H+ zone in orion_hii_pdr
		 * a problem is that, during search phase, the first temp is 4000K and the
		 * H2 abundance is larger than it will be at the illuminated face.  try to
		 * avoid turning H2 on too soon */
		h2lim = 1e-12;
	}

	/* do we want to try to calculate the CO? 
	 * >>chng 03 sep 07, add first test, so that we do not turn off
	 * heavy element network if it has ever been turned on
	 * >>chng 04 apr 23, change logic so never done when temp > 20000K */
	/* make series of tests setting co.lgCODoCalc for simplicity */
	co.lgCODoCalc = TRUE;
	/* too hot - don't bother */
	if( phycon.te > 20000. )
		co.lgCODoCalc = FALSE;

	/* molecules have been turned off */
	if( co.lgNoCOMole )
		co.lgCODoCalc = FALSE;

	/* a critical element has been turned off */
	if( dense.lgElmtOn[ipCARBON]/* *dense.lgElmtOn[ipNITROGEN]*/*dense.lgElmtOn[ipOXYGEN] ==  0 )
		co.lgCODoCalc = FALSE;

	/* >>chng 04 jul 20, do not attempt co if no O or C atoms present,
	 * since co net needs their atomic data */
	if( dense.IonLow[ipCARBON]>0 || dense.IonLow[ipOXYGEN]>0 )
		co.lgCODoCalc = FALSE;

	if( iteration!=co.iteration_co && 
		hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim )
		co.lgCODoCalc = FALSE;

	if( dense.lgElmtOn[ipSILICON] )
	{
		/* >>chng 04 jun 28, add test on atomic si, some hii region
		 * models crashed due to matrix failure when co network turned
		 * on with atomic si at very small abundance */
		if( dense.xIonDense[ipSILICON][0]/dense.gas_phase[ipSILICON] < 1e-15 )
			co.lgCODoCalc = FALSE;
	}

	/* this branch zero everything out and return */
	if( !co.lgCODoCalc )
	{
		/* these are heavy - heavy charge transfer rate that will be needed
		 * for recombination of Si+, S+, and C+ */
		mole.xMoleChTrRate[ipSILICON][1][0] = (float)(dense.xIonDense[ipMAGNESIUM][0]*
			rate_COmole.SIP_MG_SI_MGP + 
			dense.xIonDense[ipIRON][0]*rate_COmole.SIP_FE_SI_FEP);

		mole.xMoleChTrRate[ipSULPHUR][1][0] = (float)(dense.xIonDense[ipMAGNESIUM][0]*
			rate_COmole.SP_MG_S_MGP + 
			dense.xIonDense[ipIRON][0]*rate_COmole.SP_FE_S_FEP);

		mole.xMoleChTrRate[ipCARBON][1][0] = (float)(dense.xIonDense[ipMAGNESIUM][0]*
			rate_COmole.CP_MG_C_MGP + 
			dense.xIonDense[ipIRON][0]*rate_COmole.CP_FE_C_FEP);

		/* do we need to zero out the arrays and vars? */
		if( !lgMoleZeroed )
		{
			lgMoleZeroed = TRUE;
			for ( i=0; i<NUM_COMOLE_CALC; ++i )
			{
				co.hevmol[i] = 0.;
			}
			/* heating due to CO photodissociation */
			thermal.heating[0][9] = 0.;
			/* CO cooling */
			CoolHeavy.C12O16Rot = 0.;
			CoolHeavy.dC12O16Rot = 0.;
			CoolHeavy.C13O16Rot = 0.;
			CoolHeavy.dC13O16Rot = 0.;
			co.CODissHeat = 0.;
			co.rate_OH_dissoc = 0.;

			for ( i=0; i < nCORotate; i++ )
			{
				C12O16Rotate[i].PopLo = 0.;
				C13O16Rotate[i].PopLo = 0.;
				/* population of upper level */
				C12O16Rotate[i].PopHi = 0.;
				C13O16Rotate[i].PopHi = 0.;
				/* population of lower level with correction for stim emission */
				C12O16Rotate[i].PopOpc = 0.;
				C13O16Rotate[i].PopOpc = 0.;
				/* following two heat exchange excitation, deexcitation */
				C12O16Rotate[i].cool = 0.;
				C12O16Rotate[i].heat = 0.;
				C13O16Rotate[i].cool = 0.;
				C13O16Rotate[i].heat = 0.;
				/* intensity of line */
				C12O16Rotate[i].xIntensity = 0.;
				C13O16Rotate[i].xIntensity = 0.;
				/* number of photons emitted in line */
				C12O16Rotate[i].phots = 0.;
				C13O16Rotate[i].phots = 0.;
			}
			for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
			{
				for( ion=0; ion<nelem+2; ++ion )
				{
					mole.source[nelem][ion] = 0.;
					mole.sink[nelem][ion] = 0.;
				}
			}
		}
#		ifdef DEBUG_FUN
		fputs( " <->CO_solve()\n", debug_fp );
#		endif
		return;
	}

	/* this flag is used to stop calculation due to negative abundances */
	loop = 0;
	lgNegPop = FALSE;
	lgConverged = lgMolecAver("SETUP",chReason);
	do
	{
		/*oh_h_o_h2ld = co.hevmol[ipOH];*/
		CO_solve(&lgNegPop,&lgZerPop );
		/* this one takes average first, then updates molecular densities in co.hevmol,
		 * returns true if change in OH density is less than macro COTOLER_MOLAV set above*/
		lgConverged = lgMolecAver("AVER",chReason);
		if( CODEBUG )
			fprintf(ioQQQ,"codrivbug %.2f %li Neg?:%c\tZero?:%c\tOH new\t%.3e\tCO\t%.3e\tTe\t%.3e\tH2/H\t%.2e\n",
				fnzone ,
				loop ,
				TorF(lgNegPop) , 
				TorF(lgZerPop) , 
				co.hevmol[ipOH] ,
				co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]),
				phycon.te,
				hmi.H2_total/dense.gas_phase[ipHYDROGEN]);

		++loop;
	}   while (
			  /* these are a series of conditions to stop this loop - 
			   * has the limit to the number of loops been reached? */
			  (loop < LUPMAX_CODRIV) && !lgConverged  );
	
	/* say that we have found a soln before */
	lgMoleZeroed = FALSE;

	/* check whether we have converged */
	if( loop == LUPMAX_CODRIV)
	{
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "CO con1" );
		if( CODEBUG )
			fprintf(ioQQQ,"CO_drive not converged\n");
	}

	/* this flag says whether this is the first zone we are trying
	 * to converge the molecules - there will be problems during initial
	 * search so do not print anything in this case */
	lgFirstTry = (nzone==co.co_nzone && iteration==co.iteration_co);

	/* did the molecule network have negative pops? */
	if( lgNegPop )
	{
		if( conv.lgSearch && (hmi.H2_total/dense.gas_phase[ipHYDROGEN] <h2lim) &&
			(iteration==1) )
		{
			/* we are in search phase during the first iteration, 
			 * and the H2/H ratio has fallen
			 * substantially below the threshold for solving the CO network.
			 * Turn it off */
			CO_Init();
		}
		else
		{
			if( called.lgTalk && !lgFirstTry )
			{
				conv.lgConvPops = FALSE;
				fprintf( ioQQQ, 
						" PROBLEM CO_drive failed to converge1 due to negative pops, zone=%.2f,  CO/C=%.2e  OH/C=%.2e H2/H=%.2e\n", 
						fnzone, 
						co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]),
						co.hevmol[ipOH]/SDIV(dense.gas_phase[ipCARBON]),
						hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
				ConvFail( "pops" , "CO" );
			}
		}
	}

	/* this test, hit upper limit to number of iterations */
	else if( loop == LUPMAX_CODRIV )
	{
		/* do not print this if we are in the first zone where molecules are turned on */
		if( called.lgTalk && !lgFirstTry )
		{
			conv.lgConvPops = FALSE;
			fprintf( ioQQQ, 
				" PROBLEM CO_drive failed to converge2 in %li iter, zone=%.2f, CO/C=%.2e negpop?%1c reason:%s\n", 
					 loop,
					 fnzone, 
					 co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]), 
					 TorF(lgNegPop) ,
					 chReason);
			ConvFail( "pops" , "CO" );
		}
	}

	/* make sure we have not used more than all the CO */
	ASSERT(conv.lgSearch || co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]) <= 1.01 ||
		/* this is true when loop did not converge co */
		loop == LUPMAX_CODRIV );
	/*fprintf(ioQQQ,"ratiooo\t%c\t%.2f\t%f\n", TorF(conv.lgSearch),
		fnzone , co.hevmol[ipCO]/dense.gas_phase[ipCARBON] );*/

	/* these are the elements whose converge we check */
	/* this is a self consistency check, for converged soln */
	/* >>chng 04 dec 02, this test is turn  back on - don't know why it was turned off */
	if(0) {
		double sum_ion , sum_mol;
		sum_ion = dense.xIonDense[ipCARBON][0] + dense.xIonDense[ipCARBON][1];
		sum_mol = co.hevmol[ipATC] + co.hevmol[ipCP];
		if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
		{
			/*fprintf(ioQQQ,
				"non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n",
				5, fnzone, sum_ion , sum_mol);*/
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "CO con2" );
			conv.BadConvIoniz[0] = sum_ion;
			conv.BadConvIoniz[1] = sum_mol;
			if( CODEBUG )
				fprintf(ioQQQ,"CO_drive not converged\n");
		}
		sum_ion = dense.xIonDense[ipOXYGEN][0] + dense.xIonDense[ipOXYGEN][1];
		sum_mol = co.hevmol[ipATO] + co.hevmol[ipOP];
		if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
		{
			/*fprintf(ioQQQ,
				"non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n",
				7, fnzone, sum_ion , sum_mol);*/
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "CO con3" );
			conv.BadConvIoniz[0] = sum_ion;
			conv.BadConvIoniz[1] = sum_mol;
			if( CODEBUG )
				fprintf(ioQQQ,"CO_drive not converged\n");
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->CO_drive()\n", debug_fp );
#	endif
	return;
}
#undef	NCONV
#undef	COTOLER_MOLAV
#undef	CODEBUG
#undef	LUPMAX_CODRIV


