/* 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 */
/*ParsePunch parse the punch command */
/*PunchFilesInit initialize punch file pointers, called from cdInit */
/*ClosePunchFiles close all punch files */
/*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present */
#include "cddefines.h"
#include "physconst.h"
#include "elementnames.h"
#include "input.h"
#include "geometry.h"
#include "opacity.h"
#include "rfield.h"
#include "map.h"
#include "atomfeii.h"
#include "h2.h"
#include "co.h"
#include "hmi.h"
#include "version.h"
#include "grainvar.h"
#include "punch.h"
#include "parse.h"
/* check for keyword UNITS on line, then scan wavelength or energy units if present */
static void ChkUnits( char *chCard);

/* option to append instead of overwrite */
static int lgNoClobber[LIMPUN];
/* these are for some special cases, same purpose as prev no clobber */
static int lgPunConv_noclobber , lgDROn_noclobber , lgPunH2_noclobber , 
	lgPunPoint_noclobber , lgioRecom_noclobber , QHPunchFile_noclobber;

/* NB NB NB NB NB NB NB NB NB NB
 *
 * if any "special" punch commands are added (commands that copy the file pointer
 * into a separate variable, e.g. like PUNCH _DR_), be sure to add that file pointer
 * to PunchFilesInit and ClosePunchFiles !!!
 *
 * PUNCH FILE POINTERS SHOULD NEVER BE ALTERED BY ROUTINES OUTSIDE THIS MODULE !!!
 *
 * hence initializations of punch file pointers should never be included in zero() !!
 * the pointer might be lost without the file being closed
 * 
 * NB NB NB NB NB NB NB NB NB NB */


void ParsePunch(char *chCard )
{
	char chLabel[INPUT_LINE_LENGTH];
	char chFilename[INPUT_LINE_LENGTH];
	int lgEOL ;
	/* pointer to column across line image for free format number reader*/
	long int ipFFmt, 
	  i,
	  j,
	  nelem;

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

	/* while optimizing, punch commands are not executed -> return */
	if( !punch.lgOpenUnits )
	{
#		ifdef DEBUG_FUN
		fputs( " <->ParsePunch()\n", debug_fp );
#		endif
		return;
	}

	/* check that limit not exceeded */
	if( punch.npunch > LIMPUN )
	{
		fprintf( ioQQQ, 
			"The limit to the number of PUNCH options is %ld.  Increase LIMPUN in punch.h if more are needed.\nSorry.\n", 
		  LIMPUN );
		puts( "[Stop in ParsePunch]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* LAST keyword is an option to punch only on last iteration */
	punch.lgPunLstIter[punch.npunch] = lgMatch("LAST",chCard);

	/* 
	 * get file name for this punch output.
	 * GetQuote does the following -
	 * first copy original version of file name into chLabel, 
	 * string does include null termination.
	 * set filename in OrgCard and second parameter to spaces so 
	 * that not picked up below as keyword
	 * last parameter says to abort if no quote found 
	 */
	GetQuote( chLabel , chCard , TRUE );

	/* check that name is not same as opacity.opc, a special file */
	if( strcmp(chLabel , "opacity.opc") == 0 )
	{
		fprintf( ioQQQ, "ParsePunch will not allow punch file name %s, please choose another.\nSorry.\n",
			chFilename);
		puts( "[Stop in ParsePunch]" );
		cdEXIT(EXIT_FAILURE);
	}
	else if( chLabel[0]==0 )
	{
		fprintf( ioQQQ, "ParsePunch found a null file name between double quotes, please check command line.\nSorry.\n");
		puts( "[Stop in ParsePunch]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now copy to chFilename, with optional prefix first */
	/* this is optional prefix, normally a null string, set with set punch prefix command */
	strcpy( chFilename , punch.chFilenamePrefix );
	strcat( chFilename , chLabel );

	/* only open if file has not already been opened during a previous call */
	if( punch.ipPnunit[punch.npunch] == NULL )
	{
		/* punch.npunch is set to 0 in ZERO, 
		 * punch.ipPnunit is the file pointer  */

		/* open the file with the name generated above */
		if( (punch.ipPnunit[punch.npunch] = fopen( chFilename,"w" ) ) == NULL )
		{
			fprintf( ioQQQ, "ParsePunch could not open punch file %s for writing.\nSorry.\n",
				chFilename);
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* NO CLOBBER no clobber keyword is an option to append to a given file 
	 * instead of overwriting it */
	if( lgMatch("NO C",chCard) && lgMatch("CLOB",chCard) )
		lgNoClobber[punch.npunch] = TRUE;

	/* option to set no buffering for this file */
	if( lgMatch("NO BUFFER",chCard) )
		setbuf( punch.ipPnunit[punch.npunch] , NULL );

	/* put version number and title of model on output file, but only if
	 * this is requested with a "title" on the line*/
	/* >>chng 02 may 10, invert logic from before - default had been title */
	/* put version number and title of model on output file, but only if
	 * there is not a "no title" on the line*/
	if( !lgMatch("O TI",chCard) && lgMatch("TITL",chCard))
	{
		fprintf( punch.ipPnunit[punch.npunch], 
			"# %s %s\n", 
		  version.chVersion, input.chTitle );
	}

	/* usually results for each iteration are followed by a series of
	 * hash marks, ####, which fool excel.  This is an option to not print
	 * the line.  If the keyword NO HASH appears anywhere, on any punch command,
	 * the hash marks will never be produced */
	if( lgMatch("O HA",chCard) )
	{
		punch.lgHashEndIter = FALSE;
	}

	/* punch opacity must come first since elements appear as sub-keywords */
	if( lgMatch("OPAC",chCard) )
	{

		/* check for keyword UNITS on line, then scan wavelength or energy units if present,
		 * units are copied into punch.chConPunEnr */
		ChkUnits(chCard);

		strcpy( punch.chPunch[punch.npunch], "OPAC" );
		if( lgMatch("TOTA",chCard) )
		{
			/* DoPunch will call popac to parse the subcommands
			 * punch total opacity */
			strcpy( opac.chOpcTyp, "TOTL" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#nu\tTot opac\tAbs opac\tScat opac\tAlbedo\telem\n" );
		}

		else if( lgMatch("FIGU",chCard) )
		{
			/* do figure for hazy */
			strcpy( opac.chOpcTyp, "FIGU" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#nu, H, He, tot opac\n" );
		}

		else if( lgMatch("FINE",chCard) )
		{
			/* punch out the fine opacity array with only lines */
			rfield.lgPunchOpacityFine = TRUE;
			strcpy( opac.chOpcTyp, "FINE" );
			/* check for keyword UNITS on line, then scan wavelength or energy units if present,
			* units are copied into punch.chConPunEnr */
			ChkUnits(chCard);

			fprintf( punch.ipPnunit[punch.npunch], 
				"#nu\topac\n" );
			ipFFmt = 5;
			/* range option - important since so much data */
			if( lgMatch("RANGE",chCard) ) 
			{
				/* get lower and upper range, must be in Ryd */
				punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
				punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
				if( lgEOL )
				{
					fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
					puts( "[Stop in ParsePunch]" );
					cdEXIT(EXIT_FAILURE);
				}
				if( punch.punarg[0][punch.npunch] >=punch.punarg[1][punch.npunch] )
				{
					fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
					puts( "[Stop in ParsePunch]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			else
			{
				/* these mean full energy range */
				punch.punarg[0][punch.npunch] = 0.;
				punch.punarg[1][punch.npunch] = 0.;
			}
			/* optional last parameter - how many points to bring together */
			punch.punarg[2][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);

			/* default is to print every tenth one */
			if( lgEOL )
				punch.punarg[2][punch.npunch] = 10;
			if( punch.punarg[2][punch.npunch] < 1 )
			{
				fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( lgMatch("GRAI",chCard) )
		{
			/* punch grain opacity command, give optical properties of gains in calculation */
			strcpy( punch.chPunch[punch.npunch], "DUSO" );
			/* punch grain opacity command in twice, here and above in opacity */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#grain nu\ttotal(1-g)\tabst\tscat(1-g)\tscat\n" );
		}

		else if( lgMatch("BREM",chCard) )
		{
			/* punch brems opacity */
			strcpy( opac.chOpcTyp, "BREM" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#nu\tbrem opac\n" );
		}

		else if( lgMatch("SHEL",chCard) )
		{
			/* punch shells, a form of the punch opacity command for showing subshell crossections*/
			strcpy( punch.chPunch[punch.npunch], "OPAC" );

			/* punch subshell cross sections */
			strcpy( opac.chOpcTyp, "SHEL" );

			/* this is element */
			ipFFmt = 3;
			punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);

			/* this is ion */
			punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);

			/* this is shell */
			punch.punarg[2][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);

			if( lgEOL )
			{
				fprintf( ioQQQ, "There must be io unit, atom wght, ion, shell\nSorry.\n" );
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}
			fprintf( punch.ipPnunit[punch.npunch], 
				"#sub shell cross section\n" );
		}

		else
		{
			/* punch elenent opacity, produces n name.n files, one for each stage of 
			 * ionization.  the name is the 4-char version of the element's name, and
			 * n is the stage of ionization.  the file nameon the card is ignored.
			 * The code Stop in getpunchs after these files are produced. */

			/* this will be used as check that we did find an element on the command lines */
			/* nelem is -1 if an element was not found */
			if( (nelem = GetElem( chCard ) ) < 0 )
			{
				fprintf( ioQQQ, "There must be a second keyword on the opacity command.  Sorry.\n" );
				fprintf( ioQQQ, "Options are total, figure, grain, or element name.\n" );
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* copy string over */
			strcpy( opac.chOpcTyp, elementnames.chElementNameShort[nelem] );
		}
	}

	/* punch H2 has to come early since it has many suboptions */
	else if( lgMatch(" H2 ",chCard) )
	{
		/* this is in h2_io.c */
		H2_ParsePunch( chCard );
	}

	else if( lgMatch("ABUN",chCard) )
	{
		/* punch abundances */
		strcpy( punch.chPunch[punch.npunch], "ABUN" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#abundances\n" );
	}

	else if( lgMatch("AGES",chCard) )
	{
		/* punch ages */
		strcpy( punch.chPunch[punch.npunch], "AGES" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#ages\n" );
	}

	else if( lgMatch(" AGN",chCard) )
	{
		/* punch tables needed for AGN3 */
		strcpy( punch.chPunch[punch.npunch], " AGN" );
		/* this is the AGN option, to produce a table for AGN */

		/* charge exchange rate coefficients */
		if( lgMatch("CHAR",chCard) )
		{
			strcpy( punch.chPunch[punch.npunch], "CHAG" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#charge exchange rate coefficnt\n" );
		}

		else if( lgMatch("RECO",chCard) )
		{
			/* punch recombination rates for agn table */
			strcpy( punch.chPunch[punch.npunch], "RECA" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Recomrates for agn table\n" );
		}

		else if( lgMatch("OPAC",chCard) )
		{
			/* create table for appendix in AGN */
			strcpy( opac.chOpcTyp, " AGN" );
			strcpy( punch.chPunch[punch.npunch], "OPAC" );
		}

		else if( lgMatch("HECS",chCard) )
		{
			/* create table for appendix in AGN */
			strcpy( punch.chPunchArgs[punch.npunch], "HECS" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#agn he cs \n" );
		}

		else if( lgMatch("HEMI",chCard) )
		{
			/* hemis - continuum emission needed for chap 4 */
			strcpy( punch.chPunchArgs[punch.npunch], "HEMI" );
			/* routine will make its own title */
			/*fprintf( punch.ipPnunit[punch.npunch], "#agn wl contin \n" );*/

			/* check for keyword UNITS on line, then scan wavelength or energy units if present,
			 * units are copied into punch.chConPunEnr */
			ChkUnits(chCard);
		}
		else if( lgMatch("RECC",chCard) )
		{
			/* recombination cooling, for AGN */
			strcpy( punch.chPunch[punch.npunch], "HYDr" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#T\tbAS\tb1\tbB\n" );
		}
		else
		{
			fprintf( ioQQQ, " I did not recognize this .\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}

	}

	else if( lgMatch("ASSE",chCard) )
	{
		/* punch asserts */
		strcpy( punch.chPunch[punch.npunch], "ASSE" );
		/* no need to print this standard line of explanation*/
		/*fprintf( punch.ipPnunit[punch.npunch], " asserts\n" );*/
	}

	/* punch charge transfer */
	else if( lgMatch("CHAR",chCard) && lgMatch("TRAN",chCard) )
	{
		/* NB in dopunch only the first three characters are compared to find this option,
		 * search for "CHA" */
		/* punch charge transfer */
		strcpy( punch.chPunch[punch.npunch], "CHAR" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#charge exchange rate coefficnt\n" );
	}

	/* in column density need both since column is also keyword for one line per line */
	else if( lgMatch("COLU",chCard) && lgMatch("DENS",chCard) )
	{
		/* punch column density */
		strcpy( punch.chPunch[punch.npunch], "COLU" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#column densities\n" );
	}

	else if( lgMatch("COMP",chCard) )
	{
		/* punch compton, details of the energy exchange problem */
		strcpy( punch.chPunch[punch.npunch], "COMP" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#nu, comup, comdn\n" );
	}

	else if( lgMatch("COOL",chCard) )
	{
		/* cooling, actually done by routine */
		strcpy( punch.chPunch[punch.npunch], "COOL" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#temp\theat\tcool\tfrac coolant \n" );
	}

	else if( lgMatch("DYNA",chCard) )
	{

		/* punch something dealing with dynamics 
		 * in Punchdo the DYN part of key is used to call DunaPunch,
		 * with the 4th char as its second argument.  DynaPunch uses that
		 * 4th letter to decide the job */
		if( lgMatch( "ADVE" , chCard ) )
		{
			/* punch information relating to advection */
			strcpy( punch.chPunch[punch.npunch], "DYNa");
			fprintf( punch.ipPnunit[punch.npunch], 
				"#avection depth\thtot\tCoolHeat\tdCoolHeatdT\tRecomb\tPhoto\ttimestp\n" );
		}
		else
		{
			fprintf( ioQQQ, " I did not recognize a subkeyword on this PUNCH DYNAMICS command.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	else if( lgMatch("ENTH",chCard) )
	{
		/* contributors to the total enthalpy */
		strcpy( punch.chPunch[punch.npunch], "ENTH" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tTotal\tExcit\tIoniz\tBind\tKE\tther+PdV\tmag \n" );
	}

	else if( lgMatch("FEII",chCard) )
	{
		/* something to do with FeII atom - options are lines and continuum */
		/* FeII continuum, only valid if large atom is set */
		if( lgMatch("CONT",chCard) )
		{
			/* punch feii continuum */
			strcpy( punch.chPunch[punch.npunch], "FEco" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#FeII: wavelength(A)\tsummed intensity\n" );
		}
		else if( lgMatch("LINE",chCard) )
		{
			/* punch feii lines command , three optional parameters, threshold for faintest
			* line to print, lower and upper energy bounds */

			/* punch intensities from large FeII atom */
			strcpy( punch.chPunch[punch.npunch], "FEli" );

			/* short keyword makes punch half as big */
			if( lgMatch("SHOR",chCard) )
			{
				FeII.lgShortFe2 = TRUE;
			}
			else
			{
				FeII.lgShortFe2 = FALSE;
			}

			/* first optional number changes the theshold of weakest line to print*/
			i = 5;
			/* fe2thresh is intensity relative to normalization line,
			* normally Hbeta, and is set to zero in zero.c */
			FeII.fe2thresh = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				FeII.fe2thresh = 0.;
			}

			/* it is a log if negative */
			if( FeII.fe2thresh < 0. )
			{
				FeII.fe2thresh = (float)pow(10.f,FeII.fe2thresh);
			}

			/* check for energy range of lines to be punched,
			* this is to limit size of output file */
			FeII.fe2ener[0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				FeII.fe2ener[0] = 0.;
			}

			FeII.fe2ener[1] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				FeII.fe2ener[1] = 1e8;
			}
			/* if either is negative then both are logs */
			if( FeII.fe2ener[0] < 0. || FeII.fe2ener[1] < 0. )
			{
				FeII.fe2ener[0] = (float)pow(10.f,FeII.fe2ener[0]);
				FeII.fe2ener[1] = (float)pow(10.f,FeII.fe2ener[1]);
			}
			/* entered in Ryd in above, convert to wavenumbers */
			FeII.fe2ener[0] /= (float)WAVNRYD;
			FeII.fe2ener[1] /= (float)WAVNRYD;

			/* these results are actually created by the FeIIPunchLines routine
			* that lives in the FeIILevelPops file */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Results from Verner FeII atom.\n" );
		}

		else if( lgMatch("FRED",chCard) )
		{
			/* punch out some stuff for Fred's FeII atom */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#FeII rates \n" );
			strcpy( punch.chPunch[punch.npunch], "FE2f" );
		}
		else if( lgMatch("DEPA",chCard) )
		{
			/* punch out departure coef for the large FeII atom */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#FeII departure coef \n" );
			/* optional keyword all means do all levels, if not then just do subset */
			if( lgMatch(" ALL",chCard) )
			{
				/* punch all levels, calls routine FeIIPunDepart */
				strcpy( punch.chPunch[punch.npunch], "FE2D" );
			}
			else
			{
				/* punch a very few selected levels, calls routine FeIIPunDepart */
				strcpy( punch.chPunch[punch.npunch], "FE2d" );
			}
		}
		else if( lgMatch("POPU",chCard) )
		{
			/* punch out level populations for the large FeII atom */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#FeII level populations \n" );
			/* optional keyword all means do all levels, if not then just do subset */
			if( lgMatch(" ALL",chCard) )
			{
				/* punch all levels, calls routine FeIIPunDepart */
				strcpy( punch.chPunch[punch.npunch], "FE2P" );
			}
			else
			{
				/* punch a very few selected levels, calls routine FeIIPunDepart */
				strcpy( punch.chPunch[punch.npunch], "FE2p" );
			}
		}
		else
		{
			fprintf( ioQQQ, "There must be a second keyword on this command.\n" );
			fprintf( ioQQQ, "The ones I know about are departure and fred.\n" );
			fprintf( ioQQQ, "Sorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* the punch continuum command, with many options,
	 * the first 3 char of the chPunch flag will always be "CON" 
	 * with the last indicating which one */
	else if( lgMatch("CONT",chCard) )
	{
		/* this flag is checked in PrtComment to generate a caution
		 * if continuum is punched but iterations not performed */
		punch.lgPunContinuum = TRUE;

		/* check for keyword UNITS on line, then scan wavelength or energy units if present,
		 * units are copied into punch.chConPunEnr */
		ChkUnits(chCard);

		if( lgMatch("BINS",chCard) )
		{
			/* continuum binning */
			strcpy( punch.chPunch[punch.npunch], "CONB" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Continuum binning, enrOrg, Enr, width of cells\n" );
		}

		else if( lgMatch("DIFF",chCard) )
		{
			/* diffuse continuum, the locally emitted continuum stored in rfield.ConEmitLocal */
			strcpy( punch.chPunch[punch.npunch], "COND" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#energy\t ConEmitLocal \n" );
		}

		else if( lgMatch("EMIT",chCard) )
		{
			/* continuum emitted by cloud */
			strcpy( punch.chPunch[punch.npunch], "CONE" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Energy\treflec\toutward\ttotal\tline\tcont\n" );
		}

		else if( lgMatch("FINE" , chCard ) )
		{
			rfield.lgPunchOpacityFine = TRUE;
			/* fine transmitted continuum cloud */
			strcpy( punch.chPunch[punch.npunch], "CONf" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Energy\tTransmitted\n" );

			/* range option - important since so much data */
			if( lgMatch("RANGE",chCard) ) 
			{
				/* get lower and upper range, must be in Ryd */
				ipFFmt = 5;
				punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
				punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
				if( lgEOL )
				{
					fprintf(ioQQQ,"There must be two numbers, the lower and upper energies in Ryd.\nSorry.\n");
					puts( "[Stop in ParsePunch]" );
					cdEXIT(EXIT_FAILURE);
				}
				if( punch.punarg[0][punch.npunch] >=punch.punarg[1][punch.npunch] )
				{
					fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
					puts( "[Stop in ParsePunch]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			else
			{
				/* these mean full energy range */
				punch.punarg[0][punch.npunch] = 0.;
				punch.punarg[1][punch.npunch] = 0.;
			}
			/* optional last parameter - how many points to bring together */
			punch.punarg[2][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);

			/* default is to bring together ten */
			if( lgEOL )
				punch.punarg[2][punch.npunch] = 10;

			if( punch.punarg[2][punch.npunch] < 1 )
			{
				fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else if( lgMatch("GRAI",chCard) )
		{
			/* grain continuum */
			strcpy( punch.chPunch[punch.npunch], "CONG" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#energy,     graphite,  rest \n" );
		}

		else if( lgMatch("INCI",chCard) )
		{
			/* incident continuum */
			strcpy( punch.chPunch[punch.npunch], "CONC" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Incident Continuum, Enr\tnFn \n" );
		}

		else if( lgMatch("INTE",chCard) )
		{
			/* continuum interactions */
			strcpy( punch.chPunch[punch.npunch], "CONi" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Continuum interactions, inc, otslin. otscon, ConInterOut, outlin \n" );
			/* this is option for lowest energy, if nothing then zero */
			ipFFmt = 3;
			punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);
		}

		else if( lgMatch("IONI",chCard) )
		{
			/* punch ionizing continuum*/
			strcpy( punch.chPunch[punch.npunch], "CONI" );

			/* this is option for lowest energy, if nothing then zero */
			ipFFmt = 3;
			punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);

			/* this is option for smallest interaction to punch, def is 1 percent */
			punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				punch.punarg[1][punch.npunch] = 0.01f;
			}

			/* >>chng 03 may 03, add "every" option to punch this on every zone -
			 * if every is not present then only last zone is punched */
			if( lgMatch("EVER", chCard ) )
			{
				/* punch every zone */
				punch.punarg[2][punch.npunch] = 1;
			}
			else
			{
				/* only punch last zone */
				punch.punarg[2][punch.npunch] = 0;
			}

			/* put the header at the top of the file */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#cell\tnu\tflux\tflx*cs\tFinc\totsl\totsc\toutlin\toutcon\trate/tot\tintegral\tline\tcont\n" );
		}

		else if( lgMatch("OUTW",chCard) )
		{
			/* outward only continuum */
			if( lgMatch("LOCA",chCard) )
			{
				strcpy( punch.chPunch[punch.npunch], "CONo" );
				fprintf( punch.ipPnunit[punch.npunch], 
					"#Local Out   ConInterOut+line SvOt*opc pass*opc\n" );
			}
			else
			{
				strcpy( punch.chPunch[punch.npunch], "CONO" );
				fprintf( punch.ipPnunit[punch.npunch], 
					"#Out Con      OutIncid  OutConD   OutLinD   OutConS\n" );
			}
		}

		else if( lgMatch("TRAN",chCard) )
		{
			/* transmitted continuum */
			strcpy( punch.chPunch[punch.npunch], "CONT" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Transmitted Continuum\n" );
		}

		else if( lgMatch(" TWO",chCard) )
		{
			/* total two photon continua rfield.TotDiff2Pht */
			strcpy( punch.chPunch[punch.npunch], "CON2" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#energy\t n_nu\tnuF_nu \n" );
		}

		else if( lgMatch(" RAW",chCard) )
		{
			/* "raw" continua */
			strcpy( punch.chPunch[punch.npunch], "CORA" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Raw Con anu\tflux\totslin\totscon\tConRefIncid\tConEmitReflec\tConInterOut\toutlin\tConEmitOut\tline\tcont\tnLines\n" );
		}

		else if( lgMatch("REFL",chCard) )
		{
			/* reflected continuum */
			strcpy( punch.chPunch[punch.npunch], "CONR" );

			fprintf( punch.ipPnunit[punch.npunch], 
				"#Reflected      cont       line       total        albedo    ConID\n" );
		}

		else
		{
			/* this is the usual punch continuum command */
			strcpy( punch.chPunch[punch.npunch], "CON " );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Cont nu\tincident\ttrans\tDiffOut\tnet trans\treflc\ttotal\tline\tcont\tnLine\n" );
		}

	}

	/* punch information about convergence of this model 
	 * reason - why it did not converge an iteration
	 * error - zone by zone display of various convergence errors */
	else if( lgMatch("CONV",chCard) )
	{
		if( lgMatch("REAS",chCard) )
		{
			/* punch reason model declared not converged
			 * not a true punch command, since done elsewhere */
			punch.ipPunConv = punch.ipPnunit[punch.npunch];
			punch.lgPunConv = TRUE;
			fprintf( punch.ipPunConv, 
				"# reason for continued iterations\n" );
			/* this does not count as a punch option (really) */
			punch.ipPnunit[punch.npunch] = NULL;
			--punch.npunch;
			/* all punch files have the "no clobber" option */
			if( lgMatch("NO C",chCard) && lgMatch("CLOB",chCard) )
				lgPunConv_noclobber = TRUE;
		}
		else if( lgMatch("ERRO",chCard) )
		{
			/* punch zone by zone errors in pressure, electron density, and heating-cooling*/
			punch.ipPunConv = punch.ipPnunit[punch.npunch];
			punch.lgPunConv = TRUE;
			/* convergence error */
			strcpy( punch.chPunch[punch.npunch], "CNVE" );
			fprintf( punch.ipPunConv, 
				"#depth\tnPres2Ioniz\tP(cor)\tP(cur)\tP%%error\tNE(cor)\tNE(cur)\tNE%%error\tHeat\tCool\tHC%%error\n" );
		}
		else
		{
			fprintf( ioQQQ, "There must be a second keyword on this command.\n" );
			fprintf( ioQQQ, "The ones I know about are reason and error.\n" );
			fprintf( ioQQQ, "Sorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	else if( lgMatch(" DR ",chCard) )
	{
		punch.ipDRout = punch.ipPnunit[punch.npunch];
		punch.lgDROn = TRUE;

		/* set punch last flag to whatever it was above */
		punch.lgDRPLst = punch.lgPunLstIter[punch.npunch];

		fprintf( punch.ipDRout, 
			"#zone\tdepth\tdr\tdr 2 go\treason \n" );

		/* this does not count as a punch option (really) */
		punch.ipPnunit[punch.npunch] = NULL;
		punch.npunch -= 1;
		/* all punch files have the "no clobber" option */
		if( lgMatch("NO C",chCard) && lgMatch("CLOB",chCard) )
			lgDROn_noclobber = TRUE;
	}

	else if( lgMatch("ELEM",chCard) && !lgMatch("GAMMA" , chCard) )
	{
		/* option to punch ionization structure of some element */
		strcpy( punch.chPunch[punch.npunch], "ELEM" );
		/* this returns element number on c scale */
		j = GetElem(chCard);
		++j;
		/* this is the atomic number on the physical scale */
		punch.punarg[0][punch.npunch] = (float)j;
		ASSERT( j>0);
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth");
		for(i=0; i<j+1;++i )
		{
			ASSERT( j>0 );
			fprintf( punch.ipPnunit[punch.npunch], 
				"\t%2s%2li", elementnames.chElementSym[j-1],i+1);
		}
		fprintf( punch.ipPnunit[punch.npunch], "\n");
	}

	else if( lgMatch("GAMM",chCard) )
	{
		/* punch all photoionizaiton rates for all subshells */
		fprintf( punch.ipPnunit[punch.npunch], 
			"#Photoionization rates \n" );
		if( lgMatch("ELEMENT" , chCard ) )
		{
			/* element keyword, find element name and stage of ionization, 
			 * will print photoionization rates for valence of that element */
			strcpy( punch.chPunch[punch.npunch], "GAMe" );

			/* this returns element number on c scale */
			j = GetElem(chCard);
			/* this is the atomic number on the C scale */
			punch.punarg[0][punch.npunch] = (float)j;

			/* this will become the ionization stage on C scale */
			ipFFmt = 5;
			punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL) - 1;
			if( lgEOL )
				NoNumb( chCard );
			if( punch.punarg[1][punch.npunch]<0 || punch.punarg[1][punch.npunch]> j+1 )
			{
				fprintf(ioQQQ,"Bad ionization stage - please check Hazy.\nSorry.\n");
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
		else
		{
			/* no element - so make table of all rates */
			strcpy( punch.chPunch[punch.npunch], "GAMt" );
		}

	}
	else if( lgMatch("GRAI",chCard) )
	{
		/* punch grain ... options */
		if( lgMatch("OPAC",chCard) )
		{
			/* check for keyword UNITS on line, then scan wavelength or energy units, 
			 * sets punch.chConPunEnr*/
			ChkUnits(chCard);

			strcpy( punch.chPunch[punch.npunch], "DUSO" );
			/* punch grain opacity command in twice, here and above in opacity */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#grain ext, abs, scat\n" );
		}
		else if( lgMatch("PHYS",chCard) )
		{
			/* punch grain physical conditions */
			strcpy( punch.chPunch[punch.npunch], "DUSP" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#grain depth temp, charge, drift, frc heat, cool\n" );
		}
		else if( lgMatch(" QS ",chCard) )
		{
			strcpy( punch.chPunch[punch.npunch], "DUSQ" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#grain abs, scat, qs\n" );
		}
		else if( lgMatch("TEMP",chCard) )
		{
			/* punch temperatures of each grain species */
			strcpy( punch.chPunch[punch.npunch], "DUST" );
			/* cannot punch grain labels since they are not known yet */
		}
		else if( lgMatch("DRIF",chCard) )
		{
			/* punch drift velocity of each grain species */
			strcpy( punch.chPunch[punch.npunch], "DUSV" );
			/* cannot punch grain labels since they are not known yet */
		}
		else if( lgMatch("EXTI",chCard) )
		{
			/* punch grain extinction */
			strcpy( punch.chPunch[punch.npunch], "DUSE" );
			/* cannot punch grain labels since they are not known yet */
			fprintf( punch.ipPnunit[punch.npunch], 
				"#depth\tA_V(extended)\tA_V(point)\n" );
		}
		else if( lgMatch("CHAR",chCard) )
		{
			/* punch charge per grain (# elec/grain) for each grain species */
			strcpy( punch.chPunch[punch.npunch], "DUSC" );
			/* cannot punch grain labels since they are not known yet */
		}
		else if( lgMatch("HEAT",chCard) )
		{
			/* punch heating due to each grain species */
			strcpy( punch.chPunch[punch.npunch], "DUSH" );
			/* cannot punch grain labels since they are not known yet */
		}
		else if( lgMatch("POTE",chCard) )
		{
			/* punch floating potential of each grain species */
			strcpy( punch.chPunch[punch.npunch], "DUSP" );
			/* cannot punch grain labels since they are not known yet */
		}
		else if( lgMatch("H2RA",chCard) )
		{
			/* punch grain H2rate - H2 formation rate for each type of grains */
			strcpy( punch.chPunch[punch.npunch], "DUSR" );
			/* cannot punch grain labels since they are not known yet */
		}
		else
		{
			fprintf( ioQQQ, "There must be a second key; see Hazy.\nSorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	else if( lgMatch("GAUN",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], "GAUN" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#Gaunt factors.\n" );
	}
	else if( lgMatch( "HIST" , chCard ) )
	{
		/* punch pressure history of current zone */
		if( lgMatch( "PRES" , chCard ) )
		{
			/* punch pressure history - density - pressure for this zone */
			strcpy( punch.chPunch[punch.npunch], "HISp" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#itr zon\tdensity\tpres cur\tpres corr\n" );
		}
		/* punch temperature history of current zone */
		else if( lgMatch( "TEMP" , chCard ) )
		{
			/* punch pressure history - density - pressure for this zone */
			strcpy( punch.chPunch[punch.npunch], "HISt" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#itr zon\ttemperature\theating\tcooling\n" );
		}
	}

	else if( lgMatch("HTWO",chCard) )
	{
		punch.ipPunH2 = punch.ipPnunit[punch.npunch];
		punch.lgPunH2 = TRUE;
		fprintf( punch.ipPnunit[punch.npunch], 
			"#H2 creation, destruction. \n" );
		/* this does not count as a punch option (really) */
		punch.ipPnunit[punch.npunch] = NULL;
		punch.npunch -= 1;
		/* all punch files have the "no clobber" option */
		if( lgMatch("NO C",chCard) && lgMatch("CLOB",chCard) )
			lgPunH2_noclobber = TRUE;

	}

	/* QHEAT has to come before HEAT... */
	else if( lgMatch("QHEA",chCard) )
	{
		gv.QHPunchFile = punch.ipPnunit[punch.npunch];
		gv.lgQHPunLast = punch.lgPunLstIter[punch.npunch];
		fprintf( gv.QHPunchFile, 
			"#Probability distributions from quantum heating routine.\n" );
		/* this does not count as a punch option (really) */
		punch.ipPnunit[punch.npunch] = NULL;
		punch.npunch -= 1;
		/* all punch files have the "no clobber" option */
		if( lgMatch("NO C",chCard) && lgMatch("CLOB",chCard) )
			QHPunchFile_noclobber = TRUE;
	}

	else if( lgMatch("HEAT",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], "HEAT" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\thtot\tctot\theat fracs\n" );
	}

	else if( lgMatch("HUMM",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], "HUMM" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#input to DHs routine.\n" );
	}

	else if( lgMatch("HYDR",chCard) )
	{
		/* punch hydrogen physical conditions */
		if( lgMatch("COND",chCard) )
		{
			strcpy( punch.chPunch[punch.npunch], "HYDc" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#depth\tTe\tHDEN\tEDEN\tHI/H\tHII/H\tH2/H\tH2+/H\tH3+/H\tH-/H\n" );
			/* punch hydrogen ionization */
		}

		/* punch information on 21 cm excitation processes */
		else if( lgMatch("21CM",chCard) )
		{
			/* punch information about 21 cm line */
			strcpy( punch.chPunch[punch.npunch], "21cm" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#depth\tT(spin)\tT(kin)\tT(Lya)\tnLo\tnHi\tOccLya\ttau(21cm)\ttau(Lya)\n" );
		}

		else if( lgMatch("IONI",chCard) )
		{
			strcpy( punch.chPunch[punch.npunch], "HYDi" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#hydr\tzn\tgam1\tcoll ion1\tRecTot\tHRecCaB\thii/hi\tSim hii/hi\threc(esc)\tdec2grd\texc pht\texc col\trec eff\tsec ion\n" );
			/* punch hydrogen populations */
		}
		else if( lgMatch("POPU",chCard) )
		{
			strcpy( punch.chPunch[punch.npunch], "HYDp" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#depth\tHI\tHII\tHN                                                   2s       2p\n" );
		}
		else if( lgMatch("LINE",chCard) )
		{
			/* hydrogen line intensities and optical depths  */
			strcpy( punch.chPunch[punch.npunch], "HYDl" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#nup\tnlo\tE(ryd)\ttau\n" );
		}
		else
		{
			fprintf( ioQQQ, "Punch hydrogen has options: CONDitions POPUlations, and IONIzation.\nSorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	else if( lgMatch("IONI",chCard) )
	{
		if( lgMatch("RATE",chCard) )
		{
			/* punch ionization rates, search for the name of an element */
			if( (nelem = GetElem( chCard ) ) < 0 )
			{
				fprintf( ioQQQ, "There must be an element name on the ionization rates command.  Sorry.\n" );
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}
			punch.punarg[0][punch.npunch] = (float)nelem;
			strcpy( punch.chPunch[punch.npunch], "IONR" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#%s depth\teden\tdynamics.Rate\tabund\tTotIonize\tTotRecom\tSource\t ... \n",
				elementnames.chElementSym[nelem]);
		}
		else
		{
			/* punch table giving ionization means */
			strcpy( punch.chPunch[punch.npunch], "IONI" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Mean ionization distribution\n" );
		}
	}

	else if( lgMatch(" IP ",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], " IP " );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#ionization potentials, valence shell\n" );
	}

	else if( lgMatch("LEID",chCard) )
	{
		if( lgMatch( "LINE" , chCard ) )
		{
			/* final intensities of the Leiden PDR models */
			strcpy( punch.chPunch[punch.npunch], "LEIL" );
			fprintf(punch.ipPnunit[punch.npunch],"#ion\twl\tInt\trel int\n");
		}
		else
		{
			/* structure of the Leiden PDR models */
			strcpy( punch.chPunch[punch.npunch], "LEIS" );
			fprintf( punch.ipPnunit[punch.npunch], 
				/* 1-17 */
				"#Leid  depth\tA_V(extentd)\tA_V(point)\tTe\tH0\tH2\tCo\tC+\tOo\tCO\tO2\tCH\tOH\te\tHe+\tH+\tH3+\t"
				/* 18 - 30 */
				"N(H0)\tN(H2)\tN(Co)\tN(C+)\tN(Oo)\tN(CO)\tN(O2)\tN(CH)\t(OH)\tN(e)\tN(He+)\tN(H+)\tN(H3+)\t"
				/* 31 - 32 */
				"H2(Sol)\tH2(FrmGrn)\t"
				/* 33 - 46*/
				"G0(DB96)\trate(CO)\trate(C)\theat\tcool\tGrnP\tGr-Gas-Cool\tGr-Gas-Heat\tCOds\tH2dH\tH2vH\tChaT\tCR H\tMgI\tSI\t"
				"Si\tFe\tNa\tAl\tC\tC610\tC370\tC157\tC63\tC146\n" );
		}
	}

	else if( lgMatch("LINE",chCard) )
	{
		
		/* check for keyword UNITS on line, then scan wavelength or energy units, 
		 * sets punch.chConPunEnr*/
		/* >>chng 02 jan 28, following had been commented out, and punch lines
		 * failed with insane units.  What happened? */
		ChkUnits(chCard);

		/* options are punch line structure, line intensity, line array,
		 * and line data */
		if( lgMatch("STRU",chCard) ||lgMatch("EMIS",chCard) )
		{
			/* >>chng 02 jul 15 */
			/* this used to be the punch lines structure command, is now
			 * the punch lines emissivity command */
			/* this will be line emissivity vs depth */
			strcpy( punch.chPunch[punch.npunch], "LINS" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Emission line structure.\n" );
			/* read in the list of lines to examine */
			punch_line(punch.ipPnunit[punch.npunch],"READ",FALSE);
		}

		else if( lgMatch("CUMU",chCard) )
		{
			/* this will be integrated line intensity, fcn of depth */
			strcpy( punch.chPunch[punch.npunch], "LINC" );
			/* option for either relative intensity or abs luminosity */
			if( lgMatch("RELA",chCard) )
			{
				lgEOL = TRUE;
				fprintf( punch.ipPnunit[punch.npunch], 
					"#Cumulative emission line relative intensity.\n" );
			}
			else
			{
				fprintf( punch.ipPnunit[punch.npunch], 
					"#Cumulative emission line absolute intensity.\n" );
				lgEOL = FALSE;
			}
			/* read in the list of lines to examine */
			punch_line(punch.ipPnunit[punch.npunch],"READ",lgEOL);
		}

		else if( lgMatch("DATA",chCard) )
		{
			/* punch line data, done in PunchLineData */
			strcpy( punch.chPunch[punch.npunch], "LIND" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Emission line data.\n" );
		}

		else if( lgMatch("ARRA",chCard) )
		{
			/* punch energies and luminosities of transferred lines */
			strcpy( punch.chPunch[punch.npunch], "LINA" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#enr\tintensity\tid\ttype\n" );

		}

		else if( lgMatch("LABE",chCard) )
		{
			/* punch line labels */
			strcpy( punch.chPunch[punch.npunch], "LINL" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#index\tlabel\twavelength\n" );
		}

		else if( lgMatch("OPTI",chCard) )
		{
			/* punch line optical depths, done in PunchLineStuff, located in punchdo.c */
			strcpy( punch.chPunch[punch.npunch], "LINO" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#species\tenergy\topt depth\tdamp\n" );

			/* the default will be to make wavelengths line in the printout, called labels,
			 * if units appears then other units will be used instead */
			strcpy( punch.chConPunEnr[punch.npunch], 
				"labl" );

			/* check for keyword UNITS on line, then scan wavelength or energy units if present,
			 * units are copied into punch.chConPunEnr */
			if( lgMatch("UNIT",chCard) )
				ChkUnits(chCard);

			/* this is optional limit to smallest optical depths */
			ipFFmt = 5;
			punch.punarg[0][punch.npunch] = (float)pow(10.,FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL));
			/* this is default of 0.1 napier */
			if( lgEOL )
			{
				punch.punarg[0][punch.npunch] = 0.1f;
			}
		}

		else if( lgMatch("POPU",chCard) )
		{
			/* punch level populations */
			strcpy( punch.chPunch[punch.npunch], "LINP" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#population information\n" );
			/* this is optional limit to smallest population to punch */
			ipFFmt = 5;
			punch.punarg[0][punch.npunch] = (float)pow(10.,FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL));
			/* this is default of all populations */
			if( lgEOL )
			{
				punch.punarg[0][punch.npunch] = 0.f;
			}
		}

		else if( lgMatch("INTE",chCard) )
		{
			/* this will be full set of line intensities */
			strcpy( punch.chPunch[punch.npunch], "LINI" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Emission line intensities per unit inner area\n" );
			if( lgMatch("COLU",chCard) )
			{
				/* column is key to punch single column */
				strcpy( punch.chPunRltType, "column" );
			}
			else
			{
				/* array is key to punch large array */
				strcpy( punch.chPunRltType, "array " );
			}
			if( lgMatch("EVER",chCard) )
			{
				ipFFmt = 3;
				punch.LinEvery = (long int)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
				punch.lgLinEvery = TRUE;
				if( lgEOL )
				{
					fprintf( ioQQQ, 
						"There must be a second number, the number of zones to print.\nSorry.\n" );
					puts( "[Stop in ParsePunch]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			else
			{
				punch.LinEvery = geometry.nend[0];
				punch.lgLinEvery = FALSE;
			}
		}
		else
		{
			fprintf( ioQQQ, 
				"This option for PUNCH LINE is something that I do not understand.  Sorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	else if( lgMatch("LYMA",chCard) )
	{
		/* lyman alpha lya details */
		strcpy( punch.chPunch[punch.npunch], "LYMA" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tTau\tN2/N1\tTexc\tTe\tTex/Te\tdest\topacity\talbedo\n" );
	}

	else if( lgMatch(" MAP",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], "MAP " );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#te, heating, cooling.\n" );
		/* do cooling space map for specified zones
		 * if no number, or <0, do map and punch out without doing first zone
		 * does map by calling punt(" map") 
		 */
		i = 5;
		map.MapZone = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			map.MapZone = 1;
		}

		/* say output goes to special punch */
		ioMAP = punch.ipPnunit[punch.npunch];

		if( lgMatch("RANG",chCard) )
		{
			int lgLogOn;
			map.RangeMap[0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( map.RangeMap[0] <= 10. )
			{
				map.RangeMap[0] = (float)pow(10.f,map.RangeMap[0]);
				lgLogOn = TRUE;
			}
			else
			{
				lgLogOn = FALSE;
			}

			map.RangeMap[1] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgLogOn )
				map.RangeMap[1] = (float)pow(10.f,map.RangeMap[1]);

			if( lgEOL )
			{
				fprintf( ioQQQ, "There must be a zone number, followed by two temperatures, on this line.  Sorry.\n" );
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}
		
		}
	}

	else if( lgMatch("MOLE",chCard) )
	{
		/* molecules, especially for pdr calculations */
		strcpy( punch.chPunch[punch.npunch], "MOLE" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate");
		
		for(i=0; i<N_H_MOLEC; ++i )
		{
			fprintf(punch.ipPnunit[punch.npunch],"\t%s",hmi.chLab[i] );
		}
		for(i=0; i<NUM_COMOLE_CALC; ++i )
		{
			fprintf(punch.ipPnunit[punch.npunch],"\t%s",co.chLab[i] );
		}
		fprintf(punch.ipPnunit[punch.npunch],"\n");
	}

	else if( lgMatch("OPTICAL",chCard) && lgMatch("DEPTH",chCard) )
	{
		/* check for keyword UNITS on line, then scan wavelength or energy units if present,
		 * units are copied into punch.chConPunEnr */
		ChkUnits(chCard);
		if( lgMatch("FINE",chCard) )
		{
			/* punch fine continuum optical depths */
			rfield.lgPunchOpacityFine = TRUE;
			strcpy( punch.chPunch[punch.npunch], "OPTf" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#energy\tTau tot\topacity\n" );
			ipFFmt = 5;
			/* range option - important since so much data */
			if( lgMatch("RANGE",chCard) ) 
			{
				/* get lower and upper range, must be in Ryd */
				punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
				punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
				if( lgEOL )
				{
					fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
					puts( "[Stop in ParsePunch]" );
					cdEXIT(EXIT_FAILURE);
				}
				if( punch.punarg[0][punch.npunch] >=punch.punarg[1][punch.npunch] )
				{
					fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
					puts( "[Stop in ParsePunch]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			else
			{
				/* these mean full energy range */
				punch.punarg[0][punch.npunch] = 0.;
				punch.punarg[1][punch.npunch] = 0.;
			}
			/* optional last parameter - how many points to bring together */
			punch.punarg[2][punch.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  INPUT_LINE_LENGTH,&lgEOL);
			/* default is to bring together ten */
			if( lgEOL )
				punch.punarg[2][punch.npunch] = 10;
			if( punch.punarg[2][punch.npunch] < 1 )
			{
				fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
				puts( "[Stop in ParsePunch]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
		else
		{
			/* punch coarse continuum optical depths */
			strcpy( punch.chPunch[punch.npunch], "OPTc" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#energy\ttotal\tabsorp\tscat\n" );
		}

	}
	else if( lgMatch(" OTS",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], " OTS" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#otscon, lin, conOpac LinOpc\n" );
	}

	else if( lgMatch("OVER",chCard) && lgMatch(" OVE",chCard) )
	{
		/* punch overview of model results */
		strcpy( punch.chPunch[punch.npunch], "OVER" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tTe\tHtot\thden\teden\t2H_2/H\tHI\tHII\tHeI\tHeII\tHeIII\tCO/C\tC1\tC2\tC3\tC4\tO1\tO2\tO3\tO4\tO5\tO6\tAV(point)\tAV(extend)\n" );
	}

	else if( lgMatch(" PDR",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], " PDR" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tH colden\tTe\tHI/HDEN\tH2/HDEN\tH2*/HDEN\tCI/C\tCO/C\tH2O/O\tG0\tAV(point)\tAV(extend)\tTauV(point)\n" );
	}

	else if( lgMatch("PHYS",chCard) )
	{
		/* punch physical conditions */
		strcpy( punch.chPunch[punch.npunch], "PHYS" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#PhyC    Depth       Te,       HDEN,       Ne,       HTOT,      accel\n" );
	}

	else if( lgMatch("POIN",chCard) )
	{
		/* punch out the pointers */
		punch.ipPoint = punch.ipPnunit[punch.npunch];
		punch.lgPunPoint = TRUE;
		fprintf( punch.ipPnunit[punch.npunch], 
			"#pointers. \n" );
		/* this does not count as a punch option (really) */
		punch.ipPnunit[punch.npunch] = NULL;
		punch.npunch -= 1;
		/* all punch files have the "no clobber" option */
		if( lgMatch("NO C",chCard) && lgMatch("CLOB",chCard) )
			lgPunPoint_noclobber = TRUE;
	}

	else if( lgMatch("PRES",chCard) )
	{
		/* the punch pressure command */
		strcpy( punch.chPunch[punch.npunch], "PRES" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#P depth\tPcorrect\tPcurrent\tPIn+pintg\tpgas(r0)\tpgas\tpram\tprad\tinteg\twindv\tcad\tP(mag)\tconv?\n" );
	}

	else if( lgMatch("RADI",chCard) )
	{
		/* the punch radius command */
		strcpy( punch.chPunch[punch.npunch], "RADI" );
		fprintf( punch.ipPnunit[punch.npunch], "#NZONE\tradius\tr-r0\tdr\n" );
	}

	else if( lgMatch("RECO",chCard) )
	{
		if( lgMatch("COEF",chCard) )
		{
			/* recombination coefficients for everything
			 * punch.ioRecom set to zero in routine zero, non-zero value
			 * is flag to punch recombination coefficients. the output is actually
			 * produced by a series of routines, as they generate the recombination
			 * coefficients.  these include 
			 * dielsupres, helium, hydrorecom, iibod, and makerecomb*/
			punch.ioRecom = punch.ipPnunit[punch.npunch];
			punch.lgioRecom = TRUE;
			fprintf( punch.ipPnunit[punch.npunch], 
				"#recombination coefficients cm3 s-1 for current density and temperature\n" );

			/* this does not count as a punch option (really) */
			punch.ipPnunit[punch.npunch] = NULL;
			punch.npunch -= 1;
			/* all punch files have the "no clobber" option */
			if( lgMatch("NO C",chCard) && lgMatch("CLOB",chCard) )
				lgioRecom_noclobber = TRUE;
		}

		else if( lgMatch("EFFI",chCard) )
		{
			/* punch recombination efficiency */
			strcpy( punch.chPunch[punch.npunch], "RECE" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#Recomn effic H, Heo, He+\n" );
		}

		else
		{
			fprintf( ioQQQ, "No option recognized on this punch recombination command\n" );
			fprintf( ioQQQ, "Valid options are COEFFICIENTS, AGN, and EFFICIENCY\nSorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* punch results command, either as single column or wide array */
	else if( lgMatch("RESU",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], "RESU" );
		if( lgMatch("COLU",chCard) )
		{
			/* column is key to punch single column */
			strcpy( punch.chPunRltType, "column" );
		}
		else
		{
			/* array is key to punch large array */
			strcpy( punch.chPunRltType, "array " );
		}

		/* do not change following, is used as flag in getlines */
		fprintf( punch.ipPnunit[punch.npunch], 
			"#results of calculation\n" );
	}

	else if( lgMatch("SOUR",chCard) )
	{
		if( lgMatch("DEPT",chCard) )
		{
			/* print continuum source function as function of depth */
			strcpy( punch.chPunch[punch.npunch], "SOUD" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#continuum source function vs depth\n" );
		}
		else if( lgMatch("SPEC",chCard) )
		{
			/* print spectrum continuum source function at 1 depth */
			strcpy( punch.chPunch[punch.npunch], "SOUS" );
			fprintf( punch.ipPnunit[punch.npunch], 
				"#continuum source function\n" );
		}
		else
		{
			fprintf( ioQQQ, "A second keyword must appear on this line.\n" );
			fprintf( ioQQQ, "They are DEPTH and SPECTRUM.\n" );
			fprintf( ioQQQ, "Sorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}


	/* punch spectrum the new form of the punch continuum, will eventually replace the standard
	 * punch continuum command */
	else if( lgMatch("SPEC",chCard) && lgMatch("TRUM",chCard) )
	{
		/* this flag is checked in PrtComment to generate a caution
		 * if continuum is punched but iterations not performed */
		punch.lgPunContinuum = TRUE;

		/* set flag for spectrum */
		strcpy( punch.chPunch[punch.npunch], "CONN" );

		fprintf( punch.ipPnunit[punch.npunch], 
			"#Cont Enr\tincid nFn\ttrans\tdiff\tlines \n" );

		/* check for keyword UNITS on line, then scan wavelength or energy units if present,
		 * units are copied into punch.chConPunEnr */
		ChkUnits(chCard);

		/* this points to which punch new continuum this is -
		 * parameters were stored in previous set cpunch commands */
		punch.punarg[0][punch.npunch] = (float)punch.cp_npun;

		++punch.cp_npun;
		/* check that limit not exceeded */
		if( punch.cp_npun > LIMPUN )
		{
			fprintf( ioQQQ, 
				"The limit to the number of PUNCH options is %ld.  Increase LIMPUN in punch.h if more are needed.\nSorry.\n", 
			  LIMPUN );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}

	}

	else if( lgMatch("SPEC",chCard) && lgMatch("CIAL",chCard) )
	{
		/* punch special, will call routine PunSpec */
		strcpy( punch.chPunch[punch.npunch], "SPEC" );
		fprintf( punch.ipPnunit[punch.npunch], "#Special.\n" );
	}

	else if( lgMatch("TEGR",chCard) )
	{
		/* punch tegrid */
		strcpy( punch.chPunch[punch.npunch], "TEGR" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#zone, te, heat, cool.\n" );
	}

	else if( lgMatch("TEMP",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], "TEMP" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tte\tdt/dr\td^2T/dr^2\n" );
	}
	else if( lgMatch("TIME",chCard) )
	{
		/* output the elapsed time per zone */
		strcpy( punch.chPunch[punch.npunch], "TIME" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#zone\tdTime\tElapsed t\n" );
	}


	else if( lgMatch("TPRE",chCard) )
	{
		/* debug output from the temperature predictor in zonestart, 
		 * set with punch tpred command */
		strcpy( punch.chPunch[punch.npunch], "TPRE" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#zone  old temp,  guess Tnew, new temp    delta \n" );
	}

	else if( lgMatch("WIND",chCard) )
	{
		strcpy( punch.chPunch[punch.npunch], "WIND" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#radius\tdepth\tvelocity\tTot accel\tLin accel\tCon accel\tforce mult\n" );
	}

	else
	{
		fprintf( ioQQQ, 
			"ParsePunch cannot find a recognized keyword on this PUNCH command line.\nSorry.\n" );
		puts( "[Stop in ParsePunch]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* increment total number of punch commands, */
	++punch.npunch;

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

/*PunchFilesInit initialize punch file pointers, called from cdInit.
 * NB KEEP THIS ROUTINE SYNCHED UP WITH THE NEXT ONE */
void PunchFilesInit(void)
{
	long int i;
	static int lgFIRST = TRUE;

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

	/* on very first call set lgNoClobber, can be reset with no clobber on punch line */
	if( lgFIRST )
	{
		for( i=0; i < LIMPUN; i++ )
		{
			lgNoClobber[i] = FALSE;
		}
		lgPunConv_noclobber = FALSE;
		lgDROn_noclobber = FALSE;
		lgPunH2_noclobber = FALSE;
		lgPunPoint_noclobber = FALSE;
		lgioRecom_noclobber = FALSE;
		QHPunchFile_noclobber = FALSE;
		lgFIRST = FALSE;
	}

	for( i=0; i < LIMPUN; i++ )
	{
		if( !lgNoClobber[i])
		{
			punch.ipPnunit[i] = NULL;
		}
		/* this sets energy range of continuum, zero says to do full continuum */
		punch.cp_range_min[i] = 0.;
		punch.cp_range_max[i] = 0.;
		/* this means to keep native resolution of code, reset to another
		 * resolution with set cpunch resolution command */
		punch.cp_resolving_power[i] = 0.;
	}

	if( !lgPunConv_noclobber )
	{
		punch.ipPunConv = NULL;
		punch.lgPunConv = FALSE;
	}

	if( !lgDROn_noclobber )
	{
		punch.ipDRout = NULL;
		punch.lgDROn = FALSE;
	}

	if( !lgPunH2_noclobber )
	{
		punch.ipPunH2 = NULL;
		punch.lgPunH2 = FALSE;
	}

	if( !lgPunPoint_noclobber )
	{
		punch.ipPoint = NULL;
		punch.lgPunPoint = FALSE;
	}

	if( !QHPunchFile_noclobber )
		gv.QHPunchFile = NULL;

	if( !lgioRecom_noclobber )
	{
		punch.ioRecom = NULL;
		punch.lgioRecom = FALSE;
	}

	ioMAP = NULL;

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

/*ClosePunchFiles close all punch files. 
 * NB - KEEP THIS ROUTINE SYNCHED UP WITH THE PREVIOUS ONE */
void ClosePunchFiles(void)
{
	long int i;

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

	/* close all punch units cloudy opened with punch command,
	 * lgNoClobber is set true with NO CLOBBER option on punch, says to
	 * not overwrite the files */
	for( i=0; i < punch.npunch; i++ )
	{
		if( punch.ipPnunit[i] != NULL && !lgNoClobber[i])
		{
			fclose( punch.ipPnunit[i] );
			punch.ipPnunit[i] = NULL;
		}
	}

	if( punch.ipPunConv != NULL && !lgPunConv_noclobber )
	{
		fclose( punch.ipPunConv );
		punch.ipPunConv = NULL;
		punch.lgPunConv = FALSE;
	}
	if( punch.ipDRout != NULL && !lgDROn_noclobber  )
	{
		fclose( punch.ipDRout );
		punch.ipDRout = NULL;
		punch.lgDROn = FALSE;
	}
	if( punch.ipPunH2 != NULL && !lgPunH2_noclobber  )
	{
		fclose( punch.ipPunH2 );
		punch.ipPunH2 = NULL;
		punch.lgPunH2 = FALSE;
	}
	if( punch.ipPoint != NULL && !lgPunPoint_noclobber  )
	{
		fclose( punch.ipPoint );
		punch.ipPoint = NULL;
		punch.lgPunPoint = FALSE;
	}
	if( gv.QHPunchFile != NULL  && !QHPunchFile_noclobber )
	{
		fclose( gv.QHPunchFile );
		gv.QHPunchFile = NULL;
	}
	if( punch.ioRecom != NULL  && !lgioRecom_noclobber )
	{
		fclose( punch.ioRecom );
		punch.ioRecom = NULL;
		punch.lgioRecom = FALSE;
	}
	/* this file was already closed as a punch.ipPnunit, erase copy of pointer as well */
	ioMAP = NULL;

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

/*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present,
 * units are copied into punch.chConPunEnr - when doing output, the routine call
 * AnuUnit( energy ) will automatically return the energy in the right units,
 * when called to do punch output */
static void ChkUnits( char *chCard)
{

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

	/* option to set units for continuum energy in punch output */
	if( lgMatch("NITS",chCard) )
	{
		if( lgMatch("MICR",chCard) )
		{
			/* microns */
			strcpy( punch.chConPunEnr[punch.npunch], "micr" );
		}
		else if( lgMatch(" KEV",chCard) )
		{
			/* keV */
			strcpy( punch.chConPunEnr[punch.npunch], " kev" );
		}
		else if( lgMatch("CENT",chCard) )
		{
			/* centimeters*/
			strcpy( punch.chConPunEnr[punch.npunch], "cent" );
		}
		else if( lgMatch(" EV ",chCard) )
		{
			/* eV */
			strcpy( punch.chConPunEnr[punch.npunch], " ev " );
		}
		else if( lgMatch("ANGS",chCard) )
		{
			/* angstroms */
			strcpy( punch.chConPunEnr[punch.npunch], "angs" );
		}
		else if( lgMatch("WAVE",chCard) )
		{
			/* wavenumbers */
			strcpy( punch.chConPunEnr[punch.npunch], "wave" );
		}
		else if( lgMatch(" RYD",chCard) )
		{
			/* the noble rydberg */
			/* >>chng 01 sep 02, had been rdyb unlike proper ryd */
			strcpy( punch.chConPunEnr[punch.npunch], "ryd " );
		}
		else
		{
			fprintf( ioQQQ, "I did not recognize units on this line.\n" );
			fprintf( ioQQQ, "Units are keV, eV, Angstroms, Rydbergs, centimeters, and microns.\nSorry.\n" );
			puts( "[Stop in ParsePunch]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	else
	{
		strcpy( punch.chConPunEnr[punch.npunch], "ryd " );
	}

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