/*ParsePunch parse the punch command */
#include <string.h>
#include "cddefines.h"
#include "elementnames.h"
#include "pnunit.h"
#include "input.h"
#include "ipdr.h"
#include "iph2.h"
#include "levery.h"
#include "zonecnt.h"
#include "ioreco.h"
#include "optype.h"
#include "mappar.h"
#include "pop371.h"
#include "date.h"
#include "punconv.h"
#include "punpoint.h"
#include "punrlttype.h"
#include "getelem.h"
#include "getquote.h"
#include "ffmtread.h"
#include "lgmatch.h"
#include "punlin.h"
#include "physconst.h"
#include "parse.h"

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

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

	/* on first call to this routine in this coreload, initialized array
	 * to false to say that these files are not yet opened */
	if( lgFIRST )
	{
		for( j=0; j<LIMPUN; ++j )
		{
			pnunit.lgPunContinue[j] = FALSE;
			pnunit.lgUnitOpen[j] = FALSE;
			pnunit.ipPnunit[j] = NULL;
		}
		lgFIRST = FALSE;
	}

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

	/* LAST keyword is an option to only punch only on last iteration */
	if( lgMatch("LAST",chCard) )
	{
		pnunit.lgPunLstIter[pnunit.npunch] = TRUE;
	}
	else
	{
		pnunit.lgPunLstIter[pnunit.npunch] = FALSE;
	}

	/* NOCLOSE keyword is an option to keep punching into the same file */
	if( lgMatch("NOCL",chCard) )
	{
		pnunit.lgPunContinue[pnunit.npunch] = TRUE;
	}
	else
	{
		pnunit.lgPunContinue[pnunit.npunch] = FALSE;
	}

	/* 
	 * 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
	 */
	GetQuote( chLabel , chCard );


	/* only open the file if continue option is used and
	 * it has not already been opened during a previous call */
	if( !pnunit.lgUnitOpen[pnunit.npunch] )
	{
		/* pnunit.npunch is set to 0 in ZERO, 
		 * this increment prevents overwriting punch files
		 * pnunit.ipPnunit is the file pointer  */

		/* we now have the  file name, let's open it */
		/* open the file with the name generated above */
		pnunit.ipPnunit[pnunit.npunch] = fopen( chLabel,"w" );
		if( NULL == pnunit.ipPnunit[pnunit.npunch] )
		{
			fprintf( ioQQQ, " ParsePunch could not open punch file %s for writing.\n",
				chLabel);
			puts( "[Stop in getpunch]" );
			exit(1);
		}
		pnunit.lgUnitOpen[pnunit.npunch] = TRUE;
	}


	/* put version number and title of model on output file */
	fprintf( pnunit.ipPnunit[pnunit.npunch], 
		"%s %s\n", 
	  date.chVersion, input.chTitle );

	/* this must come first since elements appear as sub-keywords */
	if( lgMatch("OPAC",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "OPAC" );
		if( lgMatch("TOTA",chCard) )
		{
			/* DoPunch will call popac to parse the subcommands
			 * punch total opacity */
			strcpy( optype.chOpcTyp, "TOTL" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				"     nu      Tot opac  Abs opac  Scat opac Albedo\n" );
		}

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

		else if( lgMatch("GRAI",chCard) )
		{
			/* punch grain command, give optical properties of gains in calculation */
			strcpy( pnunit.chPunch[pnunit.npunch], "DUSO" );
			/* punch grain opacity command in twice, here and above in opacity
			 * chDst = 'OPAC' */
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				"   grain nu,  total,  absorpt,  scat\n" );
		}

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

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

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

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

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

			if( lgEOL )
			{
				fprintf( ioQQQ, " there must be io unit, atom wght, ion, shell\n" );
				puts( "[Stop in getpunch]" );
				exit(1);
			}
			fprintf( pnunit.ipPnunit[pnunit.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 */
			lgHit = FALSE;
			for( nelem=0; nelem<LIMELM; ++nelem )
			{
				/* look for the name of one of the 30 elements */
				if( lgMatch(elementnames.chElementNameShort[nelem],chCard) )
				{
					strcpy( optype.chOpcTyp, elementnames.chElementNameShort[nelem] );
					lgHit = TRUE;
					break;
				}
			}

			/* lgHit was set true if an element was found, error if not */
			if( !lgHit )
			{
				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 getpunch]" );
				exit(1);
			}
		}
	}

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

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

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

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

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

	else if( lgMatch("COOL",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "COOL" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " heat, cool, frac coolant \n" );
	}

	else if( lgMatch("CONT",chCard) )
	{
		/* option to set units for continuum energy in punch output */
		if( lgMatch("NITS",chCard) )
		{
			if( lgMatch("MICR",chCard) )
			{
				strcpy( pnunit.chConPunEnr[pnunit.npunch], "micr" );
			}
			else if( lgMatch(" KEV",chCard) )
			{
				strcpy( pnunit.chConPunEnr[pnunit.npunch], " kev" );
			}
			else if( lgMatch(" EV ",chCard) )
			{
				strcpy( pnunit.chConPunEnr[pnunit.npunch], " ev " );
			}
			else if( lgMatch("ANGS",chCard) )
			{
				strcpy( pnunit.chConPunEnr[pnunit.npunch], "angs" );
			}
			else if( lgMatch("RYDB",chCard) )
			{
				strcpy( pnunit.chConPunEnr[pnunit.npunch], "rydb" );
			}
			else
			{
				fprintf( ioQQQ, " I did not recognize units on this line.\n" );
				fprintf( ioQQQ, " Units are keV, eV, Angstroms, microns\n" );
				puts( "[Stop in getpunch]" );
				exit(1);
			}
		}
		else
		{
			strcpy( pnunit.chConPunEnr[pnunit.npunch], "ryd " );
		}

		if( lgMatch("BINS",chCard) )
		{
			/* continuum binning */
			strcpy( pnunit.chPunch[pnunit.npunch], "CONB" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Continuum binning, energy, width of cells\n" );
		}

		else if( lgMatch("DIFF",chCard) )
		{
			/* diffuse continuum */
			strcpy( pnunit.chPunch[pnunit.npunch], "COND" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Diffuse Continuum\n" );
		}

		else if( lgMatch("EMIT",chCard) )
		{
			/* emitted continuum */
			strcpy( pnunit.chPunch[pnunit.npunch], "CONE" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				"     Enr      reflec     outward     total line cont\n" );
		}

		else if( lgMatch("INCI",chCard) )
		{
			/* incident continuum */
			strcpy( pnunit.chPunch[pnunit.npunch], "CONC" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Incident Continuum, Enr nFn \n" );
		}

		else if( lgMatch("INTE",chCard) )
		{
			/* continuum interactions */
			strcpy( pnunit.chPunch[pnunit.npunch], "CONi" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Continuum interactions, inc, otslin. otscon, outcon, outlin \n" );
			/* this is option for lowest energy, if nothing then zero */
			ipFFmt = 3;
			pnunit.punarg[0][pnunit.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  76,&lgEOL);
		}

		else if( lgMatch("IONI",chCard) )
		{
			/* punch ionizing continuum:
			 * this is option for lowest energy, if nothing then zero */
			ipFFmt = 3;
			pnunit.punarg[0][pnunit.npunch] = (float)FFmtRead(chCard,&ipFFmt,
			  76,&lgEOL);

			/* this is option for smallest interaction to punch, def is 1 percent */
			pnunit.punarg[1][pnunit.npunch] = (float)FFmtRead(chCard,&ipFFmt,76,&lgEOL);
			if( lgEOL )
			{
				pnunit.punarg[1][pnunit.npunch] = 0.01f;
			}
			strcpy( pnunit.chPunch[pnunit.npunch], "CONI" );
			/* put the header at the top of the file */
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				"IonizCont nu Totl flux  flx*cs,  Finc otsl otsc out  rate/tot, integral  line cont\n" );
		}

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

		else if( lgMatch("TRAN",chCard) )
		{
			/* transmitted continuum */
			strcpy( pnunit.chPunch[pnunit.npunch], "CONT" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Transmitted Continuum\n" );
		}

		else if( lgMatch(" RAW",chCard) )
		{
			/* "raw" continua */
			strcpy( pnunit.chPunch[pnunit.npunch], "CORA" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				"Raw Con anu    flux,     otslin     otscon,    refcon     outcon     outlin     ContDiff \n" );
		}

		else if( lgMatch("REFL",chCard) )
		{
			/* reflected continuum */
			strcpy( pnunit.chPunch[pnunit.npunch], "CONR" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Reflected      cont       line       total        albedo    ConID\n" );
		}

		else
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "CON " );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Cont nu    incident   trans    DiffOut  net trans   reflc     total   line cont\n" );
		}

	}
	else if( lgMatch("CONV",chCard) )
	{
		/* punch reason model declared not converged
		 * not a true punch command, since done elsewhere */
		PunConv.ipPunConv = pnunit.ipPnunit[pnunit.npunch];
		PunConv.lgPunConv = TRUE;
		fprintf( PunConv.ipPunConv, " reason for continued iterations\n" );
		/* this does not count as a punch option (really) */
		pnunit.npunch -= 1;
	}

	else if( lgMatch(" DR ",chCard) )
	{
		ipdr.ipDRout = pnunit.ipPnunit[pnunit.npunch];
		ipdr.lgDROn = TRUE;
		/* set punch last flag to whatever it was above */
		ipdr.lgDRPLst = pnunit.lgPunLstIter[pnunit.npunch];
		fprintf( ipdr.ipDRout, " Choice of zone thickness. \n" );
		/* this does not count as a punch option (really) */
		pnunit.lgUnitOpen[pnunit.npunch] = FALSE;
		pnunit.npunch -= 1;
	}

	else if( lgMatch("ELEM",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "ELEM" );
		GetElem(chCard,&j);
		/* this will be the atomic number */
		pnunit.punarg[0][pnunit.npunch] = (float)j;
		fprintf( pnunit.ipPnunit[pnunit.npunch], " element structure, nelem=%3ld\n", 
		  j );
	}

	else if( lgMatch("FEII",chCard) )
	{
		if( lgMatch("FRED",chCard) )
		{
			/* punch out some stuff for Fred's FeII atom */
			fprintf( pnunit.ipPnunit[pnunit.npunch], " FeII rates \n" );
			strcpy( pnunit.chPunch[pnunit.npunch], "FE2f" );
		}
		else if( lgMatch("DEPA",chCard) )
		{
			/* punch out departure coef for the large FeII atom */
			fprintf( pnunit.ipPnunit[pnunit.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( pnunit.chPunch[pnunit.npunch], "FE2D" );
			}
			else
			{
				/* punch a very few selected levels, calls routine FeIIPunDepart */
				strcpy( pnunit.chPunch[pnunit.npunch], "FE2d" );
			}
		}
		else if( lgMatch("POPU",chCard) )
		{
			/* punch out level populations for the large FeII atom */
			fprintf( pnunit.ipPnunit[pnunit.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( pnunit.chPunch[pnunit.npunch], "FE2P" );
			}
			else
			{
				/* punch a very few selected levels, calls routine FeIIPunDepart */
				strcpy( pnunit.chPunch[pnunit.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 getpunch]" );
			exit(1);
		}
	}

	else if( lgMatch("GAMM",chCard) )
	{
		/* punch all photoionizaiton rates for all subshells */
		fprintf( pnunit.ipPnunit[pnunit.npunch], " Photoionization rates \n" );
		strcpy( pnunit.chPunch[pnunit.npunch], "GAMM" );

	}
	else if( lgMatch("GRAI",chCard) )
	{
		/* punch grains
		 * chDst = 'NONE' */
		if( lgMatch("OPAC",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "DUSO" );
			/* punch grain opacity command in twice, here and above in opacity
			 * chDst = 'OPAC' */
			fprintf( pnunit.ipPnunit[pnunit.npunch], " grain ext, abs, scat\n" );
		}
		else if( lgMatch("PHYS",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "DUSP" );
			/* chDst = 'PHYS' */
			fprintf( pnunit.ipPnunit[pnunit.npunch], " grain temp, charge, drift, frc heat, cool\n" );
		}
		else if( lgMatch(" QS ",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "DUSQ" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], " grain abs, scat, qs\n" );
		}
		else
		{
			fprintf( ioQQQ, " There must be a second key; they are  OPAC, _QS_, and PHYS\n" );
			puts( "[Stop in getpunch]" );
			exit(1);
		}
	}

	else if( lgMatch("GAUN",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "GAUN" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " Gaunt factors.\n" );
	}

	else if( lgMatch("HTWO",chCard) )
	{
		iph2.ipPunH2 = pnunit.ipPnunit[pnunit.npunch];
		iph2.lgPunH2 = TRUE;
		fprintf( pnunit.ipPnunit[pnunit.npunch], " H2 creation, destruction. \n" );
		/* this does not count as a punch option (really) */
		pnunit.npunch -= 1;

	}
	else if( lgMatch("HEAT",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "HEAT" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " htot, ctot, heat fracs\n" );
	}

	/* information about the helium atom */
	else if( lgMatch("HELI",chCard) )
	{
		/* singlet populations and rates */
		if( lgMatch("SING",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "HEL1" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], "he1    level,   simple,  clhe2ov1  smhe2ov1   rec,     ioniz,    hei3 \n" );
		}

		/* triplet populations and rates */
		else if( lgMatch("TRIP",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "HEL3" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], "HeTrip  level,    simple,    rec,     ioniz,    coll \n" );
		}
		else
		{
			fprintf( ioQQQ, " did not find a keyword on this  line.  Keys are SING, _ION, and TRIP.\n" );
			puts( "[Stop in getpunch]" );
			exit(1);
		}
	}

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

	else if( lgMatch("HYDR",chCard) )
	{
		/* punch hydrogen physical conditions */
		if( lgMatch("COND",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "HYDc" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], " Hyd, DR,     Te,     HDEN,    EDEN,   HI/H,    HII/H,   H2/H,     H2+/H,   H3+/H,   H-/H\n" );
			/* punch hydrogen ionization */
		}
		else if( lgMatch("IONI",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "HYDi" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], "   gam1    RecEff   HRecRad    hii/hi Shii/hi   hrec(13) dec grd  GamEff   ColEff\n" );
			/* punch hydrogen populations */
		}
		else if( lgMatch("POPU",chCard) )
		{
			strcpy( pnunit.chPunch[pnunit.npunch], "HYDp" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], "   depth      HI      HII       HN                                                   2s       2p\n" );
		}
		else
		{
			fprintf( ioQQQ, " Punch hydrogen has options: CONDitions POPUlations and IONIzation.\n" );
			puts( "[Stop in getpunch]" );
			exit(1);
		}
	}

	else if( lgMatch("IONI",chCard) )
	{
		/* punch table giving ionization means */
		strcpy( pnunit.chPunch[pnunit.npunch], "IONI" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], 
			" Mean ionization distribution\n" );
	}

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

	else if( lgMatch("LINE",chCard) )
	{
		/* options are punch line structure, line intensity, line array,
		 * and line data */
		if( lgMatch("STRU",chCard) )
		{
			/* this will be line structure, emissivity vs depth */
			strcpy( pnunit.chPunch[pnunit.npunch], "LINS" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" Emission line structure.\n" );
			/* read in the list of lines to examine */
			punlin(pnunit.ipPnunit[pnunit.npunch],"READ",FALSE);
		}

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

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

		else if( lgMatch("ARRA",chCard) )
		{
			/* punch energies and luminosities of transferred lines */
			strcpy( pnunit.chPunch[pnunit.npunch], "LINA" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], " Emission line array\n" );
		}

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

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

	else if( lgMatch("LYMA",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "LYMA" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " depth Tau, N2/N1, Texc, Te, Tex/Te dest ots\n" );
	}

	else if( lgMatch(" MAP",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "MAP " );
		fprintf( pnunit.ipPnunit[pnunit.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 */
		i = 5;
		MapPar.MapZone = (long)FFmtRead(chCard,&i,76,&lgEOL);
		if( lgEOL )
		{
			MapPar.MapZone = 1;
		}

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

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

			MapPar.RangeMap[1] = (float)FFmtRead(chCard,&i,76,&lgEOL);
			if( lgLogOn )
				MapPar.RangeMap[1] = (float)pow(10.,MapPar.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]" );
				exit(1);
			}
		
		}
	}

	else if( lgMatch("MOLE",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "MOLE" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " dr, h2/h CO/H, H2O/H, OH/H, CH/HO2/H\n" );
	}

	else if( lgMatch("OPTI",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "OPTI" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " total, absorp, scat taus\n" );

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

	else if( lgMatch("OVER",chCard) )
	{
		/* punch overview of model results */
		strcpy( pnunit.chPunch[pnunit.npunch], "OVER" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], 
			"  depth     Te     Htot    hden   eden   HI   HII   HeI  HeII  HeIII   C1    C2    C3    C4    O1  O2    O3    O4    O5    O6\n" );
	}

	else if( lgMatch(" PDR",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], " PDR" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], 
			" depth    H colden    AV        Te      HI/HDEN   H2/HDEN   CI/C      CO/C       H2O/O\n" );
	}

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

	else if( lgMatch("POIN",chCard) )
	{
		/* punch out the pointers */
		PunPoint.ipPoint = pnunit.ipPnunit[pnunit.npunch];
		PunPoint.lgPunPoint = TRUE;
		fprintf( pnunit.ipPnunit[pnunit.npunch], " pointers. \n" );
		/* this does not count as a punch option (really) */
		pnunit.npunch -= 1;
	}

	else if( lgMatch("PRES",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "PRES" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], 
			" P depth,   p(G+L)    PIn+pintg  pgas(r0)   pgas     pram       prad      integ,    windv  conv?\n" );
	}

	else if( lgMatch("RADI",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "RADI" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " NZONE, radius, r-r0, dr.\n" );
	}

	else if( lgMatch("RECO",chCard) )
	{
		if( lgMatch("COEF",chCard) )
		{
			/* recombination coefficients for everything
			 * 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*/
			ioRecom = pnunit.ipPnunit[pnunit.npunch];
			lgioRecom = TRUE;
			fprintf( pnunit.ipPnunit[pnunit.npunch], 
				" recombination coefficients cm3 s-1 for current density and temperature\n" );

			/* this does not count as a punch option (really) */
			pnunit.npunch -= 1;
		}

		else if( lgMatch("EFFI",chCard) )
		{
			/* punch recombination efficiency */
			strcpy( pnunit.chPunch[pnunit.npunch], "RECE" );
			fprintf( pnunit.ipPnunit[pnunit.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 and EFFICIENCY\n" );
			puts( "[Stop in getpunch]" );
			exit(1);
		}
	}

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

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

	else if( lgMatch("SOUR",chCard) )
	{
		if( lgMatch("DEPT",chCard) )
		{
			/* print continuum source function as function of depth */
			strcpy( pnunit.chPunch[pnunit.npunch], "SOUD" );
			fprintf( pnunit.ipPnunit[pnunit.npunch], " continuum source function vs depth\n" );
		}
		else if( lgMatch("SPEC",chCard) )
		{
			/* print spectrum continuum source function at 1 depth */
			strcpy( pnunit.chPunch[pnunit.npunch], "SOUS" );
			fprintf( pnunit.ipPnunit[pnunit.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 getpunch]" );
			exit(1);
		}
	}

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

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

	else if( lgMatch("TEMP",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "TEMP" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " DEPTH, te, dt/dr.\n" );
	}

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

	else if( lgMatch("VERN",chCard) )
	{
		/* punch verner command , three optional parameters, threshold for faintest
		 * line to print, lower and upper energy bounds */

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

		/* 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,76,&lgEOL);
		if( lgEOL )
		{
			FeII.fe2thresh = 0.;
		}

		/* it is a log if negative */
		if( FeII.fe2thresh < 0. )
		{
			FeII.fe2thresh = (float)pow(10.,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,76,&lgEOL);
		if( lgEOL )
		{
			FeII.fe2ener[0] = 0.;
		}

		FeII.fe2ener[1] = (float)FFmtRead(chCard,&i,76,&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.,FeII.fe2ener[0]);
			FeII.fe2ener[1] = (float)pow(10.,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 pop371 file */
		fprintf( pnunit.ipPnunit[pnunit.npunch], " Results from Verner FeII atom.\n" );
	}

	else if( lgMatch("WIND",chCard) )
	{
		strcpy( pnunit.chPunch[pnunit.npunch], "WIND" );
		fprintf( pnunit.ipPnunit[pnunit.npunch], " radius, DEPTH,vel, accel, force mul\n" );
	}

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

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

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