/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*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 "stopcalc.h"
#include "abund.h"
#include "geometry.h"
#include "dense.h"
#include "plot.h"
#include "rfield.h"
#include "dynamics.h"
#include "magnetic.h"
#include "trace.h"
#include "h2.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 "conv.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 */

	float a, 
	  ar1, 
	  teset, 
	  z;
	double thickness;

#	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 %s\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");
		}
		setbuf(ioQQQ,NULL);

		/* 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 )
		{
			/* accept both forms of feii */
			if( lgMatch("FEII",chCard) || lgMatch("FE II",chCard) )
			{
				/* parse the atom feii command - routine in atom_feii.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[rfield.nspec][0] = 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);
		}

		/* replace with CMB command, both will continue to be parsed */
		else if( (strcmp(chKey4,"CMB ") == 0) || (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 */
			ParseCMB(z,&nqh,&ar1);
		}

		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 density, log of rate relative to background, or density of rel. electrons,
			 * quantity is log unless keyword linear appears */
			ParseCosmicRays( chCard );
		}

		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 */
			ParseCrashDo(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 )
		{
			/* >>chng 05 dec 21, change diel command to set diel recom command */
			fprintf( ioQQQ, " The DIELectronic command has been replaced with the SET DIELectronic recombination command.\n" );
			fprintf( ioQQQ, " Please have a look at Hazy.\n Sorry.\n\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;
			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" );
			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.,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], "LUMI" );

			/* full range of continuum will be used */
			rfield.range[nqh][0] = rfield.emm;
			rfield.range[nqh][1] = rfield.egamry;
			rfield.totpow[nqh] = (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]);
				optimize.vincr[optimize.nparm] = 0.1f;
				optimize.nvarxt[optimize.nparm] = 1;
				optimize.nparm += 1;
			}

			/* finally increment number of continua */
			++nqh;
		}

		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 at 1 Ryd rather than 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, " This command should not be necessary.\n" );
				fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this simulation.\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,"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.,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[nqh][0]);
				optimize.vparm[2][optimize.nparm] = (float)log10(rfield.range[nqh][1]);*/
				/* >>change to linear option */
				optimize.vparm[1][optimize.nparm] = (float)rfield.range[nqh][0];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[nqh][1];
				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[rfield.nspec][0] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				/* default width is 0.05 */
				rfield.cutoff[rfield.nspec][0] = 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);
			}

			strcpy( rfield.chRSpec[nqh], "4 PI" );
			rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				NoNumb(chCard);
			}
			if( lgMatch("LINE",chCard) )
			{
				/* option for linear input parameter */
				rfield.totpow[nqh] = log10(rfield.totpow[nqh]);
			}

			strcpy( rfield.chSpNorm[nqh], "LUMI" );

			/* the solar option - number is log total solar luminosity */
			if( lgMatch("SOLA",chCard) )
			{
				/* option to use log of total lumin in solar units */
				rfield.range[nqh][0] = rfield.emm;
				rfield.range[nqh][1] = rfield.egamry;
				/* log of solar luminosity in erg s-1 */
				rfield.totpow[nqh] += 33.5827f;
			}
			else
			{
				/* check if range is present and parse it if it is - ParseRangeOption in readsun 
				 * this includes TOTAL keyword for total luminosity */
				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[nqh][0];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[nqh][1];
				optimize.nvarxt[optimize.nparm] = 3;
				optimize.nparm += 1;
			}
			++nqh;
		}

		else if( strcmp(chKey4,"MAGN") == 0 )
		{
			/* parse the magnetic field command, routine in magnetic.c */
			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.,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[nqh][0];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[nqh][1];
				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[nqh][0];
				optimize.vparm[2][optimize.nparm] = (float)rfield.range[nqh][1];
				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) || lgMatch("FE II",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 )
		{
			/* line has turbulent velocity in km/s */
			i = 5;
			DoppVel.TurbVel = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* true if no number on line - check for equipartition */
			if( lgEOL )
			{
				if( lgMatch("EQUI",chCard) )
				{
					/* turbulence equipartition option */
					DoppVel.lgTurbEquiMag = TRUE;
				}
				else
					NoNumb(chCard);

			}

			if( lgMatch(" LOG",chCard) )
			{
				DoppVel.TurbVel = (float)pow(10.f,DoppVel.TurbVel);
			}
			/* convert from km/s to cm/s */
			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.f,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 )
	{
		if( called.lgTalk )
		{
			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. )
	{
		if( called.lgTalk )
		{
			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 )
	{
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " Specifying both a density power law and constant pressure is impossible.\n" );
		}
		lgStop = TRUE;
	}

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

	/* check if the combination of stopping column density and H density are physically plausible */
	thickness = vfmin( StopCalc.tauend,
		StopCalc.colpls,
		StopCalc.colnut,
		StopCalc.col_h2,
		StopCalc.col_h2_nut , 
		StopCalc.HColStop , 
		FEND );
	if( thickness < COLUMN_INIT )
	{
		/* a stop column density was specified - check on physical thickness this corresponds to */
		thickness /= (dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
		if( thickness > 1e25 )
		{
			fprintf( ioQQQ, 
				">>> The specified column density and hydrogen density correspond to a thickness of %.2e cm.\n",
				thickness);
			fprintf( ioQQQ, 
				">>> This seems large to me.\n");
			if( lgBit32 )
				fprintf(ioQQQ,">>> This is a 32-bit CPU- a very large radius may cause overflow.\n\n");
		}
	}

	/* >>chng 04 nov 27, test on these two */
	if( !rfield.lgDoLineTrans && rfield.lgOpacityFine )
	{
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " NO LINE TRANSER set but fine opacities still computed.\n" );
			fprintf( ioQQQ, " Turning off fine opacities.\n\n" );
		}
		rfield.lgOpacityFine = FALSE;
	}

	/* >>chng 04 nov 27, test of these two */
	if( h2.lgH2ON && (!rfield.lgDoLineTrans || !rfield.lgOpacityFine) )
	{
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " Large H2 molecule turned on but line transfer and fine opacities are not.\n" );
			fprintf( ioQQQ, " Turning on line transfer and fine opacities.\n\n" );
		}
		rfield.lgOpacityFine = TRUE;
		rfield.lgDoLineTrans = TRUE;
	}

	/* print warning if temperature set below default but cosmic rays not turned on */
	if( (hextra.cryden == 0.) && called.lgTalk  && 
		(StopCalc.tend < TEND || 
		/* thermal.ConstTemp def is zero, set pos when used */
		(thermal.ConstTemp > 0. && thermal.ConstTemp < TEND ) ) )
	{
		fprintf( ioQQQ, "\n"
			            " >>> The simulation is going into neutral gas but cosmic rays are not included.\n" );
		fprintf( ioQQQ, " >>> Ion-molecule chemistry will not occur without a source of ionization.\n" );
		fprintf( ioQQQ, " >>> The chemistry network may collapse deep in molecular regions.\n" );
		fprintf( ioQQQ, " >>> Consider adding galactic background cosmic rays with the COSMIC RAYS BACKBOUND command.\n" );
		fprintf( ioQQQ, " >>> You may be making a big mistake.\n\n\n\n\n" );
	}

	/* print warning if temperature set below default but C > O */
	if( called.lgTalk  &&  dense.gas_phase[ipCARBON]/SDIV( dense.gas_phase[ipOXYGEN]) >= 1. )
	{
		fprintf( ioQQQ, "\n"
			            " >>> The simulation is going into possibly molecular gas but the carbon/oxygen abundance ratio is greater than unity.\n" );
		fprintf( ioQQQ, " >>> The chemistry network is designed for environments with C/O < 1.\n" );
		fprintf( ioQQQ, " >>> The chemistry network may (or may not) collapse deep in molecular regions where CO is fully formed.\n" );
		fprintf( ioQQQ, " >>> This is a limitation of the current chemistry network.\n\n\n\n\n" );
	}

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

	/* set coverage - limit number of iterations and zones */
	if( iterations.lgConverge_set)
	{
		iterations.itermx = MIN2( iterations.itermx , iterations.lim_iter );
		for( j=0; j < ITRDIM; j++ )
		{
			geometry.nend[j] = MIN2( geometry.nend[j] , iterations.lim_zone );
		}
	}

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