/*cdDrive main routine to call cloudy under all circumstances) */
/*cdReasonGeo wrte why the model stopped and type of geometry on io file */
/*cdWarnings write all warnings entered into comment stack */
/*cdEms obtain the local emissivity for a line, for the last computed zone */
/*cdColm get the column density for a constituent  */
/*cdLine get the predicted line intensity, also index for line in stack */
/*cdLine_ip get the predicted line intensity, using index for line in stack */
/*cdDLine get the predicted emergent line intensity, also index for line in stack */
/*cdCautions print out all cautions after calculation, on arbitrary io unit */
/*cdLastTemp routine to query results and return temperature of last zone */
/*cdTimescales returns thermal, recombination, and H2 foramtion timescales */
/*cdSurprises print out all surprises on arbitrary unit number */
/*cdNotes print stack of notes about current calculation */
/*cdGetPres routine to query results and return pressure of last zone */
/*cdTalk tells the code whether to print results or be silent */
/*cdOutp redirect output to arbitrary Fortran unit number */
/*cdRead routine to read in command lines when cloudy used as subroutine */
/*cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit */
/*cdIonFrac get ionization fractions for a constituent */
/*cdTemp get mean electron temperature for any element */
/*cdGetCooling routine to query results and return cooling of last zone */
/*cdGetHeating routine to query results and return heating of last zone */
/*cdNoExec call this routine to tell code not to actually execute */
/*cdDate - puts date of code into string */
/*cdVersion produces string that gives version number of the code */
/*cdExecTime any routine can call this, and find the time since cdInit was called */
/*cdPrintCommands( FILE *) prints all input commands into file */
/*cdDrive main routine to call cloudy under all circumstances) */
/*cdNwcns get the number of cautions and warnings, to tell if calculation is ok */
/*debugLine provides a debugging hook into the main line array  */
/*cdPrtWL print line wavelengths in Angstroms in the standard format - just a wrapper */
/*cdEms_ip obtain the local emissivity for a line with known index */

#include "cddefines.h"
#include "vary.h"
#include "trace.h"
#include "converge.h"
#include "norm.h"
#include "pressure.h"
#include "heat.h"
#include "linesave.h"
#include "prt.h"
#include "colden.h"
#include "hevmolec.h"
#include "radius.h"
#include "elementnames.h"
#include "mean.h"
#include "phycon.h"
#include "called.h"
#include "input.h"
#include "version.h"
#include "cooling.h"
#include "optimr.h"
#include "dooptimize.h"
#include "timesc.h"
#include "cloudy.h"
#include "aaaa.h"
#include "warnings.h"
#include "cddrive.h"

/* >>chng 02 may 02, change SMALLFLOAT (1.e-35f) -> DELTA (1e-20f)
 * to prevent FP overflows in wavelength comparisons, PvH */
#define DELTA 1.e-20f

/*************************************************************************
 *
 * cdDrive - main routine to call cloudy - returns 0 if all ok, 1 if problems
 *
 ************************************************************************/

int cdDrive(void )
{
	int lgBAD;

#	ifdef DEBUG_FUN
	fputs( "<+>cdDrive()\n", debug_fp );
#	endif
	/*********************************
	 * main routine to call cloudy   *
	 *********************************/

	/* this is set FALSE when code loaded, set TRUE when cdInit called,
	 * this is check that cdInit was called first */
	if( !lgcdInitCalled )
	{
		printf(" cdInit was not called first - this must be the first call.\n");
		puts( "[Stop in cdDrive]" );
		cdEXIT(1);
	}

	if( trace.lgTraceInput )
	{
		fprintf( ioQQQ, 
			"cdDrive: lgOptimr=%1i lgVaryOn=%1i lgNoVary=%1i input.nSave:%li\n",
			optimr.lgOptimr , vary.lgVaryOn , vary.lgNoVary, input.nSave );
	}

	/* should we call cloudy, or the optimization driver? */
	/* possible to have VARY on line without OPTIMIZE being set */
	if( optimr.lgOptimr && vary.lgVaryOn )
	{
		vary.lgVaryOn = TRUE;
	}
	else
	{
		vary.lgVaryOn = FALSE;
	}

	/* this is option to have command line saying "no vary"
	 * to turn off optimizer */
	if( vary.lgNoVary )
		vary.lgVaryOn = FALSE;

	if( vary.lgVaryOn )
	{
		if( trace.lgTraceInput )
			fprintf( ioQQQ, "cdDrive: calling DoOptimize\n");
		/* option to drive optimizer set if OPTIMIZE was in input stream */
		lgBAD = DoOptimize();
	}
	else
	{
		if( trace.lgTraceInput )
			fprintf( ioQQQ, "cdDrive: calling cloudy\n");
		/* optimize did not occur, only compute one model, call cloudy */
		lgBAD = cloudy();
	}

	/* reset flag saying that cdInit has not been called */
	lgcdInitCalled = FALSE;

	if( called.lgBusted || lgBAD )
	{
		/* lgBusted set true if something wrong, so return lgBAD false. */
#		ifdef DEBUG_FUN
		fputs( " <->cdDrive()\n", debug_fp );
#		endif
		return(1);
	}
	else
	{
		/* everything is ok, return 0 */
#		ifdef DEBUG_FUN
		fputs( " <->cdDrive()\n", debug_fp );
#		endif
		return(0);
	}
}

/*************************************************************************
 *
 * cdReasonGeo wrte why the model stopped and type of geometry on io file 
 *
 ************************************************************************/


/*cdReasonGeo wrte why the model stopped and type of geometry on io file */
void cdReasonGeo(FILE * ioOUT)
{

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

	/*this is the reason the calculation stopped*/
	fprintf( ioOUT, "%s", warnings.chRgcln[0] );
	fprintf( ioOUT , "\n" );
	/* this is the geometry */
	fprintf( ioOUT, "%s", warnings.chRgcln[1] );
	fprintf( ioOUT , "\n" );

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


/*************************************************************************
 *
 * cdWarnings write all warnings entered into comment stack 
 *
 ************************************************************************/

/*cdWarnings write all warnings entered into comment stack */

void cdWarnings(FILE *ioPNT )
{
	long int i;

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

	if( warnings.nwarn > 0 )
		{
			for( i=0; i < warnings.nwarn; i++ )
			{
				/* these are all warnings that were entered in comment */
				fprintf( ioPNT, "%s", warnings.chWarnln[i] );
				fprintf( ioPNT, "\n" );
			}
		}


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


/*************************************************************************
 *
 * cdCautions print out all cautions after calculation, on arbitrary io unit 
 *
 ************************************************************************/

/*cdCautions print out all cautions after calculation, on arbitrary io unit */

void cdCautions(FILE * ioOUT)
{
	long int i;

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

	if( warnings.ncaun > 0 )
	{
		for( i=0; i < warnings.ncaun; i++ )
		{
			fprintf( ioOUT, "%s", warnings.chCaunln[i] );
			fprintf( ioOUT, "\n" );
		}
	}

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

/*************************************************************************
 *
 * cdTimescales returns thermal, recombination, and H2 foramtion timescales 
 *
 ************************************************************************/

void cdTimescales(
	/* the thermal cooling timescale */
	double *TTherm , 
	/* the hydrogen recombination timescale */
	double *THRecom , 
	/* the H2 formation timescale */
	double *TH2 )
{

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

	/* these were all evaluated in AgeCheck, which was called by PrtComment */

	/* thermal or cooling timescale */
	*TTherm = timesc.ttherm;

	/* the hydrogen recombination timescale */
	*THRecom = timesc.threc;

	/* the H2 formation timescale */
	*TH2 = timesc.BigH2MoleForm;

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


/*************************************************************************
 *
 * cdSurprises print out all surprises on arbitrary unit number 
 *
 ************************************************************************/

/*cdSurprises print out all surprises on arbitrary unit number */

void cdSurprises(FILE * ioOUT)
{
	long int i;

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

	if( warnings.nbang > 0 )
	{
		for( i=0; i < warnings.nbang; i++ )
		{
			fprintf( ioOUT, "%s", warnings.chBangln[i] );
			fprintf( ioOUT, "\n" );
		}
	}


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


/*************************************************************************
 *
 * cdNotes print stack of notes about current calculation 
 *
 ************************************************************************/

/*cdNotes print stack of notes about current calculation */

void cdNotes(FILE * ioOUT)
{
	long int i;

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

	if( warnings.nnote > 0 )
		{
			for( i=0; i < warnings.nnote; i++ )
			{
				fprintf( ioOUT, "%s", warnings.chNoteln[i] );
				fprintf( ioOUT, "\n" );
			}
		}

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

/*************************************************************************
 *
 * cdGetCooling routine to query results and return cooling of last zone 
 *
 ************************************************************************/

/*cdGetCooling routine to query results and return cooling of last zone */
double cdGetCooling(void) /* return cooling for last zone */
{
	return (cooling.ctot);
}

/*************************************************************************
 *
 * cdVersion - puts version number of code into string 
 * incoming string must have at least 8 char and will become null
 * terminated string
 *
 ************************************************************************/

void cdVersion(char chString[] ) 
{
	/* first call routine that sets var to make sure we get valid numbers */
	AAAA();
	strcpy( chString , version.chVersion );
	return ;
}

/*************************************************************************
 *
 * cdDate - puts date of code into string 
 * incoming string must have at least 8 char and will become null
 * terminated string
 *
 ************************************************************************/

void cdDate(char chString[] ) 
{
	/* first call routine that sets var to make sure we get valid numbers */
	AAAA();
	strcpy( chString , version.chDate );
	return ;
}


/*************************************************************************
 *
 * cdGetHeating routine to query results and return heating of last zone
 *
 ************************************************************************/

/*cdGetHeating routine to query results and return heating of last zone */

double cdGetHeating(void) /* return heating for last zone */
{
	return (heat.htot);
}


/*************************************************************************
 *
 * cdNoExec call this routine to tell code not to actually execute
 *
 ************************************************************************/

/*cdNoExec call this routine to tell code not to actually execute */
#include "noexec.h"

void cdNoExec(void)
{

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

	/* this is an option to have the code read in all input quanties
	 * ot to NOT execute the actual model - to check on input parameters
	 * */
	noexec.lgNoExec = TRUE;


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


/*************************************************************************
 *
 * cdSetExecTime routine to initialize variables keeping track of time at start of calculation
 *
 ************************************************************************/

/*cdExecTime any routine can call this, and find the time since cdInit was called */
#include <time.h>
/* will be used to save initial time */
static clock_t before;
/* this must be true when cdExeTime is called, and is set true
 * when SetExecTime is called */
static int lgCalled=FALSE;

/* this routine is called by cdinit when everything is initialized,
 * so that every time cdExecTime is called the elapsed time is returned */
void cdSetExecTime(void)
{
	/* save startup time */
	before = clock();
	lgCalled = TRUE;
}

/*************************************************************************
 *
 * cdExecTime any routine can call this, and find the time since cdInit was called
 *
 ************************************************************************/

double cdExecTime(void)
{

	/* check that we were properly initialized */
	if( lgCalled )
	{
		/* this is elapsed time in seconds */
		return(  (double)(clock() - before) / (double)CLOCKS_PER_SEC )  ;
	}
	else
	{
		/* this is a big problem, we were asked for the elapsed time but
		 * the timer was not initialized by calling SetExecTime */
		fprintf( ioQQQ, " cdExecTime was called before SetExecTime, impossible.\n" );
		fprintf( ioQQQ, "Sorry.\n" );
		puts( "[Stop in cdExecTime]" );
		cdEXIT(1);
	}
}

/*************************************************************************
 *
 * cdPrintCommands prints all input commands into file
 *
 ************************************************************************/

/* cdPrintCommands( FILE *)
 * prints all input commands into file */
void cdPrintCommands( FILE * ioOUT )
{
	long int i;
	fprintf( ioOUT, " Input commands follow:\n" );
	fprintf( ioOUT, "c ======================\n" );

	for( i=0; i <= input.nSave; i++ )
	{
		fprintf( ioOUT, "%s\n", input.chCardSav[i] );
	}
	fprintf( ioOUT, "c ======================\n" );
}


/*************************************************************************
 *
 * cdEms obtain the local emissivity for a line, for the last computed zone
 *
 ************************************************************************/

long int cdEmis(
	/* return value will be index of line within stack,
	 * negative of number of lines in the stack if the line could not be found*/
	/* 4 char null terminated string label */
	char *chLabel,
	/* line wavelength */
	float wavelength, 
	/* the vol emissivity of this line in last computed zone */
	double *emiss )
{
	/* use following to store local image of character strings */
	char chCARD[81];
	char chCaps[5];
	long int j;
	double a;
	float errorwave;

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

	/* routine returns the emissivity in the desired line
	 * only used internally by the code, to do punch lines structure */

	strcpy( chCARD, chLabel );

	/* make sure chLabel is all caps */
	caps(chCARD);/* convert to caps */

	/* get the error assocated with 4 significant figures */
	if( wavelength > 0. )
	{
		a = log10( wavelength );
		if( wavelength>=1. )
		{
			a -= floor(a);
		}
		else
		{
			a += floor(a);
		}
		a = pow(10.,a);
		/* a is now 1 - 10, so this corresponds to 4 sig fig */
		errorwave = (float)(0.00105/a);
	}
	else
	{
		errorwave = 1e-4f;
	}

	for( j=0; j < LineSave.nsum; j++ )
	{
		/* change chLabel to all caps to be like input chAmLab */
		cap4(chCaps , (char*)LineSv[j].chALab);

		/* check wavelength and chLabel for a match */
		if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength)<errorwave && 
			strcmp(chCaps,chCARD) == 0 )
		{
			/* match, so set emiss to emissivity in line */
			*emiss = LineSv[j].emslin;
			
#			ifdef DEBUG_FUN
			fputs( " <->cdEms()\n", debug_fp );
#			endif
			/* and announce success by returning line index within stack */
			return j ;
		}
	}

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

	/* we fell through without finding the line - return false */
	return -LineSave.nsum;
}


/*************************************************************************
 *
 * cdEms_ip obtain the local emissivity for a line with known index
 *
 ************************************************************************/

void cdEmis_ip(
	/* index of the line in the stack */
	long int ipLine, 
	/* the vol emissivity of this line in last computed zone */
	double *emiss )
{

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

	/* routine returns the emissivity in the desired line
	 * used to do punch lines structure */

	/* match, so set emiss to emissivity in line */
	assert( ipLine >= 0 && ipLine < LineSave.nsum );
	*emiss = LineSv[ipLine].emslin;

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

	return;
}

/*************************************************************************
 *
 * cdColm get the column density for a constituent 
 *
 ************************************************************************/

int cdColm(
	/* return value is zero if all ok, 1 if errors happened */
	/* 4-char + eol string that is first 
	 * 4 char of element name as spelled by cloudy */
	char *chLabel,	

	/* integer stage of ionization, 1 for atom, 0 for CO or H2 */
	long int ion,

	/* the theoretical column density derived by the code */
	double *theocl )	
{
	long int nelem;
	/* use following to store local image of character strings */
	char chCARD[81];

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

	strcpy( chCARD, chLabel );
	/* convert element label to all caps */
	caps(chCARD);

	/* zero ionization stage has special meaning.  The quanties recognized are
	 * two molecules, "H2  ", "OH  ", and "CO  " 
	 * "CII*" excited state C+ */
	if( ion < 0 )
	{
		fprintf( ioQQQ, " cdColm called with insane ion, =%li\n", 
		  ion );
#		ifdef DEBUG_FUN
		fputs( " <->cdColm()\n", debug_fp );
#		endif
		return(1);
	}
	else if( ion == 0 )
	{
		/* this case molecular column density */
		/* want the molecular hydrogen column density */
		if( strcmp( chCARD , "H2  " )==0 )
		{
			*theocl = coldenCom.colden[ipCOLH2];
		}

		/* carbon monoxide column density */
		else if( strcmp( chCARD , "CO  " )==0 )
		{
			*theocl = hevmolec.hevcol[ipCO];
		}

		/* OH column density */
		else if( strcmp( chCARD , "OH  " )==0 )
		{
			*theocl = hevmolec.hevcol[ipOH];
		}

		/* CII^* column density, population of J=3/2 upper level of split ground term */
		else if( strcmp( chCARD , "CII*" )==0 )
		{
			*theocl = coldenCom.C2Colden[1];
		}

		/* clueless as to what was meant - bail */
		else
		{
			fprintf( ioQQQ, " cdColm called with unknown element chLabel, =%4.4s\n", 
			  chLabel );
#			ifdef DEBUG_FUN
			fputs( " <->cdColm()\n", debug_fp );
#			endif
			return(1);
		}
	}
	else
	{
		/* this case, ionization stage of some element */
		/* find which element this is */
		nelem = 0;
		while( nelem < LIMELM && 
			strncmp(chCARD,elementnames.chElementNameShort[nelem],4) != 0 )
		{
			++nelem;
		}

		/* this is true if we have one of the first 30 elements in the label */
		if( nelem < LIMELM )
		{

			/* sanity check - does this ionization stage exist? */
			if( ion > nelem + 2 )
			{
				fprintf( ioQQQ, 
				  " cdColm asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 
				  ion, chLabel );
				return(1);
			}

			/* the column density */
			*theocl = IonMeans.xIonMeans[0][nelem][ion];
		}
		else
		{
			fprintf( ioQQQ, 
			  " cdColm did not understand this combinatino of ion %4ld and element %4.4s.\n", 
			  ion, chLabel );
			return(1);
		}
	}

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


/*************************************************************************
 *
 *cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit 
 *
 ************************************************************************/

void cdErrors(FILE *ioOUT)
{
	long int nc, 
	  nn, 
	  npe, 
	  ns, 
	  nte, 
	  nw ,
	  nIone,
	  nEdene;
	int lgAbort;

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

	/* first check for number of warnings, cautions, etc */
	cdNwcns(&lgAbort,&nw,&nc,&nn,&ns,&nte,&npe, &nIone, &nEdene );

	/* only say something is one of these problems is nonzero */
	if( nw || nc || nte || npe ||	nIone || nEdene || lgAbort )
	{
		/* say the title of the model */
		fprintf( ioOUT, "%75.75s\n", input.chTitle );

		if( lgAbort )
			fprintf(ioOUT," Calculation ended with abort!\n");

		/* print warnings on the io unit */
		if( nw != 0 )
		{
			cdWarnings(ioOUT);
		}

		/* print cautions on the io unit */
		if( nc != 0 )
		{
			cdCautions(ioOUT);
		}

		if( nte != 0 )
		{
			fprintf( ioOUT , "Te failures=%4ld\n", nte );
		}

		if( npe != 0 )
		{
			fprintf( ioOUT , "Pressure failures=%4ld\n", npe );
		}

		if( nIone != 0 )
		{
			fprintf( ioOUT , "Ionization failures=%4ld\n", nte );
		}

		if( nEdene != 0 )
		{
			fprintf( ioOUT , "Electron density failures=%4ld\n", npe );
		}
	}

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

/*************************************************************************
 *
 *cdGetPres routine to query results and return pressure of last zone 
 *
 ************************************************************************/

void cdGetPres(
		double *PresTotal,  /* total pressure, all forms*/
		double *PresGas, /* gas pressure */
		double *PresRad) /* radiation pressure */
{

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

	*PresGas = pressure.PressureGas;
	*PresRad = pressure.PressureRadiation;
	*PresTotal = pressure.PressureGas + pressure.PressureRadiation;

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

/*************************************************************************
 *
 *cdLastTemp routine to query results and return temperature of last zone 
 *
 ************************************************************************/


double cdLastTemp(void)
{

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


#	ifdef DEBUG_FUN
	fputs( " <->cdLastTemp()\n", debug_fp );
#	endif
	return( phycon.te);
}

/*************************************************************************
 *
 *cdIonFrac get ionization fractions for a constituent
 *
 ************************************************************************/

int cdIonFrac(
	/* four char string, null terminzed, giving the element name */
	char *chLabel, 
	/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
	long int IonStage, 
	/* will be fractional ionization */
	double *fracin, 
	/* how to weight the average, must be "VOLUME" or "RADIUS" */
	char *chWeight ,
	/* if true then weighting also has electron density, if false then only volume or radius */
	int lgDensity ) 
	/* return value is 0 if element was found, non-zero if failed */
{

	int lgVol;
	long int ip, 
		ipaaa, /* used as index within aaa vector*/
		nelem;
	float aaa[LIMELM + 1];
	/* use following to store local image of character strings */
	char chCARD[81];

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

	strcpy( chCARD, chWeight );
	/* make sure chWeight is all caps */
	caps(chCARD);/* convert to caps */

	/*caps(chWeight);*/

	if( strcmp(chCARD,"RADIUS") == 0 )
	{
		lgVol = FALSE;
	}

	else if( strcmp(chCARD,"VOLUME") == 0 )
	{
		lgVol = TRUE;
	}

	else
	{
		fprintf( ioQQQ, " cdIonFrac: chWeight=%6.6s makes no sense to me\n", 
		  chWeight );
		*fracin = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(1);
	}

	/* first ensure that chLabel is all caps */
	strcpy( chCARD, chLabel );
	/* make sure chWeight is all caps */
	caps(chCARD);/* convert to caps */

	/* find which element this is, nelem is on physical scale, H is 1 */
	nelem = 1;
	while( nelem <= LIMELM && 
		strcmp(chCARD,elementnames.chElementNameShort[nelem-1]) != 0 )
	{
		nelem += 1;
	}
	/* after this loop nelem is atomic number of element, H is 1 */

	/* if element not recognized and above loop falls through, nelem is LIMELM+1 
	 * nelem counter is on physical scale, H = 1 Zn = 30 */
	if( nelem > LIMELM )
	{
		fprintf( ioQQQ, " cdIonFrac called with unknown element chLabel, =%4.4s\n", 
		  chLabel );
		
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(1);
	}

	/* sanity check - does this ionization stage exist? 
	 * both IonStage and nelem are on physical scales, H atoms are 1 1 */

	/* ipaaa will be used as pointer within the aaa array that contains average values,
	 * done this way to prevent lint from falsing in acess to aaa array */
	ipaaa = IonStage - 1;

	/*if( IonStage > nelem + 1 || IonStage < 1 || IonStage > LIMELM+1 )*/
	if( ipaaa > nelem || ipaaa < 0 || ipaaa > LIMELM )
	{
		fprintf( ioQQQ, " cdIonFrac asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 
		  IonStage, chLabel );
		*fracin = -1.;
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(1);
	}

	/* get either volume or radius average, aaa is filled in from 0 */
	if( lgVol )
	{
		/* aaa is dim'd LIMELM+1 so largest arg is LIMELM
		 * 'i' means ionization, not temperature 
		 * nelem is on the physical scale */
		/* last arg says whether to include electron density */
		VolMean('i', nelem,&ip,aaa,lgDensity);
		*fracin = pow(10.f,aaa[IonStage-1]);
	}
	else
	{
		/* last arg says whether to include electron density */
		RadMean('i', nelem,&ip,aaa,lgDensity);
		*fracin = pow(10.f,aaa[IonStage-1]);
	}

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

	/* we succeeded - say so */
	return(0);
}

/*************************************************************************
 *
 * debugLine provides a debugging hook into the main line array 
 *
 ************************************************************************/

 /* following routine provides a debugging hook into the main line array 
 * loops over whole array and finds every line that matches length,
 * the wavelength, the argument to the function
 * put breakpoint inside if test */
long debugLine( float wavelength )
{
	long j, kount;
	double a;
	float errorwave;
	kount = 0;

	/* get the error assocated with 4 significant figures */
	if( wavelength > 0. )
	{
		a = log10( wavelength );
		if( wavelength>=1. )
		{
			a -= floor(a);
		}
		else
		{
			a += floor(a);
		}
		a = pow(10.,a);
		/* a is now 1 - 10, so this corresponds to 4 sig fig */
		errorwave = (float)(0.00105/a);
	}
	else
	{
		errorwave = 1e-4f;
	}

	for( j=0; j < LineSave.nsum; j++ )
	{
		/* check wavelength and chLabel for a match */
		if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength) < errorwave )
		{
			printf("%s\n", LineSv[j].chALab);
			++kount;
		}
	}
	printf(" hits = %li\n", kount );
	return(kount);
}


/*************************************************************************
 *
 *cdLine get the predicted line intensity, also index for line in stack 
 *
 ************************************************************************/

 /* returns array index for line in array stack if we found the line, 
 * or FALSE==0 if we did not find the line */
long int cdLine(char *chLabel, 
		/* wavelength of line in angstroms, not format printed by code */
	  float wavelength, 
	  /* linear intensity relative to normalization line*/
	  double *relint, 
	  /* log of luminosity or intensity of line */
	  double *absint )
{
	char chCaps[5], 
	  chFind[5];
	long int ipobs, 
	  j;
	double a;
	float errorwave;

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

	/* this is zero when cdLine called with cdNoExec called too */
	if( LineSave.nsum == 0 )
	{
		*relint =  0.;
		*absint = 0.;
		return 0;
	}
	assert( norm.ipNormWavL >= 0);
	assert( LineSave.nsum > 0);

	/* change chLabel to all caps */
	cap4(chFind,chLabel);

	/* get the error assocated with 4 significant figures */
	if( wavelength > 0. )
	{
		a = log10( wavelength+FLT_EPSILON);
		if( wavelength>=1. )
		{
			a -= floor(a);
		}
		else
		{
			a += floor(a);
		}
		a = pow(10.,a);
		/* a is now 1 - 10, so this corresponds to 4 sig fig */
		errorwave = (float)(0.00105/a);
	}
	else
	{
		errorwave = 1e-4f;
		a = 0.;
	}

	{
		/*@-redef@*/
		enum{DEBUG=FALSE};
		/*@+redef@*/
		if( DEBUG && fabs(wavelength-1000.)<0.01 )
		{
			fprintf(ioQQQ,"cddrive wl %.4e a %.3e error %.3e\n",
				wavelength, a, errorwave );
		}
	}

	/* now go through entire line stack, do not do 0, which is H-beta and
	 * in stack further down - this is to stop query for H-beta from returning
	 * 0, the flag for line not found */
	for( j=1; j < LineSave.nsum; j++ )
	{

		/* check wavelength and chLabel for a match */
		/* DELTA since often want wavelength of zero */
		if( fabs(LineSv[j].wavelength-wavelength)/MAX2(DELTA,wavelength) < errorwave )
		{

			/* change chLabel to all caps to be like input chALab */
			cap4(chCaps , (char*)LineSv[j].chALab);

			/* now see if labels agree */
			if( strcmp(chCaps,chFind) == 0 )
			{
				/* match, so set pointer */
				ipobs = j;

				/* does the normalization line have a positive intensity*/
				if( LineSv[norm.ipNormWavL].sumlin > 0. )
				{
					*relint = LineSv[ipobs].sumlin/LineSv[norm.ipNormWavL].sumlin*
					  norm.ScaleNormLine;
				}
				else
				{
					*relint = 0.;
				}

				/* return log of current line intensity if it is positive */
				if( LineSv[ipobs].sumlin > 0. )
				{
					*absint = log10(LineSv[ipobs].sumlin) + radius.Conv2PrtInten;
				}
				else
				{
					/* line intensity is actually zero, return small number */
					*absint = -37.;
				}
				
#				ifdef DEBUG_FUN
				fputs( " <->cdLine()\n", debug_fp );
#				endif
				/* we found the line, return pointer to its location */
				return ipobs;
			}
		}
	}

	*absint = 0.;
	*relint = 0.;

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

	/* if we fell down to here we did not find the line */
	/* >>chng 00 sep 02, return negative of total number of lines as debugging aid */
	return -LineSave.nsum;
}

/*************************************************************************
 *
 *cdLine_ip get the predicted line intensity, using index for line in stack 
 *
 ************************************************************************/

void cdLine_ip(long int ipLine, 
	  /* linear intensity relative to normalization line*/
	  double *relint, 
	  /* log of luminosity or intensity of line */
	  double *absint )
{

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

	/* this is zero when cdLine called with cdNoExec called too */
	if( LineSave.nsum == 0 )
	{
		*relint =  0.;
		*absint = 0.;
		return;
	}
	assert( norm.ipNormWavL >= 0);
	assert( LineSave.nsum > 0);

	/* does the normalization line have a positive intensity*/
	if( LineSv[norm.ipNormWavL].sumlin > 0. )
	{
		*relint = LineSv[ipLine].sumlin/LineSv[norm.ipNormWavL].sumlin*
		  norm.ScaleNormLine;
	}
	else
	{
		*relint = 0.;
	}

	/* return log of current line intensity if it is positive */
	if( LineSv[ipLine].sumlin > 0. )
	{
		*absint = log10(LineSv[ipLine].sumlin) + radius.Conv2PrtInten;
	}
	else
	{
		/* line intensity is actually zero, return small number */
		*absint = -37.;
	}

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

	return;
}

/*************************************************************************
 *
 *cdDLine get the predicted emergent line intensity, also index for line in stack 
 *
 ************************************************************************/

 /* returns array index for line in array stack if we found the line, 
 * or FALSE==0 if we did not find the line */
long int cdDLine(char *chLabel, 
		/* wavelength of line as printed by code*/
	  float wavelength, 
	  /* linear intensity relative to normalization line*/
	  double *relint, 
	  /* log of luminosity or intensity of line */
	  double *absint )
{
	char chCaps[5], 
	  chFind[5];
	long int ipobs, 
	  j;
	double a;
	float errorwave;

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

	/* this is zero when cdDLine called with cdNoExec called too */
	if( LineSave.ndsum == 0 )
	{
		*relint =  0.;
		*absint = 0.;
		return 0;
	}
	assert( norm.ipNormWavL >= 0);
	assert( LineSave.nsum > 0);

	/* change chLabel to all caps */
	cap4(chFind,chLabel);

	/* get the error assocated with 4 significant figures */
	if( wavelength > 0. )
	{
		a = log10( wavelength );
		if( wavelength>=1. )
		{
			a -= floor(a);
		}
		else
		{
			a += floor(a);
		}
		a = pow(10.,a);
		/* a is now 1 - 10, so this corresponds to 4 sig fig */
		errorwave = (float)(0.00105/a);
	}
	else
	{
		errorwave = 1e-4f;
	}

	/* now go through entire line stack, do not do 0, which is H-beta and
	 * in stack further down - this is to stop query for H-beta from returning
	 * 0, the flag for line not found */
	for( j=1; j < LineSave.ndsum; j++ )
	{

		/* check wavelength and chLabel for a match */
		/* DELTA since often want wavelength of zero */
		if( fabs(LineDSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength) < errorwave )
		{

			/* change chLabel to all caps to be like input chALab */
			cap4(chCaps , (char*)LineDSv[j].chSMDLab);

			/* now see if labels agree */
			if( strcmp(chCaps,chFind) == 0 )
			{
				/* match, so set array index */
				ipobs = j;

				/* does the normalization line have a positive intensity*/
				if( LineDSv[norm.ipEmerNormWavL].smdlin > 0. )
				{
					*relint = LineDSv[ipobs].smdlin/LineDSv[norm.ipEmerNormWavL].smdlin*
					  norm.ScaleNormLine;
				}
				else
				{
					*relint = 0.;
				}

				/* return log of current line intensity if it is positive */
				if( LineDSv[ipobs].smdlin > 0. )
				{
					*absint = log10(LineDSv[ipobs].smdlin) + radius.Conv2PrtInten;
				}
				else
				{
					/* line intensity is actually zero, return small number */
					*absint = -37.;
				}
				
#				ifdef DEBUG_FUN
				fputs( " <->cdDLine()\n", debug_fp );
#				endif
				/* we found the line, return pointer to its location */
				return ipobs;
			}
		}
	}

	*absint = 0.;
	*relint = 0.;

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

	/* if we fell down to here we did not find the line */
	/* >>chng 00 sep 02, return negative of total number of lines as debugging aid */
	return -LineSave.nsum;
}

/*************************************************************************
 *
 *cdNwcns get the number of cautions and warnings, to tell if calculation is ok
 *
 ************************************************************************/

void cdNwcns(
  /* abort status, this better be false */
  int *lgAbort ,
  /* the number of warnings, cautions, notes, and surprises */
  long int *NumberWarnings, 
  long int *NumberCautions, 
  long int *NumberNotes, 
  long int *NumberSurprises, 
  /* the number of temperature convergence failures */
  long int *NumberTempFailures, 
  /* the number of pressure convergence failures */
  long int *NumberPresFailures,
  /* the number of ionzation convergence failures */
  long int *NumberIonFailures, 
  /* the number of electron density convergence failures */
  long int *NumberNeFailures )
{

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

	/* this would be set true if code ended with abort - very very bad */
	*lgAbort = called.lgBusted;
	/* this sub is called after comment, to find the number of various comments */
	*NumberWarnings = warnings.nwarn;
	*NumberCautions = warnings.ncaun;
	*NumberNotes = warnings.nnote;
	*NumberSurprises = warnings.nbang;

	/* these are counters that were incremented during convergence failures */
	*NumberTempFailures = conv.nTeFail;
	*NumberPresFailures = conv.nPreFail;
	*NumberIonFailures = conv.nIonFail;
	*NumberNeFailures = conv.nNeFail;

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

/*************************************************************************
 *
 * cdOutp redirect output to arbitrary openned file
 *
 ************************************************************************/

void cdOutp( FILE *ioOut)
{

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

	/* ioQQQ is pointer to output fiile */
	ioQQQ = ioOut;

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

/*************************************************************************
 *
 *cdTalk tells the code whether to print results or be silent
 *
 ************************************************************************/

void cdTalk(int lgTOn)
{

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

	called.lgTalk = lgTOn;

	if( lgTOn )
	{
		/* means talk has not been forced off */
		called.lgTalkForce = FALSE;
	}
	else
	{
		/* means talk is forced off */
		called.lgTalkForce = TRUE;
	}

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

/*************************************************************************
 *
 * cdTemp get mean electron temperature for any element 
 *
 ************************************************************************/

/* This routine finds the mean electron temperature for any ionization stage 
 * It returns 0 if it could find the species, 1 if it could not find the species.
 * The first argument is a null terminated 4 char string that gives the element
 * name as spelled by Cloudy.  
 * The second argument is ion stage, 1 for atom, 2 for first ion, etc
 * This third argument will be returned as result,
 * Last parameter is either "VOLUME" or "RADIUS" to give weighting 
 *
 * if the ion stage is zero then the element label will have a special meaning.
 * The string "21CM" is will return the 21 cm temperature.
 * The string "H2  " will return the temperature weighted wrt the H2 density */

int cdTemp(
	/* four char string, null terminzed, giving the element name */
	char *chLabel, 
	/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
	 * 0 means that chLabel is a special case */
	long int IonStage, 
	/* will be temperature */
	double *TeMean, 
	/* how to weight the average, must be "VOLUME" or "RADIUS" */
	char *chWeight ) 
	/* return value is o if things are ok and element was found, 
	 * non-zero if element not found or there are problems */
{

	int lgVol;
	long int ip, 
		ipaaa, /* used as pointer within aaa vector*/
		nelem;
	float aaa[LIMELM + 1];
	/* use following to store local image of character strings */
	char chWGHT[81] , chELEM[81];

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

	/* make sure chWeight is all caps */
	strcpy( chWGHT, chWeight );
	caps(chWGHT);/* convert to caps */

	/* ensure that chLabel is all caps */
	strcpy( chELEM, chLabel );
	caps(chELEM);/* convert to caps */

	/* now see if volume or radius weighting */
	if( strcmp(chWGHT,"RADIUS") == 0 )
	{
		lgVol = FALSE;
	}
	else if( strcmp(chWGHT,"VOLUME") == 0 )
	{
		lgVol = TRUE;
	}
	else
	{
		fprintf( ioQQQ, " cdTemp: chWeight=%6.6s makes no sense to me, the options are RADIUS and VOLUME.\n", 
		  chWeight );
		*TeMean = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(1);
	}

	if( IonStage == 0 )
	{
		/* special case, element is 21cm mean */
		if( strcmp(chELEM,"21CM") == 0 )
		{
			if( lgVol )
			{
				/* this is mean of inverse temperature weighted over volume */
				if( IonMeans.HarMeanTempVolume[1] > SMALLFLOAT )
				{
					*TeMean = IonMeans.HarMeanTempVolume[0]/IonMeans.HarMeanTempVolume[1];
				}
				else
				{
					*TeMean = 0.;
				}
				if( *TeMean > SMALLFLOAT )
				{
					*TeMean = 1./ *TeMean;
				}
				else
				{
					*TeMean = 0.;
				}
			}
			else
			{
				/* this is mean of inverse temperature weighted over radius */
				if( IonMeans.HarMeanTempRadius[1] > SMALLFLOAT )
				{
					*TeMean = IonMeans.HarMeanTempRadius[0]/IonMeans.HarMeanTempRadius[1];
				}
				else
				{
					*TeMean = 0.;
				}
				if( *TeMean > SMALLFLOAT )
				{
					*TeMean = 1./ *TeMean;
				}
				else
				{
					*TeMean = 0.;
				}
			}
		}
		/* special case, mean wrt H_2 */
		else if( strcmp(chELEM,"H2  ") == 0 )
		{
			if( lgVol )
			{
				/* mean temp with H2 over volume */
				if( IonMeans.H2MeanTempVolume[1] > SMALLFLOAT )
				{
					*TeMean = IonMeans.H2MeanTempVolume[0] / IonMeans.H2MeanTempVolume[1];
				}
				else
				{
					*TeMean = 0.;
				}
			}
			else
			{
				/* mean temp with H2 over radius */
				if( IonMeans.H2MeanTempRadius[1] > SMALLFLOAT )
				{
					*TeMean = IonMeans.H2MeanTempRadius[0] / IonMeans.H2MeanTempRadius[1];
				}
				else
				{
					*TeMean = 0.;
				}
			}
		}
		else
		{
			fprintf( ioQQQ, " cdTemp called with ion=0 and unknown quantity, =%4.4s\n", 
			  chLabel );
			fprintf( ioQQQ, " I know about 21CM and H2__ \n");
			
#			ifdef DEBUG_FUN
			fputs( " <->cdTemp()\n", debug_fp );
#			endif
			return(1);
		}

#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		/* say things are ok */
		return(0);
	}

	/* find which element this is */
	nelem = 0;
	while( nelem < LIMELM && 
		strcmp(chELEM,elementnames.chElementNameShort[nelem]) != 0 )
	{
		++nelem;
	}
	/* after this loop nelem is atomic number of element, H is 1 */

	/* if element not recognized and above loop falls through, nelem is LIMELM+1 
	 * nelem counter is on physical scale, H = 1 Zn = 30 */
	if( nelem >= LIMELM )
	{
		fprintf( ioQQQ, " cdTemp called with unknown element chLabel, =%4.4s\n", 
		  chLabel );
		
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(1);
	}

	/* sanity check - does this ionization stage exist? 
	 * both IonStage on physical scale, nelem on c */

	/* ipaaa will be used as pointer within the aaa array that contains average values,
	 * done this way to prevent lint from falsing in acess to aaa array */
	ipaaa = IonStage - 1;

	/*if( IonStage > nelem + 1 || IonStage < 1 || IonStage > LIMELM+1 )*/
	if( ipaaa > nelem+1 || ipaaa < 0 || ipaaa > LIMELM )
	{
		fprintf( ioQQQ, " cdTemp asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 
		  IonStage, chLabel );
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(1);
	}

	/* get either volume or radius average, aaa is filled in from 0 */
	if( lgVol )
	{
		/* aaa is dim'd LIMELM+1 so largest arg is LIMELM */
		VolMean('t', nelem+1,&ip,aaa,FALSE);
		*TeMean = pow(10.f,aaa[ipaaa]);
	}
	else
	{
		RadMean('t', nelem+1,&ip,aaa,FALSE);
		*TeMean = pow(10.f,aaa[ipaaa]);
	}

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

/*************************************************************************
 *
 *cdRead routine to read in command lines when cloudy used as subroutine 
 *
 ************************************************************************/

int cdRead(char *chInputLine)	/* the string containing the commands */
		/* returns the number of lines available in command stack 
		 * this is limit to how many more commands can be read */
{
	char *chEOL , /* will be used to seach for end of line symbols */
		chCARD[81],
		chLocal[81];

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

	/* this is set FALSE when code loaded, set TRUE when cdInit called,
	 * this is check that cdInit was called first */
	if( !lgcdInitCalled )
	{
		printf(" cdInit was not called first - this must be the first call.\n");
		puts( "[Stop in cdRead]" );
		cdEXIT(1);
	}

	/* totally ignore any card starting with a #, *, space, //, or % */
	if( (chInputLine[0]=='#') || 
		 (chInputLine[0]=='*') || 
		 (chInputLine[0]=='%') ||
		 (chInputLine[0]=='\n')||
		 (chInputLine[0]==' ') || 
		 (strncmp(chInputLine,"//", 2 )==0 ) )
	{
		/* return value is number of lines that can still be stuffed in */
		return(NKRD - input.nSave);
	}

	/***************************************************************************
	* validate a location to store this line image, then store the version     *
	* that has been truncated from special end of line characters              *
	* stored image will have <=80 valid characters                             *
	***************************************************************************/

	/* this will now point to the next free slot in the line image save array 
	 * this is where we will stuff this line image */
	++input.nSave;

	if( input.nSave >= NKRD )
	{
		/* too many input commands were entered - bail */
		fprintf( ioQQQ, 
			" Too many line images entered to cdRead.  The limit is%5d\n", 
		  NKRD );
		fprintf( ioQQQ, 
			" The limit to the number of allowed input lines is%4d.  This limit was exceeded.  Sorry.\n", 
		  NKRD );
		fprintf( ioQQQ, 
			" This limit is set by the variable NKRD which appears in parameter statements throughout the code.\n" );
		fprintf( ioQQQ, 
			" Increase it everywhere it appears.\n" );
		puts( "[Stop in cdRead]" );
		cdEXIT(1);
	}

	/* now kill any part of line image after special end of line character,
	 * this stops info user wants ignored from entering after here */
	if( (chEOL = strchr(chInputLine , '\n' ) ) !=NULL )
	{
		*chEOL = '\0';
	}
	if( (chEOL = strchr(chInputLine , '%' ) ) !=NULL )
	{
		*chEOL = '\0';
	}
	if( (chEOL = strchr(chInputLine , ';' ) ) !=NULL )
	{
		*chEOL = '\0';
	}
	if( (chEOL = strstr(chInputLine , "//" ) ) !=NULL )
	{
		*chEOL = '\0';
	}

	/* check that there are less than 80 characters to the null character 
	 * first find the end of string mark */
	chEOL = strchr(chInputLine , '\0' );

	/* do something if no end of string, or if string longer than 
	 * LINELENGTH characters */
	if( (chEOL==NULL) || (chEOL - chInputLine)>LINELENGTH )
	{
		/* remember, in C [LINELENGTH] is the 81th character */
		chInputLine[LINELENGTH] = '\0' ;
		/* this will generate several comments later on, so we know we had too long a line */
		lgCardSavTooLong = TRUE;
	}

	/* now do it again, since we now want to make sure that there is a trailing space
	 * if the line is shorter than 80 char */
	chEOL = strchr(chInputLine , '\0' );
	assert( chEOL !=NULL );

	/* pad with a space if short enough,
	 * if not short enough for this to be done, then up to user to create correct input */
	strcpy( chLocal , chInputLine );
	if( chEOL-chInputLine <80 )
	{
		strcat( chLocal , " " );
	}

	/* copy string over the chCARD, then convert this to caps to chk vary*/
	strcpy( chCARD, chLocal );
	caps(chCARD);/* convert to caps */

	/* save string in master array for later use in readr */
	strcpy( input.chCardSav[input.nSave], chLocal );

	/* check whether this is a trace command, turn on printout if so */
	if( strncmp(chCARD,"TRACE",5) == 0 )
	{
		trace.lgTraceInput = TRUE;
	}

	/* print input lines if trace specified */
	if( trace.lgTraceInput )
	{
		fprintf( ioQQQ,"cdRead=%s=\n",input.chCardSav[input.nSave] );
	}

	/* now check whether VARY is specified */
	if( lgMatch("VARY",chCARD) )
	{
		/* a vary flag was on the line image */
		vary.lgVaryOn = TRUE;
	}

	/* now check whether line is "no vary" command */
	if( strncmp(chCARD,"NO VARY",7) == 0 )
	{
		vary.lgNoVary = TRUE;
	}

	if( strncmp(chCARD,"OPTI",4) == 0 )
	{
		/* optimize command read in */
		optimr.lgOptimr = TRUE;
	}

#	ifdef DEBUG_FUN
	fputs( " <->cdRead()\n", debug_fp );
#	endif
	return( NKRD - input.nSave );
}


/* print line wavelengths in Angstroms in the standard format - just a wrapper */
void cdPrtWL( FILE *io , float wl )
{
	prt_wl( io , wl );
	return;

}