/* This file is part of Cloudy and is copyright (C) 1978-2004 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_step */
/*CO_drive - public routine, calls CO_step to converge molecules */
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "hmi.h"
#include "mole.h"
#include "converge.h"
#include "phycon.h"
#include "trace.h"
#include "thermal.h"
#include "called.h"
#include "coolheavy.h"
#include "co.h"

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

/* 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_step,
 * returns TRUE if converged */
static int lgMolecAver(
	/* job == "SETUP", set things up during first call,
	 * job == "AVER" take average of arrays */
	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;

#	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 */
		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-03f, 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, 6.79E+04f, 1.58E+01f, 1.79E+02f, 2.62E-10f,
				1.79E+03f, 1.00E+00f, 1.13E+05f, 3.56E-04f, 2.25E+00f,
				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 */
	
			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 );
			ASSERT( 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
				&&	ipOCNP  ==  31);
			ASSERT( 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 );
			ASSERT( 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
				&&  ipCP    ==  61
				&&  ipOP    ==  62 );
			ASSERT( ipSIP   ==  63
				&&  ipNP    ==  64
				&&  ipSP    ==  65
				&&  ipATC   ==  66
				&&  ipATO   ==  67
				&&  ipATSI  ==  68
				&&  ipATN   ==  69
				&&  ipATS   ==  70
				&&  NUM_COMOLE_CALC==71 );

			/* 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];
		
			/* 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," 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," lgMolecAver iteration %li zone %li resetting since first call on new iteration. H2/N=%.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," 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];
		/*fprintf(ioQQQ,"DEBUG sil mol\t%.2f\t%.3e\t%.3e\n",
			fnzone , 
			co.hevmol[ipATSI],
			co.hevmol[ipSIP]);*/

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

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

	/* 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]) );
	}

#	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;
	int lgDoCalc;
	double h2lim;
	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;
	}

	/* 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 lgDoCalc for simplicity */
	lgDoCalc = TRUE;
	/* too hot - don't bother */
	if( phycon.te > 20000. )
		lgDoCalc = FALSE;

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

	/* a critical element has been turned off */
	if( dense.lgElmtOn[ipCARBON]/* *dense.lgElmtOn[ipNITROGEN]*/*dense.lgElmtOn[ipOXYGEN] ==  0 )
		lgDoCalc = 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 )
		lgDoCalc = FALSE;

	if( iteration!=co.iteration_co && 
		hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim )
		lgDoCalc = 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 )
			lgDoCalc = FALSE;
	}

	/* this branch zero everything out and return */
	if( !lgDoCalc )
	{
		/* 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_step()\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_step(&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 %li Neg?:%c\tZero?:%c\tOH new\t%.3e\tCO\t%.3e\tTe\t%.3e\tH2/H\t%.2e\n",
					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 converge due to negative pops, zone=%4ld,  CO/C=%.2e  OH/C=%.2e H2/H=%.2e\n", 
						nzone, 
						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" );
			}
		}
	}
#	if 0
	/* >>chng 03 sep 23, don't comment on this - it happens when CO fully formed */
	/* did the molecule network have zero pops? */
	else if ( lgZerPop )
	{
		if ( called.lgTalk && !lgFirstTry && (hmi.H2_total/dense.gas_phase[ipHYDROGEN] > h2lim))
		{
			fprintf( ioQQQ, 
					 " PROBLEM CO_drive failed to converge due to zero pops, zone=%4ld,  CO/C=%.2e  OH/C=%.2e \n", 
					 nzone, 
					 co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]),
					 co.hevmol[ipOH]/SDIV(dense.gas_phase[ipCARBON]) );
		}
	}
#	endif

	/* 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 converge in %li iter, zone=%ld, CO/C=%.2e negpop?%1c reason:%s\n", 
					 loop,
					 nzone, 
					 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 */
	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

