/* 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 */
/*ParseCommands main command line parser, decode command, then call other routines to read */
#include "cddefines.h"
#include "physconst.h"
#include "optimize.h"
#include "doppvel.h"
#include "abund.h"
#include "geometry.h"
#include "dense.h"
#include "plot.h"
#include "ionbal.h"
#include "rfield.h"
#include "dynamics.h"
#include "magnetic.h"
#include "trace.h"
#include "rt.h"
#include "thermal.h"
#include "opacity.h"
#include "atomfeii.h"
#include "called.h"
#include "extinc.h"
#include "wind.h"
#include "hextra.h"
#include "iterations.h"
#include "radius.h"
#include "input.h" 
#include "assertresults.h"
#include "phycon.h"
#include "fudgec.h"
#include "version.h"
#include "neutrn.h"
#include "converge.h"
#include "tabden.h"
#include "fabden.h"
#include "parse.h"
/*lint -e661 creation of  out of bound pointer */

void ParseCommands(void)
{
	char chCard[INPUT_LINE_LENGTH], 
	  chKey2[3], 
	  chKey3[4], 
	  chKey4[5], 
	  chKey5[6];

	int lgDSet, 
	  lgEOF, 
	  lgEOL, 
	  lgNu2;
	int lgStop ,
		lgStop_not_enough_info;

	char chK1;

	long int i, 
	  j, 
	  nqh;
	long nInitFile=0;/* used to count number of init files read in */

	double var;

	float a, 
	  ar1, 
	  teset, 
	  z;

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

	/* following says abundances are solar  */
	abund.lgAbnSolar = TRUE;

	/* there are no plots desired yet */
	plotCom.lgPlotON = FALSE;
	plotCom.nplot = 0;

	/* this flag remembers whether grains have ever been turned on.  It is passed
	 * to routine ParseAbundances - there grains will be turned on with commands
	 * such as abundances ism, unless grains were previously set 
	 * with a grains command.  this way abundances will not clobber explictly set
	 * grains. */
	lgDSet = FALSE;

	radius.Radius = 0.;
	radius.rinner = 0.;
	nqh = 0;
	rfield.nspec = 0;

	/* will be set true if no buffering command entered */
	input.lgSetNoBuffering = FALSE;

	/* initialize some code variables in case assert command used in input stream */
	InitAssertResults();

	for( i=0; i < LIMSPC; i++ )
	{
		strcpy( rfield.chRSpec[i], "UNKN" );
	}
	optimize.nparm = 0;

	/* this is an option to turn on trace printout on the nth
	 * call from the optimizer */
	if( optimize.lgTrOpt )
	{
		/* nTrOpt was set with "optimize trace" command,
		 * is iteration to turn on trace */
		if( optimize.nTrOpt == optimize.nOptimiz )
		{
			trace.lgTrace = TRUE;
			called.lgTalk = TRUE;
			trace.lgTrOvrd = TRUE;
			fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
		}
	}

	/* say this is a beta version if we are talking and it is the truth */
	if( version.nBetaVer > 0 && called.lgTalk )
	{
		fprintf( ioQQQ, 
			"\n                             This is %s (beta %ld) of Cloudy, and is intended for testing only.\n\n", 
			version.chVersion, version.nBetaVer );
	}

	if( called.lgTalk )
	{
		/* this code prints pretty lines at top of output box */
		for( i=0; i<57 ; i++ ) 
		{
			fprintf(ioQQQ,"%c",' ');
		}
		fprintf( ioQQQ, " Cloudy %7.7s\n\n", version.chVersion);

		for( i=0; i<23 ; i++ ) 
		{
			fprintf(ioQQQ,"%c",' ');
		}

		/* now print box and date of version, before printing commands */
		fprintf( ioQQQ, "**************************************");
		fprintf( ioQQQ, "%7.7s", version.chDate);
		fprintf( ioQQQ, "**************************************\n");

		for( i=0; i<23 ; i++ ) 
		{
			fprintf(ioQQQ,"%c",' ');
		}
		fprintf( ioQQQ, "*");

		for( i=0; i<81 ; i++ ) 
		{
			fprintf(ioQQQ,"%c",' ');
		}
		fprintf( ioQQQ, "*\n");
	}

	/* read in commands and print them   */

	/* following signals which read is in progress,
	 * 1 is forward read, of input commands
	 * -1 is reverse read from bottom of line array, of cloudy.ini file */
	input.iReadWay = 1;

	/* initialize array reader, this sub does nothing but set
	 * initial value of a variable, depending on value of iReadWay */
	input_init();

	/* read until eof or blank line, then rturn control back to main program */

	/*
	 * input_readarray puts the next stored line image into chCard
	 * it will be <=INPUT_LINE_LENGTH char long, with eol at end
	 * it will be in mixed upper and lower case
	 * lgEOF is set true if we hit eof and no more lines present
	 * code moving over to info contained in structures,
	 *
	 * input.chOrgCard == original image of command line
	 * input.chCARDCAPS == original image converted to caps
	 */
	input_readarray(chCard,&lgEOF);

	/* line starting with blank is eof, this checks we are not at eof */
	while( !lgEOF && chCard[0] != ' ' )
	{
		/* echo the line before we cap it, but only if it does
		 * not contain the keyword HIDE */
		/* >>chng 04 jan 21, add HIDE option, mostly for print quiet command */
		if( called.lgTalk && !lgMatch("HIDE",input.chCARDCAPS) )
		{
			fprintf( ioQQQ, "                       * ");
			i=0;
			while( chCard[i]!='\0' )
			{
				fprintf(ioQQQ,"%c",chCard[i]);
				++i;
			}

			while( i<80 )
			{
				fprintf(ioQQQ,"%c",' ');
				++i;
			}
			fprintf(ioQQQ,"*\n");
		}

		/* change chCard to all caps */
		caps( chCard );
		
		/* now set up several partial keys */
		chK1 = chCard[0];

		/* first two characters into chKey2 */
		strncpy( chKey2 , chCard , 2 );
		chKey2[2] = '\0';

		/* first three characters into chKey3 */
		strncpy( chKey3 , chCard , 3 );
		chKey3[3] = '\0';

		/* first four characters into chKey4 */
		strncpy( chKey4 , chCard , 4 );
		chKey4[4] = '\0';

		/* first four characters into chKey4 */
		strncpy( chKey5 , chCard , 5 );
		chKey5[5] = '\0';

		/* check whether VARY is on line */
		if( lgMatch("VARY",chCard) )
		{
			optimize.lgVarOn = TRUE;
			if( optimize.nparm + 1 > LIMPAR )
			{
				fprintf( ioQQQ, " Too many VARY lines entered; the limit is%4ld\n", 
				  LIMPAR );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else
		{
			optimize.lgVarOn = FALSE;
		}
 
		/*if( strncmp(chCard,"C ",2) == 0 )*/
		if( chCard[0]=='C' && (chCard[1]==' ' || chCard[1]== 0) )
		{
			/* this is null test since lines starting with "C " are comments,
			 * the char after the c can be a space or a newline */
			i = 1;
		}

		/* start to look for keywords */
		else if( strcmp(chKey4,"ABSO") == 0 )
		{
			/* enter luminosity in absolute magnitudes, in reads2 */
			ParseAbsMag(chCard,&nqh);
		}

		else if( strcmp(chKey3,"AGE") == 0 )
		{
			/* enter age of cloud so we can check for time-steady reads2 */
			ParseAge(chCard);
		}

		else if( strcmp(chKey4,"AGN ") == 0 )
		{
			/* enter generic style agn continuum, in reads2 */
			ParseAgn(chCard);
		}

		else if( strcmp(chKey4,"ABUN") == 0 )
		{
			/* chemical abundances, readsun, lgDSet is set true if grains command
			 * ever entered, so that abundances command does not override grains command */
			ParseAbundances(chCard,lgDSet);
			/* abundances no longer solar */
			abund.lgAbnSolar = FALSE;
		}

		else if( strcmp(chKey4,"APER") == 0 )
		{
			/* aperture command to simulate pencil beam or long slit */

			/* if the "BEAM" or "SLIT" option is specified then only part 
			 * of the geometry is observed, and intensities
			 * should not be weighted by r^2.  There are two limiting cases, SLIT in which
			 * the slit is longer than the diameter of the nebula and the contibution to the
			 * detected luminosity goes as r^1, and BEAM when the contribution is r^0, 
			 * the same as plane parallel 
			 * 
			 * default value of geometry.iEmissPower is 2 (set in zero.c) for full geometry  
			 */
			if( lgMatch("SLIT",chCard) )
			{
				/* long slit is case where slit is longer than diameter, so emissivity contributes
				 * r^1 to the observed luminosity */
				geometry.iEmissPower = 1;
			}
			else if( lgMatch("BEAM",chCard) )
			{
				/* long slit is case where slit is longer than diameter, so emissivity contributes
				 * r^1 to the observed luminosity */
				/* this is the default or SHORT case, where r^0 is dependence */
				geometry.iEmissPower = 0;
			}
			else
			{
				fprintf( ioQQQ, " One of the keywords SLIT or BEAM must appear.\n" );
				fprintf( ioQQQ, " Sorry.\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( strcmp(chKey4,"ASSE") == 0 )
		{
			/* assert that code predicts certain results, in assertresults.h */
			ParseAssertResults();
		}

		else if( strcmp(chKey4,"ATOM") == 0 )
		{
			if( lgMatch("FEII",chCard) )
			{
				/* parse the atom feii command - routine in FeIILevelPops.c */
				ParseAtomFeII(chCard);
			}

			else if( lgMatch("H-LI",chCard) )
			{
				/* parse the atom h-like command */
				ParseAtomHLike(chCard);
			}

			else if( lgMatch("HE-L",chCard) )
			{
				/* stuff to do with helium iso-sequence, in reads2 */
				ParseAtomHeLike(chCard);
			}

			else if( lgMatch("ROTO",chCard) || lgMatch(" CO ",chCard) )
			{
				/* parameters for the rigid rotor CO molecules (yes, its not an atom)
				 * command is ATOM ROTOR or ATOM CO */
				ParseAtomCO(chCard);
			}

			else if( lgMatch(" H2 ",chCard) )
			{
				/* parameters for the rigid rotor CO molecules (yes, its not an atom)
				 * command is ATOM ROTOR, routine is stored in parse_atomh2 */
				ParseAtomH2(chCard);
			}

			else
			{
				fprintf( ioQQQ, " I could not recognize a keyword on this atom command.\n");
				fprintf( ioQQQ, " The available keys are FeII and H-Like.\n");
				fprintf( ioQQQ, " Sorry.\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( strcmp(chKey4,"BACK") == 0 )
		{
			/* cosmic background, in parse_backgrd */
			ParseBackgrd(&nqh,chCard,&ar1);
		}

		else if( strcmp(chKey4,"BLAC") == 0 )
		{
			/* black body, in reads2 */
			ParseBlackbody(chCard,&nqh,&ar1);

			/* vary option */
			if( optimize.lgVarOn )
			{
				/* no luminosity options on vary */
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "BLACKbody=%f" );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/* log of temp stored here  */
				/* >>chng 02 feb 11, log was not here, don't know what happened to it */
				optimize.vparm[0][optimize.nparm] = (float)log10(rfield.slope[rfield.nspec-1]);
				/* the increment in the first steps away from the original value */
				optimize.vincr[optimize.nparm] = 0.5;
				optimize.nparm += 1;
			}
		}

		else if( strcmp(chKey4,"BREM") == 0 )
		{
			/* brems continuum from central object */
			strcpy( rfield.chSpType[rfield.nspec], "BREMS" );
			i = 5;
			rfield.slope[rfield.nspec] = 
				(float)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}

			/* temperature is interpreted as log if <=10, linear otherwise*/
			if( rfield.slope[rfield.nspec] <= 10. )
			{
				rfield.slope[rfield.nspec] =  
					pow(10.,rfield.slope[rfield.nspec]);
			}
			rfield.cutoff[0][rfield.nspec] = 0.;

			/* option for vary keyword */
			if( optimize.lgVarOn )
			{
				/* only one parameter */
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "BREMS, T=%f" );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/* log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] =  (float)rfield.slope[rfield.nspec];
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
			++rfield.nspec;
			if( rfield.nspec >= LIMSPC )
			{
				/* too many continua were entered */
				fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( strcmp(chKey4,"CASE") == 0 )
		/* do hydrogen case b */
		{
			/* set flag saying we are doing case b */
			opac.lgCaseB = TRUE;

			/* scan in an optional optical depth in lya */
			i = 5;
			opac.tlamin = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				/* set tau-La to 10**9 if not specified 
				opac.tlamin = 1e9;*/
				/* >>chng 01 sep 24, change from 1e9 default to 1e5, to get
				 * more realistic conditions in HII region.  Very large tau
				 * caused extreme Lya behavior */
				opac.tlamin = 1e5;
			}
			else
			{
				opac.tlamin = (float)pow(10.f,opac.tlamin);
			}

			/* Hummer and Storey case B; no collisions from n=1, 2 (usually in) */
			if( lgMatch("HUMM",chCard) )
			{
				opac.lgCaseB_HummerStorey = TRUE;
			}

			/* the NO PHOTOIONIZATION option, turns off excited state photoionization */
			if( lgMatch("O PH",chCard) )
			{
				opac.lgCaseB_no_photo = TRUE;
			}

			/* the NO PDEST option, turns off line dest by background opacities */
			if( lgMatch("O PDE",chCard) )
			{
				/* turn off destruction off line dest by background opacities */
				opac.lgCaseB_no_pdest = TRUE;
			}
		}

		else if( strcmp(chKey4,"CEXT") == 0 )
		{
			/* add "extra" cooling, power law temp dependence */
			thermal.lgCExtraOn = TRUE;
			i = 5;
			thermal.CoolExtra = (float)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			thermal.cextpw = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		}

		else if( strcmp(chKey4,"COMP") == 0 )
		{
			/* compile ascii version of stellar atmosphere continua in volk */
			ParseCompile(chCard);
		}

		else if( strcmp(chKey4,"CONS") == 0 )
		{
			/* constant temperature, pressure, density, or gas pressure
			 * in readsun */
			ParseConstant(chCard);
		}

		else if( strcmp(chKey4,"CORO") == 0 )
		{
			/* coronal equilibrium; set constant temperature to number on line
			 *  in readsun */
			ParseCoronal(chCard,&nqh,&ar1);
		}

		else if( strcmp(chKey4,"COSM") == 0 )
		{
			/* cosmic ray heating, log of density of rel. electrons */
			i = 5;
			hextra.cryden = (float)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			if( lgEOL )
			{
				if( lgMatch("BACK",chCard) )
				{
					/* >>chng 04 mar 10, background is now 
					 * >>refer	cr	ion	Williams, J.P., Bergin, E.A., Caseli, P., Myers, P.C., & Plume, R. 1998, ApJ, 503, 689 */
					/* galactic background cosmic ray density to produce
					 * secondary ionization rate quoted by Tielens and Hollenbach */
					/* hextra.cryden = 2e-9f;*/
					/* >>chng 99 jun 24, slight change to value
					 * quoted by 
					 * >>refer	cosmic ray	ionization rate	McKee, C.M., 1999, astro-ph 9901370
					 * this will produce a total
					 * secondary ionization rate of 2.5e-17 s^-1, as tested in 
					 * tsuite secondary.in.  If each ionization produces 2.4 eV of heat,
					 * the background heating rate should be 9.6e-29 * n*/
					/* >>chng 00 nov 28, changed density to 4.9e-9 to reproduce TH85a
					 * when photoionization is turned off. 
					 >>refer	cosmic ray	ionization rate	Tielens, A.G.G.M., & Hollenbach, D., 1998, ApJ, 291, 722
					 */
					/* hextra.cryden = 7.07e-9f;*/
					/* this value reproduces the TH cr ionization rate when the factor
					 * of 0.46 is included.  This will directly go onto the h ionization rate
					 * without the factor of 0.46 there.  this is necessary for the more
					 * general case where cr ionization is actually self-consistently determined
					 * from rate hot electrons injected into the plasma */
					/*hextra.cryden = 2.25e-9f;*/
					hextra.cryden = hextra.background;
					hextra.crtemp = 2.6e9f;
				}
				else
				{
					fprintf(ioQQQ,
						" Either the log of the cosmic ray density, or the keyword BACKGROUND, must appear.\n");
					NoNumb(chCard);
				}
			}
			else
			{
				/*  optional power law density  */
				hextra.crpowr = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
				/*  option to specify a temp for non-rel electrons */
				hextra.crtemp = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
				if( lgEOL )
				{
					/* relavistic limit (Balbus and McKee) */
					hextra.crtemp = 2.6e9;
				}
				else
				{
					var = pow(10.f,hextra.crtemp);
					hextra.crtemp = (float)MIN2(var,2.6e9);
				}
			}
			/* this is current cosmic ray density divided by background */
			hextra.cryden_ov_background = hextra.cryden / hextra.background;
		}

		else if( strcmp(chKey4,"COVE") == 0 )
		{
			/* covering factor for gas */
			i = 5;
			/* The geometric covering factor accounts for how much of 4\pi is
			 * covered by gas, and so linearly multiplies the predictd intensities */
			geometry.covgeo = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

			if( lgEOL )
			{
				NoNumb(chCard);
			}

			/* if negative then log, convert to linear */
			if( geometry.covgeo <= 0. )
			{
				geometry.covgeo = (float)pow(10.f,geometry.covgeo);
			}

			if( geometry.covgeo > 1. )
			{
				fprintf( ioQQQ, " A covering factor greater than 1 makes no physical sense.  Sorry.\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* radiative transfer covering factor will be equal to the geometric one */
			geometry.covrt = geometry.covgeo;
		}

		else if( strcmp(chKey4,"CRAS") == 0 )
		{
			/* any of several tests to cause the code to crash */
			CrashDo(chCard);
		}

		else if( strcmp(chKey4,"CYLI") == 0 )
		{
			/* height of cylinder in cm */
			i = 5;
			radius.lgCylnOn = TRUE;
			radius.CylindHigh = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			if( lgEOL )
			{
				NoNumb(chCard);
			}
		}

		else if( strcmp(chKey4,"DIEL") == 0 )
		{
			/* change various factors for dielectronic recombination */
			if( lgMatch("KLUD",chCard) )
			{
				/* option to turn off guess of diel rec coef for 3rd row elements */
				i = 3;
				/* this is first call, no number, lgEOL true but zero returned, the intended effect*/
				ionbal.GuessDiel[0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
				for( j=1; j<4; ++j )
				{
					ionbal.GuessDiel[j] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
					/* here lgEOL means to use previous values */
					if( lgEOL )
						ionbal.GuessDiel[j] = ionbal.GuessDiel[j-1];
				}
			}

			else if( lgMatch("BURG"  , chCard) )
			{
				/* modify suppression of burgess dielectronic recombinations */
				if( lgMatch(" ON ",chCard) )
				{
					/* leave suppression on - this is default state */
					ionbal.lgSupDie[0] = TRUE;
				}

				else if( lgMatch(" OFF",chCard) )
				{
					/* turn suppression off */
					ionbal.lgSupDie[0] = FALSE;
				}

				else
				{
					fprintf( ioQQQ, " flag ON or OFF must appear.\n" );
					puts( "[Stop in ParseCommands]" );
					cdEXIT(EXIT_FAILURE);
				}
			}

			else if( lgMatch("NUSS"  , chCard) )
			{
				/* modify suppression of nussbaumer and storey dielectronic recomb */
				if( lgMatch(" ON ",chCard) )
				{
					/* turn suppression on */
					ionbal.lgSupDie[1] = TRUE;
				}
				else if( lgMatch(" OFF",chCard) )
				{
					/* leave suppression off - this is default state */
					ionbal.lgSupDie[1] = FALSE;
				}
				else
				{
					fprintf( ioQQQ, " flag ON or OFF must appear.\n" );
					puts( "[Stop in ParseCommands]" );
					cdEXIT(EXIT_FAILURE);
				}
			}

			else if( lgMatch("STEVE" , chCard ) ) 
			{
				/* turn on steve option, must be one number, the log of the coefficient,
				 * which will be multiplied by the charge of the ion before recombination */
				if( lgMatch( " OFF" , chCard ) )
				{
					ionbal.lg_guess_coef = FALSE;
				}
				else
				{
					ionbal.lg_guess_coef = TRUE;
				}
				/* option to add log normal noise */
				if( lgMatch("NOISE" , chCard ) ) 
				{
					i = 3;
					ionbal.guess_noise = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
					/* Seed the random-number generator with current time so that
					 * the numbers will be different every time we run. */
					srand( (unsigned)time( NULL ) );
				}
				else
				{
					ionbal.guess_noise = 0.;
				}
			}
			else
			{
				fprintf( ioQQQ, " key KLUDge, BURGess, STEVE, or NUSSbaumer must appear.\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( strcmp(chKey4,"DIFF") == 0 )
		{
			/* set method of transferring diffuse fields,
			 * default is outward only, cdDffTrns label "OU2", set in zero.c */
			if( lgMatch(" OTS",chCard) )
			{
				if( lgMatch("SIMP",chCard) )
				{
					/* this is simple ots, which never allows photons to escape */
					strcpy( rfield.chDffTrns, "OSS" );
				}
				else
				{
					/* this is regular ots, which allows photons to escape */
					strcpy( rfield.chDffTrns, "OTS" );
				}
				rfield.lgOutOnly = FALSE;
			}
			else if( lgMatch(" OUT",chCard) )
			{
				rfield.lgOutOnly = TRUE;
				i = 4;
				j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
				if( lgEOL )
				{
					/* this is the default set in zero */
					strcpy( rfield.chDffTrns, "OU2" );
				}
				else
				{
					if( j > 0 && j < 10 )
					{
						sprintf( rfield.chDffTrns, "OU%1ld", j );
					}
					else
					{
						fprintf( ioQQQ, " must be between 1 and 9 \n" );
						puts( "[Stop in ParseCommands]" );
						cdEXIT(EXIT_FAILURE);
					}
				}
			}

			else
			{
				fprintf( ioQQQ, " There should have been OUTward or OTS on this line.  Sorry.\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( strcmp(chKey4,"DIST") == 0 )
		{
			/* distance to the object */
			i = 5;
			radius.distance = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			/* default is for quantity to be log of distance, linear keyword
			 * overrides this - is LINEAR is not on line then exp */
			if( !lgMatch("LINE",chCard ) )
			{
				radius.distance = pow(10., radius.distance );
			}
			/* default is radius in cm - if parsecs appears then must 
			 * convert to cm */
			if( lgMatch("PARS",chCard ) )
			{
				radius.distance *= PARSEC;
			}
		}

		else if( strcmp(chKey4,"DLAW") == 0 )
		{
			/* either use fabden routine, or read in table of points */
			ParseDLaw(chCard);
		}

		else if( strcmp(chKey4,"DOUB") == 0 )
		{
			/* double optical depth scale after each iteration */
			rt.DoubleTau = 2.;
		}

		else if( strcmp(chKey4,"DRIV") == 0 )
		{
			/* call one of several drivers, readsun */
			ParseDriveCmnd(chCard);
		}

		else if( strcmp(chKey4,"EDEN") == 0 )
		{
			/* option to add "extra" electrons, to test compton limit
			 *  for very low T(star) - option is log of eden */
			i = 5;
			dense.EdenExtra = (float)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			/* warn that this model is meaningless */
			phycon.lgPhysOK = FALSE;
		}

		else if( strcmp(chKey4,"ELEM") == 0 )
		{
			/* element command;
			 * scale or abundance options, to change abundance of specific element
			 * read option to change order of elements
			 * in reads2.f */
			ParseElement(chCard);
		}

		else if( strcmp(chKey4,"ENER") == 0 )
		{
			/* energy density (degrees K) of this continuum source */
			i = 5;
			nqh += 1;
			if( nqh > LIMSPC )
			{
				/* too many continua were entered */
				fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* silly, but calms down the lint */
			ASSERT( nqh <= LIMSPC );
			strcpy( rfield.chRSpec[nqh-1], "SQCM" );
			teset = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			/* set radius to very large value if not already set */
			/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
			if( !radius.lgRadiusKnown )
			{
				/* set radius to large value */
				ar1 = (float)radius.rdfalt;
				radius.Radius = pow(10.f,radius.rdfalt);
			}

			/* convert temp to log, recognize linear option */
			if( lgMatch("LINE",chCard) || teset > 10. )
			{
				/* option for linear temperature, must store log */
				teset = (float)log10(teset);
			}

			if( teset > 5. )
			{
				fprintf( ioQQQ, " This intensity may be too large.  The code may crash due to overflow.  Was log intended?\n" );
			}

			/* teset is not log of temp, now get log of total luminosity */
			strcpy( rfield.chSpNorm[nqh-1], "LUMI" );

			/* full range of continuum will be used */
			rfield.range[0][nqh-1] = rfield.emm;
			rfield.range[1][nqh-1] = rfield.egamry;
			rfield.totpow[nqh-1] = (4.*teset - 4.2464476 + 0.60206);

			/* vary option */
			if( optimize.lgVarOn )
			{
				strcpy( optimize.chVarFmt[optimize.nparm], "ENERGY DENSITY %f log " );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				optimize.vparm[0][optimize.nparm] = (float)log10(rfield.totpow[nqh-1]);
				optimize.vincr[optimize.nparm] = 0.1f;
				optimize.nvarxt[optimize.nparm] = 1;
				optimize.nparm += 1;
			}
		}

		else if( strcmp(chKey4,"EXTI") == 0 )
		{
			/* extinguish ionizing continuum by absorbing column AFTER
			 * setting luminosity or Q(H).  First number is the column
			 * density (log), second number is leakage (def=0%)
			 * last number is lowest energy (ryd), last two may be omitted
			 * from right to left 
			 * 
			 * extinction is actually done in extin, which is called by ContSetIntensity */

			i = 5;
			extinc.excolm = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			/* >>chng 01 dec 19, add linear option */
			/* default is for the number to be the log of the column. 
			 * there is a linear option for the column or optical depth,
			 * if linear does not occur then log, and convert to linear */
			if( !lgMatch("LINE" , chCard ) )
			{
				/* >>chng 03 jan 02, make comment if first value seems too large */
				if( extinc.excolm>35. )
				{
					fprintf(ioQQQ,
						" The first parameter on this command line is the log of either the column density or optical depth.\n");
					fprintf(ioQQQ,
						" The value seems pretty big to me - please check it.\n");
					/* flush it since we will probably crash */
					fflush(ioQQQ);
				}
				extinc.excolm = (float)pow(10.f,extinc.excolm);
			}

			/* option to set leakage - default is 0. */
			extinc.exleak = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* optional leakage is zero */
			if( lgEOL )
			{
				extinc.exleak = 0.;
			}
			/* negative leaks are logs */
			if( extinc.exleak < 0. )
			{
				extinc.exleak = (float)pow(10.f,extinc.exleak);
			}
			if( extinc.exleak > 1. )
			{
				/* but leaks greater than 1 are not allowed */
				if( called.lgTalk )
				{
					fprintf( ioQQQ, " A leakage of%9.0f%% was entered - this must be less than 100%%\n", 
					  extinc.exleak*100. );
				}
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* option to set lowest energy for absorber */
			extinc.exlow = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				extinc.exlow = 0.99946f;
			}

			else
			{
				if( extinc.exlow <= 0. )
					extinc.exlow = (float)pow(10.f,extinc.exlow);
				if( extinc.exlow < 0.99946 )
				{
					fprintf( ioQQQ, " Energy less than 1 Ryd!!\n" );
				}
			}

			/* >>chng 01 dec 19, add optical depth option, to replace column density */
			if( lgMatch("OPTI" , chCard ) )
			{
				/* convert the optical depth into the proper column density */
				extinc.excolm /= (float)(extinc.cnst_col2optdepth*
					pow(extinc.exlow,extinc.cnst_power) );
			}

		}

		else if( strcmp(chKey4,"FAIL") == 0 )
		{
			/* reset number of temp failures allowed, default=20 */
			i = 5;

			/* save previous value */
			j = conv.LimFail;
			conv.LimFail = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}

			/* >>chng 01 mar 14, switch logic on maps */
			/* ' map' flag, make map when failure, default is no map,
			 * second check is so that no map does not trigger a map */
			if( lgMatch(" MAP",chCard) && !lgMatch(" NO ",chCard) )
			{
				conv.lgMap = TRUE;
			}

			/* complain if failures was increased above default */
			if( conv.LimFail > j )
			{
				fprintf( ioQQQ, " Cloudy should not have problems making this command necessary.\n" );
				fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this model.\n" );
			}

		}

		else if( strcmp(chKey4,"FILL") == 0 )
		{
			/* filling factor, power law exponent (default=1., 0.) */
			i = 5;
			a = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}

			if( a <= 0. )
			{
				/* number less than or equal to 0, must have been entered as log */
				geometry.FillFac = (float)pow(10.f,a);
			}
			else
			{
				/* number greater than zero, must have been the real thing */
				geometry.FillFac = a;
				if( geometry.FillFac > 1.0 )
				{
					if( called.lgTalk )
					{
						fprintf( ioQQQ, " Filling factor > 1, reset to 1\n" );
					}
					geometry.FillFac = 1.;
				}
			}
			geometry.fiscal = geometry.FillFac;

			/* option to have power law dependence, filpow set to 0 in zerologic */
			geometry.filpow = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

			/* vary option */
			if( optimize.lgVarOn )
			{
				strcpy( optimize.chVarFmt[optimize.nparm], "FILLING FACTOR= %f power=%f" );

				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				optimize.vparm[0][optimize.nparm] = (float)log10(geometry.FillFac);
				optimize.vincr[optimize.nparm] = 0.5;

				/* power law dependence here, but cannot be varied */
				optimize.vparm[1][optimize.nparm] = geometry.filpow;
				optimize.nvarxt[optimize.nparm] = 2;

				/* do not allow filling factor to go positive */
				optimize.varang[optimize.nparm][0] = -1e38f;
				optimize.varang[optimize.nparm][1] = 0.f;
				optimize.nparm += 1;
			}
		}

		else if( strcmp(chKey4,"FIRE") == 0 )
		{
			/* cosmic thermal background radiation, argument is redshift */
			i = 5;
			/* if no number on line then (float)FFmtRead returns z=0; i.e., now */
			z = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* in readsun */
			ParseFireBall(z,&nqh,&ar1);
		}

		else if( strcmp(chKey4,"FLUC") == 0 )
		{
			/* rapid density fluctuations, in readsun */
			ParseFluc(chCard);
		}

		else if( strcmp(chKey4,"F(NU") == 0 )
		{
			/* this is the specific flux at nu
			 *  following says F(nu) not nuF(nu) */
			lgNu2 = FALSE;
			/* in reads2 */
			ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2);
		}

		else if( strcmp(chKey4,"FORC") == 0 )
		{
			/* set temperature of first zone, but don't keep constant
			 *  log if < 10 */
			i = 5;
			thermal.ConstTemp = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}

			if( lgMatch(" LOG",chCard) || (thermal.ConstTemp <= 10. && 
			  !lgMatch("LINE",chCard)) )
			{
				thermal.ConstTemp = (float)pow(10.f,thermal.ConstTemp);
			}

			/* low energy bounds of continuum array do not permit too-low a temp */
			if( thermal.ConstTemp < 3. )
			{
				fprintf( ioQQQ, " TE reset to 3K: entered number too small.\n" );
				thermal.ConstTemp = 3.;
			}
		}

		else if( strcmp(chKey4,"FUDG") == 0 )
		{
			/* enter a fudge factor */
			i = 5;
			for( j=0; j < NFUDGC; j++ )
			{
				fudgec.fudgea[j] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
				if( !lgEOL )
					fudgec.nfudge = j+1; /* this will be fortran counting for now */
			}
		}

		else if( strcmp(chKey4,"GLOB") == 0 )
		{
			/* globule with density increasing inward
			 * parameters are outer density, radius of globule, and density power */
			ParseGlobule(chCard);
		}

		else if( (strcmp(chKey4,"GRAI") == 0 )||( strcmp(chKey4,"PGRA") == 0 ) )
		{
			/* read parameters dealing with grains, in reads2 */
			ParseGrain(chCard,&lgDSet);
		}

		else if( strcmp(chKey4,"HDEN") == 0 )
		{
			/* parse the hden command to set the hydrogen density, in reads2 */
			ParseHDEN(chCard);
		}

		else if( strcmp(chKey4,"HELI") == 0 )
		{
			/* stuff to do with helium, in reads2 */
			ParseAtomHeLike(chCard);
		}

		else if( strcmp(chKey4,"HEXT") == 0 )
		{
			/* "extra" heating rate, so that first= log(erg(cm-3, s-1))
			 * second optional number is scale radius, so that HXTOT = TurbHeat*SEXP(DEPTH/SCALE)
			 * if missing then constant heating.
			 * third option is depth from shielded face, to mimic irradiation from both sides*/
			i = 5;
			hextra.TurbHeat = (float)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			a = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				hextra.turrad = 0.;
			}
			else
			{
				hextra.turrad = (float)pow(10.f,a);
			}
			a = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				hextra.turback = 0.;
			}
			else
			{
				hextra.turback = (float)pow(10.f,a);
			}
		}

		else if( strcmp(chKey4,"HIGH") == 0 )
		{
			/* approach equilibrium from high te */
			thermal.lgTeHigh = TRUE;
		}

		else if( strncmp( chCard ,"HYDROGEN",8) == 0 )
		{
			fprintf(ioQQQ," This command has been replaced with the ATOM H-LIKE command.\n");
			fprintf(ioQQQ," I will parse the command for now, but may not in the future.\n");
			/* stuff to do with hydrogen atom, in reads2 */
			ParseAtomHLike(chCard);
		}

		else if( strcmp(chKey4,"ILLU") == 0 )
		{
			/* illumination angle */
			i = 5;
			geometry.AngleIllum = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			/* default is degrees, but radian key also there */
			if( lgMatch("RADI",chCard) )
			{
				/* entered as radians, convert to degrees first */
				geometry.AngleIllum /= (float)(PI/180.);
			}
			/* we now have degrees, make sure between 0 and 90 */
			if( geometry.AngleIllum < 0. || geometry.AngleIllum >= 90. )
			{
				fprintf( ioQQQ, " Angle of illumination must be between 0 and 90.\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* we really want 1. / COS( theta ) with theta in radians */
			geometry.AngleIllum = (float)(geometry.AngleIllum*PI/180.);
			geometry.AngleIllum = (float)(1./cos(geometry.AngleIllum));
		}

		else if( strcmp(chKey4,"INIT") == 0 )
		{

			/* read cloudy.ini initialization file
			 * need to read in file every time, since some vars
		     * get reset in zero - would be unsafe not to read in again */
			ParseInit(chCard);

			/* check that only one init file was in the input stream -
			 * we cannot now read more than one */
			++nInitFile;
			if( nInitFile > 1 )
			{
				fprintf( ioQQQ, 
					" This is the second init file, I can only handle one.\nSorry.\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* we will put the data for the ini file at the end of the array of
			 * line images, this is the increment for stuffing the lines in - negative */
			input.iReadWay = -1;

			/* initialize array reader, this sub does nothing but set
			 * initial value of a variable, depending on value of iReadWay */
			input_init();
		}

		else if( strcmp(chKey5,"INTEN") == 0 )
		{
			/* intensity of this continuum source */
			i = 5;
			if( nqh >= LIMSPC )
			{
				/* too many continua were entered */
				fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* silly, but calms down the lint */
			ASSERT( nqh < LIMSPC );
			strcpy( rfield.chRSpec[nqh], "SQCM" );
			rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}

			/* set radius to very large value if not already set */
			/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
			if( !radius.lgRadiusKnown )
			{
				/* set radius to large value */
				ar1 = (float)radius.rdfalt;
				radius.Radius = pow(10.f,radius.rdfalt);
			}
			if( lgMatch("LINE",chCard) )
			{
				/* silly, but calms down the lint */
				ASSERT( nqh < LIMSPC );
				/* option for linear input parameter */
				rfield.totpow[nqh] = log10(rfield.totpow[nqh]);
			}
			strcpy( rfield.chSpNorm[nqh], "LUMI" );
			/* ParseRangeOption in readsun */
			ParseRangeOption(nqh,chCard);

			/* vary option */
			if( optimize.lgVarOn )
			{
				/* create the format we will use to write the command */
				/* >>chng 03 apr 30, had used log for range option, but if one or other
				 * was 1 ryd, then became 0 in the log, which was misinterpreted as
				 * one of the continuum bounds */
				/*strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY%f log range %f %f" );*/
				strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY %f range %f %f " );
				/* array index for this command */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				optimize.vparm[0][optimize.nparm] = (float)rfield.totpow[nqh];
				optimize.vincr[optimize.nparm] = 0.5;
				/* range option, but cannot be varied */
				/*optimize.vparm[1][optimize.nparm] = (float)log10(rfield.range[0][nqh]);
				optimize.vparm[2][optimize.nparm] = (float)log10(rfield.range[1][nqh]);*/
				/* >>change to linear option */
				optimize.vparm[1][optimize.nparm] = (float)rfield.range[0][nqh];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[1][nqh];
				optimize.nvarxt[optimize.nparm] = 3;
				++optimize.nparm;
			}
			++nqh;
		}

		else if( strcmp(chKey5,"INTER") == 0 )
		{
			/* interpolate on input tables for continuum, set of power  laws used
			 * input ordered pairs nu( ryd or log(hz) ), log(fnu)
			 * additional lines begin CONTINUE
			 * first check that this is the one and only INTERP command
			 * in readsun */
			ParseInterp(chCard,&lgEOF);
		}

		else if( strcmp(chKey4,"IONI") == 0 )
		{
			/* inter ionization parameter U=Q/12 R*R N C;
			 * defined per hydrogen, not per electron (as before)
			 * radius must also be entered if spherical, but not needed if plane */
			ParseIonPar(&nqh,chCard,&ar1);
		}

		else if( strcmp(chKey4,"ITER") == 0 )
		{
			/* iterate to get optical depths and diffuse field properly */
			i = 5;
			iterations.itermx = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) - 1;
			iterations.itermx = MAX2(iterations.itermx,1);
			if( iterations.itermx > ITRDIM - 1 )
			{
				fprintf( ioQQQ, " No more than %ld can be performed because of vector sizes; reduce ITER or increase ITRDIM in code.\n", 
				  ITRDIM );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
			if( lgMatch("CONV",chCard) )
			{
				/* option to keep iterating until it converges, checks on convergence
				 * are done in update, and checked again in prtcomment */
				conv.lgAutoIt = TRUE;
				/* above would have been legitimite setting of ITERMX, set to default 10 */
				if( lgEOL )
				{
					iterations.itermx = 10 - 1;
				}
				a = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
				/* change convergence criteria, preset in zero */
				if( !lgEOL )
				{
					conv.autocv = a;
				}
			}
		}

		else if( strcmp(chKey4,"L(NU") == 0 )
		{
			/* this is the specific luminositiy at nu
			 * following says n(nu) not nuF(nu) */
			lgNu2 = FALSE;
			/* in reads2 */
			ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2);
		}

		else if( strcmp(chKey4,"LASE") == 0 )
		{
			/* mostly one frequency (a laser) to check gamma's,*/
			/* say the continuum type */
			strcpy( rfield.chSpType[rfield.nspec], "LASER" );

			/* scan off the laser's energy */
			i = 5;
			rfield.slope[rfield.nspec] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);

			/* negative energies are logs */
			if( rfield.slope[rfield.nspec] <= 0. )
			{
				rfield.slope[rfield.nspec] = 
				pow(10.,rfield.slope[rfield.nspec]);
			}
			if( lgEOL )
			{
				NoNumb(chCard);
			}

			/* next number is optional relative width of laser */
			rfield.cutoff[0][rfield.nspec] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				/* default width is 0.05 */
				rfield.cutoff[0][rfield.nspec] = 0.05;
			}

			/* increment counter making sure we still have space enough */
			++rfield.nspec;
			if( rfield.nspec >= LIMSPC )
			{
				/* too many continua were entered */
				fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( strcmp(chKey4,"LUMI") == 0 )
		{
			/* luminosity of this continuum source */
			i = 5;
			if( nqh >= LIMSPC )
			{
				/* too many continua were entered */
				fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* silly, but calms down the lint */
			ASSERT( nqh < LIMSPC );
			strcpy( rfield.chRSpec[nqh], "4 PI" );
			rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			if( lgMatch("LINE",chCard) )
			{
				/* silly, but calms down the lint */
				ASSERT( nqh < LIMSPC );
				/* option for linear input parameter */
				rfield.totpow[nqh] = log10(rfield.totpow[nqh]);
			}
			/* silly, but calms down the lint */
			ASSERT( nqh < LIMSPC );
			strcpy( rfield.chSpNorm[nqh], "LUMI" );
			if( lgMatch("SOLA",chCard) )
			{
				/* silly, but calms down the lint */
				ASSERT( nqh < LIMSPC );
				/* option to use log of total lumin in solar units */
				rfield.range[0][nqh] = rfield.emm;
				rfield.range[1][nqh] = rfield.egamry;
				rfield.totpow[nqh] += 33.5827f;
			}
			else
			{
				/* ParseRangeOption in readsun */
				ParseRangeOption(nqh,chCard);
			}
			/* vary option */
			if( optimize.lgVarOn )
			{
				/* >>chng 03 apr 30, had used log for range option, but if one or other
				 * was 1 ryd, then became 0 in the log, which was misinterpreted as
				 * one of the continuum bounds */
				strcpy( optimize.chVarFmt[optimize.nparm], "LUMINOSITY %f range %f %f " );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				optimize.vparm[0][optimize.nparm] = (float)rfield.totpow[nqh];
				optimize.vincr[optimize.nparm] = 0.5;
				/* range option in, but cannot be varied */
				optimize.vparm[1][optimize.nparm] = (float)rfield.range[0][nqh];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[1][nqh];
				optimize.nvarxt[optimize.nparm] = 3;
				optimize.nparm += 1;
			}
			++nqh;
		}

		else if( strcmp(chKey4,"MAGN") == 0 )
		{
			/* parse the magnetic field command */
			ParseMagnet( chCard );
		}

		else if( strcmp(chKey4,"MAP ") == 0 )
		{
			/* do cooling space map for specified zones, in reads2 */
			ParseMap(chCard);
		}

		else if( strcmp(chKey4,"META") == 0 )
		{
			/* read depletion for metals, all elements heavier than He
			 * in reads2 */
			ParseMetal(chCard);
			/* abundances no longer solar */
			abund.lgAbnSolar = FALSE;
		}

		else if( strcmp(chKey4,"NEUT") == 0 )
		{
			/* heating and ionization due to fast neutrons, arg is total luminosity
			 * in neutrons rel to boblmetric lumin; sec number is efficiency */
			i = 5;
			neutrn.frcneu = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			if( neutrn.frcneu > 0. )
			{
				neutrn.frcneu = (float)log10(neutrn.frcneu);
			}

			neutrn.lgNeutrnHeatOn = TRUE;
			neutrn.effneu = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				neutrn.effneu = 1.0;
			}
			else
			{
				if( neutrn.effneu <= 0. )
					neutrn.effneu = (float)pow(10.f,neutrn.effneu);
			}
		}

		else if( strcmp(chKey3,"NO ") == 0 )
		{
			/* don't do something, in readsun */
			ParseDont(chCard);
		}

		else if( strcmp(chKey4,"NORM") == 0 )
		{
			/* normalise lines to this rather than h-b, sec number is scale factor */
			ParseNorm(chCard);
		}

		else if( strcmp(chKey4,"NUF(") == 0 )
		{
			/* flux density of this continuum source, at optional frequency
			 *  expressed as product nu*f_nu */
			lgNu2 = TRUE;
			/* in reads2 */
			ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2);
		}

		else if( strcmp(chKey4,"NUL(") == 0 )
		{
			/* specific luminosity density of this continuum source, at opt freq
			 * expressed as product nu*f_nu */
			lgNu2 = TRUE;
			/* in reads2 */
			ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2);
		}

		else if( strcmp(chKey4,"OPTI") == 0 )
		{
			/* option to optimize the model by varying certain parameters
			 * in reads2 */
			ParseOptimize(chCard);
		}

		else if( strcmp(chKey4,"PHI(") == 0 )
		{
			/* enter phi(h), the number of h-ionizing photons/cm2 */
			i = 5;
			if( nqh >= LIMSPC )
			{
				/* too many continua were entered */
				fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* silly, but calms down the lint */
			ASSERT( nqh < LIMSPC );
			strcpy( rfield.chRSpec[nqh], "SQCM" );
			strcpy( rfield.chSpNorm[nqh], "PHI " );
			rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			/* set radius to very large value if not already set */
			/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
			if( !radius.lgRadiusKnown )
			{
				/* set radius to large value */
				ar1 = (float)radius.rdfalt;
				radius.Radius = pow(10.f,radius.rdfalt);
			}
			/* check if continuum so intense that crash is likely */
			if( rfield.totpow[nqh] > 29. )
			{
				fprintf( ioQQQ, " Is the flux for this continuum correct?\n" );
				fprintf( ioQQQ, " It appears too bright to me.\n" );
			}
			/* ParseRangeOption in readsun xx parse_rangeoption*/
			ParseRangeOption(nqh,chCard);
			/* vary option */
			if( optimize.lgVarOn )
			{
				/* >>chng 03 apr 30, had used log for range option, but if one or other
				 * was 1 ryd, then became 0 in the log, which was misinterpreted as
				 * one of the continuum bounds */
				strcpy( optimize.chVarFmt[optimize.nparm], "phi(h) %f range %f %f" );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				optimize.vparm[0][optimize.nparm] = (float)rfield.totpow[nqh];
				optimize.vincr[optimize.nparm] = 0.5;
				/* range option in, but cannot be varied */
				optimize.vparm[1][optimize.nparm] = (float)rfield.range[0][nqh];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[1][nqh];
				optimize.nvarxt[optimize.nparm] = 3;
				optimize.nparm += 1;
			}
			++nqh;
		}

		else if( strcmp(chKey4,"PLOT") == 0 )
		{
			/* make plot of log(nu.f(nu)) vs log(nu) or opacity
			 * in file plot */
			ParsePlot(chCard);
		}

		else if( strcmp(chKey4,"POWE") == 0 )
		{
			/* power law with cutoff, in reads2 */
			ParsePowerlawContinuum(chCard);
		}

		else if( strcmp(chKey4,"PRIN") == 0 )
		{
			/* adjust print schedule, in readsun */
			ParsePrint(chCard);
		}

		else if( strcmp(chKey4,"PUNC") == 0 )
		{
			/* punch something, in punch */
			ParsePunch(chCard);
		}

		else if( strcmp(chKey4,"Q(H)") == 0 )
		{
			/* log of number of ionizing photons */
			i = 5;
			if( nqh >= LIMSPC )
			{
				/* too many continua were entered */
				fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
				puts( "[Stop in ParseCommands]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* silly, but calms down the lint */
			ASSERT( nqh < LIMSPC );
			strcpy( rfield.chRSpec[nqh], "4 PI" );
			strcpy( rfield.chSpNorm[nqh], "Q(H)" );
			rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( rfield.totpow[nqh] > 100. && called.lgTalk )
			{
				fprintf( ioQQQ, " Is this reasonable?\n" );
			}
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			/* ParseRangeOption in readsun */
			ParseRangeOption(nqh,chCard);
			/* vary option */
			if( optimize.lgVarOn )
			{
				/* >>chng 03 apr 30, had used log for range option, but if one or other
				 * was 1 ryd, then became 0 in the log, which was misinterpreted as
				 * one of the continuum bounds */
				strcpy( optimize.chVarFmt[optimize.nparm], "Q(H) %f range %f %f" );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				optimize.vparm[0][optimize.nparm] = (float)rfield.totpow[nqh];
				optimize.vincr[optimize.nparm] = 0.5;
				/* range option in, but cannot be varied */
				optimize.vparm[1][optimize.nparm] = (float)rfield.range[0][nqh];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[1][nqh];
				optimize.nvarxt[optimize.nparm] = 3;
				++optimize.nparm;
			}
			/* increment number of luminosity sources specified */
			++nqh;
		}

		else if( strcmp(chKey4,"RATI") == 0 )
		{
			/* enter a continuum luminosity as a ratio of
			 * nuFnu for this continuum relative to a previous continuum
			 * format; first number is ratio of second to first continuum
			 * second number is energy for this ratio
			 * if third numbewr on line, then 2nd number is energy of
			 * first continuum, while 3rd number is energy of second continuum
			 * in reads2 */
			ParseRatio(chCard,&nqh);
		}

		else if( strcmp(chKey4,"RADI") == 0 )
		{
			/* log of inner and outer radii, default second=infinity,
			 * if R2<R1 then R2=R1+R2
			 * there is an optional keyword, "PARSEC" on the line
			 * to use PC as units, reads2 */
			ParseRadius(chCard,&ar1);
		}

		else if( strcmp(chKey4,"ROBE") == 0 )
		{
			/* this is the Roberto Terlivich command, no telling if it still works */
			radius.dRadSign = -1.;
		}

		else if( strcmp(chKey4,"SET ") == 0 )
		{
			/* set something, in reads2 */
			ParseSet(chCard);
		}

		else if( strcmp(chKey4,"SPEC") == 0 )
		{
			/* special key, can do anything */
			cdEXIT(EXIT_FAILURE);
		}

		else if( strcmp(chKey4,"SPHE") == 0 )
		{
			/* compute a spherical model, diffuse field from other side in
			 * in reads2 */
			ParseSphere(chCard);
		}

		else if( strcmp(chKey4,"STOP") == 0 )
		{
			/* stop model at desired zone, temperature, column density or tau-912
			 * in readsun */
			ParseStop(chCard);
		}

		else if( strcmp(chKey4,"TABL") == 0 )
		{
			/* interpolate on input tables for continuum, set of power  laws used
			 * input stored in big BLOCK data
			 * first check that this is the one and only INTERP command
			 * in readsun */
			ParseTable(&nqh,chCard,&ar1);
		}

		else if( strcmp(chKey4,"TAUM") == 0 )
		{
			/* minimum optical depths for lines (normally 1e-20)
			 *  (used in STARTER to set lines up) */
			i = 5;
			opac.taumin = (float)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
			if( lgEOL )
			{
				NoNumb(chCard);
			}
		}

		else if( strcmp(chKey4 , "TEST" ) == 0 )
		{
			char chStuff[INPUT_LINE_LENGTH];
			/* generate test input stream for rapid testing of code */

			/* hydrogen density */
			sprintf( chStuff , "HDEN 4 " );
			ParseHDEN(chStuff);

			/* make a constant temperature model */
			sprintf( chStuff , "CONSTANT TEMPER 4  " );
			ParseConstant(chStuff);

			/* continuum to include full energy range */
			sprintf( chStuff , "TABLE AGN  " );
			ParseTable(&nqh,chStuff,&ar1);

			/* set ionization parameter */
			sprintf( chStuff , "IONIZATION PARAMETER -2  " );
			ParseIonPar(&nqh,chStuff,&ar1);

			/* use old abundances */
			sprintf( chStuff , "ABUNDANCES OLD SOLAR 84  " );
			ParseAbundances(chStuff,lgDSet);

			/* >>chng 02 apr 19, add this */
			/* set total lyman continuum depth - this is to prevent caution
			 * that lyman continuum was thin but expected to be thick */
			sprintf( chStuff , "STOP LYMAN OPTICAL -4  " );
			ParseStop(chStuff);

			/* use largest possible hydrogen atom */
			if( lgMatch("LARG",chCard) )
			{
				sprintf( chStuff , "ATOM H-LIKE ELEMENT HYDROGEN LEVELS LIMIT  " );
				ParseAtomHLike(chStuff);
			}

			/* option to turn on the H2 molecule */
			if( lgMatch(" H2 ",chCard) )
			{
				/* this will both enable the molecule and also set the threshold for actually
				 * computing it to a very low value, so that it is actually used in ionized gas */
				sprintf( chStuff , "ATOM H2 LIMIT -20 " );
				ParseAtomH2(chStuff);
			}

			/* stop in second zone, so we do use the zone increment logic */
			sprintf( chStuff , "STOP ZONE 2  " );
			ParseStop(chStuff);

			/* set thickness */
			sprintf( chStuff , "SET DR 0  " );
			ParseSet(chStuff);

			/* this is option to also turn on large FeII atom */
			if( lgMatch("FEII",chCard) )
			{
				sprintf( chStuff , "ATOM FEII VERNER " );
				ParseAtomFeII(chStuff);
			}

			/* create series of assert commands */
			sprintf( input.chCARDCAPS , "ASSERT HYDROGEN 1 IONIZATION -3.040 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT HELIUM 2 IONIZATION -1.067 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT CARBON 2 IONIZATION -2.301 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT CARBON 3 IONIZATION -0.653 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT CARBON 4 IONIZATION -0.348 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT CARBON 5 IONIZATION -0.490 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT OXYGEN 3 IONIZATION  -0.800 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT OXYGEN 4 IONIZATION -0.180 " );
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT OXYGEN 5 IONIZATION -0.770 " );
			ParseAssertResults();

			/* >>chng 02 apr 19, from 0.7258 to 0.946, due to adding lyman cont depth */
			sprintf( input.chCARDCAPS , "ASSERT LINE \"CA B\" 4861 0.946 " );
			/* must have copy of this in chOrgCard, which is used by the routine to get lab */
			strcpy( input.chOrgCard , input.chCARDCAPS);
			ParseAssertResults();

			/* >>chng 02 apr 19, from 2.4603 to 3.25, due to adding lyman cont depth */
			sprintf( input.chCARDCAPS , "ASSERT LINE \"O  3\" 5007 3.25 " );
			/* must have copy of this in chOrgCard, which is used by the routine to get lab */
			strcpy( input.chOrgCard , input.chCARDCAPS);
			ParseAssertResults();

			sprintf( input.chCARDCAPS , "ASSERT HTOT -15.011" );
			ParseAssertResults();

		}

		else if( strcmp(chKey4,"TITL") == 0 )
		{
			/* read in title of model starting in col 5 */
			strcpy( input.chTitle , &input.chOrgCard[5] );
		}

		else if( strcmp(chKey4,"TLAW") == 0 )
		{
			/* some temperature vs depth law */
			ParseTLaw(chCard);
		}

		else if( strcmp(chKey4,"TOLE") == 0 )
		{
			fprintf(ioQQQ,
				"This command has been replaced with the SET TEMPERATURE TOLERANCE command.\n");
			fprintf(ioQQQ,
				"I will parse the command for now, but may not in a future version.\n");
			/* tolerance in heating cooling match, number is error/total */
			i = 5;
			conv.HeatCoolRelErrorAllowed = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			if( conv.HeatCoolRelErrorAllowed <= 0. )
			{
				conv.HeatCoolRelErrorAllowed = (float)pow(10.f,conv.HeatCoolRelErrorAllowed);
			}
		}

		else if( strcmp(chKey4,"TRAC") == 0 )
		{
			/* turn on trace output, in reads2 */
			ParseTrace(chCard);
		}

		else if( strcmp(chKey4,"TURB") == 0 )
		{
			/* turbulent velocity in km/sec */
			i = 5;
			DoppVel.TurbVel = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}

			if( lgMatch(" LOG",chCard) )
			{
				DoppVel.TurbVel = (float)pow(10.f,DoppVel.TurbVel);
			}
			/* convert to cm/sec */
			DoppVel.TurbVel *= 1e5;

			/* option to have turbulence be dissipative, keyword is dissipate,
			 * argument is log of scale length in cm */
			if( lgMatch("DISS",chCard) )
			{
				DoppVel.DispScale = (float)pow(10., FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) );
				if( lgEOL )
				{
					NoNumb(chCard);
				}
			}

			/* vary option */
			if( optimize.lgVarOn )
			{
				/* only one parameter to vary */
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "TURBULENCE= %f log" );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/* turbulent velocity */
				optimize.vparm[0][optimize.nparm] = (float)log10(DoppVel.TurbVel/
				  1e5);
				optimize.vincr[optimize.nparm] = 0.5;
				optimize.nparm += 1;
			}
		}

		else if( strcmp(chKey4,"WIND") == 0 )
		{
			/* NB - advection and wind commands are now a single command */
			/* parse the wind command, in dynamics.c */
			ParseDynaWind(chCard);
		}

		else if( strcmp(chKey2,"XI") == 0 )
		{
			/* inter x-ray ionization parameter xi=L/r^2 n;
			 * defined per hydrogen
			 * radius must also be entered if spherical, but not needed if plane */
			ParseIonPar(&nqh,chCard,&ar1);
		}

		else if( (chK1 != '#' && chK1 != '*') && chK1 != '%' )
		{
			/* end of keys - if we get here key is unrecognized---------- */
			fprintf( ioQQQ, "  Unrecognized command. Key=\"%4.4s\".  This is routine ParseCommands.\n", 
			  chKey4 );
			fprintf( ioQQQ, " The line image was==%s==\n", 
			  chCard );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseCommands]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* get next line image, this is tail of while loop that will 
		 * look for lgEOF or blank card */
		input_readarray(chCard,&lgEOF);
	}

	/*------------------------------------------------------------------- */
	/* fall through - hit lgEOF or blank line */

	/* >>chng 00 mar 31, removed dummy call to ParseQuantHeat, PvH */

	for( i=0; i < INPUT_LINE_LENGTH; i++ )
	{
		chCard[i] = ' ';
	}
	chCard[INPUT_LINE_LENGTH-1] = '\0';

	if( called.lgTalk )
	{
		fprintf( ioQQQ, "                       * %80.80s*\n", 
		  chCard );
		fprintf( ioQQQ, "                       ***********************************************************************************\n\n\n\n" );
	}

	/* lgStop_not_enough_info says that not enough info for model, so stop */
	lgStop_not_enough_info = FALSE;
	lgStop = FALSE;

	/* NSAVE is number of lines saved, =0 if no commands entered */
	if( input.nSave < 0 )
	{
		fprintf( ioQQQ, " No commands were entered - whats up?\n" );
		puts( "[Stop in ParseCommands]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
	{
		dense.gas_phase[ipHYDROGEN] = (float)log10(fabden(radius.Radius,radius.depth));
	}
	else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
	{
		/* >>chng 96 nov 29, added tabden, based on Kevin Volk code */
		dense.gas_phase[ipHYDROGEN] =(float)log10(tabden(radius.Radius,radius.depth));
	}

	/* perform sanity checks for some of the things read in, in reads2 
	readck(nqh);*/
	/* check hden within range. */
	if( dense.gas_phase[ipHYDROGEN] == -99. )
	{
		fprintf( ioQQQ, " Hydrogen density MUST be specified.\n" );
		lgStop_not_enough_info = TRUE;
		lgStop = TRUE;
	}

	/* fixup some variables */
	dense.gas_phase[ipHYDROGEN] = (float)pow(10.,dense.gas_phase[ipHYDROGEN]);

	radius.rinner = radius.Radius;

	/* mass loss rate for wind model, normally WINDV0=0 */
	wind.emdot = (float)(dense.gas_phase[ipHYDROGEN]*wind.windv0);

	/* iterate to convergence and wind models are mutually exclusive */
	if( wind.windv0 > 0. && conv.lgAutoIt )
	{
		fprintf( ioQQQ, " >>>Due to the nature of the Sobolev approximation, it makes no sense to converge a windy model.\n" );
		fprintf( ioQQQ, " >>>Iterate to convergence is turned off\n\n\n" );
		conv.lgAutoIt = FALSE;
		iterations.itermx = 0;
	}

	/* iterate to convergence and case b do not make sense together */
	/* WJH 22 May 2004: unless we are doing i-front dynamics (negative wind.windv0) */
	if( opac.lgCaseB && conv.lgAutoIt && wind.windv0 >= 0. )
	{
		fprintf( ioQQQ, " >>>Case B is an artifical test, it makes no sense to converge this model.\n" );
		fprintf( ioQQQ, " >>>Iterate to convergence is turned off\n\n\n" );
		conv.lgAutoIt = FALSE;
		iterations.itermx = 0;
	}

	/* specifying a density power and constant pressure makes no sense */
	if( dense.DensityPower!=0. && strcmp( dense.chDenseLaw, "CPRE" )==0 )
	{
		fprintf( ioQQQ, " Specifying both a density power law and constant pressure is impossible.\n" );
		lgStop = TRUE;
	}

	if( !rfield.lgIonizReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 )
	{
		fprintf( ioQQQ, " NO REEVALUATE IONIZATION can only be used with constant density.\n" );
		fprintf( ioQQQ, " Resetting to reevaluate ionization.\n\n" );
		rfield.lgIonizReevaluate = TRUE;
	}

	/* this is an option to turn on trace printout on the nth
	 * call from the optimizer - only allow trace if
	 * this is the case and nOptimiz 1 below nTrOpt */
	if( optimize.lgTrOpt )
	{
		/* nTrOpt was set with "optimize trace" command,
		 * is iteration to turn on trace */
		if( optimize.nTrOpt != optimize.nOptimiz )
		{
			trace.lgTrace = FALSE;
			/* following overrides turning on trace elsewhere */
			trace.lgTrOvrd = FALSE;
		}
		else
		{
			trace.lgTrace = TRUE;
			called.lgTalk = TRUE;
			trace.lgTrOvrd = TRUE;
			fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
		}
	}

	if( called.lgTalk && dense.gas_phase[ipHYDROGEN] < 1e-4 )
	{
		fprintf( ioQQQ, " >>>Is this value of the hydrogen density reasonable?\n" );
		fprintf( ioQQQ, " >>>It seems pretty low to me.\n\n\n" );
	}
	else if( called.lgTalk && dense.gas_phase[ipHYDROGEN] > 1e15 )
	{
		fprintf( ioQQQ, " >>>Is this value of the hydrogen density reasonable?\n" );
		fprintf( ioQQQ, " >>>It seems pretty high to me.\n\n\n" );
	}

	/* is the model going to crash because of extreme density? */
	if( called.lgTalk && !lgStop )
	{
		double ck;
		ck = 1e-37;
		ck /= 1e30;
		if( ck == 0. && ((dense.gas_phase[ipHYDROGEN] < -5.3) || (dense.gas_phase[ipHYDROGEN] > 19.2)) )
		{
			fprintf( ioQQQ, 
				" Simulation may crash because of extreme density.  "
				"Do not use a 32-bit machine, or make everything double precision.\n\n" );
		}
	}

	if( rfield.nspec == 0 )
	{
		fprintf( ioQQQ, " Type of continuum MUST be specified.\n" );
		lgStop = TRUE;
		lgStop_not_enough_info = TRUE;
	}

	if( nqh == 0 )
	{
		fprintf( ioQQQ, " Luminosity of continuum MUST be specified.\n" );
		lgStop = TRUE;
		lgStop_not_enough_info = TRUE;
	}

	/* set radius to very large value if not already set */
	/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
	/* NO! this is the one occurance where it must remain == 0 since flag not set 
	 * when some intensity commands sets r to very large value */
	if( radius.Radius == 0. )
	{
		fprintf( ioQQQ, " Starting radius MUST be specified.\n" );
		lgStop = TRUE;
		lgStop_not_enough_info = TRUE;
	}

	if( rfield.nspec != nqh )
	{
		fprintf( ioQQQ, " There were not the same number of continuum shapes and luminosities entered.\n" );
		lgStop = TRUE;
	}

	if( lgStop_not_enough_info )
	{
		fprintf( ioQQQ, " Insufficient information for model, STOP.\n" );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in ParseCommands]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( lgStop )
	{
		puts( "[Stop in ParseCommands]" );
		cdEXIT(EXIT_FAILURE);
	}

#	ifdef DEBUG_FUN
	fputs( " <->ParseCommands()\n", debug_fp );
#	endif
	return;
}
/*lint +e661 creation of  out of bound pointer */
