/home66/gary/public_html/cloudy/c08_branch/source/parse_punch.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  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 "cddrive.h"
#include "physconst.h"
#include "elementnames.h"
#include "input.h"
#include "geometry.h"
#include "prt.h"
#include "optimize.h"
#include "rfield.h"
#include "hcmap.h"
#include "atomfeii.h"
#include "h2.h"
#include "mole.h"
#include "hmi.h"
#include "version.h"
#include "grainvar.h"
#include "parse.h"
#include "grid.h"
#include "punch.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 - default was to overwrite,
 * >>chng 07 feb 02 invert logic, default not to overwrite, but clobber
 * keyword says to overwrite */
static bool lgNoClobber[LIMPUN];

/* these are for some special cases, same purpose as previous no clobber */
static bool lgPunConv_noclobber , lgDROn_noclobber , 
        lgPunPoint_noclobber , lgioRecom_noclobber , lgQHPunchFile_noclobber,
        lgTraceConvergeBase_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 */

/* punch file header headers - these are written into the string chHeader when
 * the command is parsed
 * punch.lgPunHeader determines whether header is punched 
 */


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

#define MAX_HEADER_SIZE 2000

        char chHeader[MAX_HEADER_SIZE], 
                chTemp[MAX_HEADER_SIZE];


        DEBUG_ENTRY( "ParsePunch()" );

        /* initialize chHeader string with nonsense, compare at end to see if we have any actual headers. */
        sprintf( chHeader, "ArNdY38dZ9us4N4e12SEcuQ" );

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

        /* if this is first call to this routine we will need to init some vars */
        {
                static bool lgMustInit=true;
                if( lgMustInit )
                {
                        lgMustInit = false;

                        /* following are for punch LineList option
                         * nLineList is number of em lines, -1 if not defined */
                        punch.nLineList = (long int*)MALLOC((size_t)(LIMPUN*sizeof(long int)) );
                        /* lgLineListRatio flag saying to take ratios of pairs of lines */
                        punch.lgLineListRatio = (bool*)MALLOC((size_t)(LIMPUN*sizeof(bool)) );
                        /* wlLineList is set of emission lines for LineList */
                        punch.wlLineList = (realnum **)MALLOC((size_t)(LIMPUN*sizeof(realnum *)) );
                        /* chLineListLabel is label for line list */
                        punch.chLineListLabel = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) );
                        /* remaining dimensions will be created if punch LineList command entered 
                         * but init nLineList to -1 to signal not set */
                        for( i=0; i<LIMPUN; ++i )
                        {
                                punch.nLineList[i] = -1;
                                punch.lgFITS[i] = false;
                        }

                        /* following are parallel but for punch averages option 
                         * nAverageList is number of averages, -1 if not defined */
                        punch.nAverageList = (long int*)MALLOC((size_t)(LIMPUN*sizeof(long int)) );
                        /* chAverageSpeciesLabel is label for species of each average */
                        punch.chAverageSpeciesLabel = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) );
                        /* chAverageType is label for species of type of average */
                        punch.chAverageType = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) );
                        /* nAverageIonList is set of ion stages for average */
                        punch.nAverageIonList = (int **)MALLOC((size_t)(LIMPUN*sizeof(int *)) );
                        /* nAverage2ndPar is second parameters of each average */
                        punch.nAverage2ndPar = (int **)MALLOC((size_t)(LIMPUN*sizeof(int *)) );
                        /* remaining dimensions will be created if punch averages command entered 
                         * but init nLineList to -1 to signal not set */
                        for( i=0; i<LIMPUN; ++i )
                        {
                                punch.nAverageList[i] = -1;
                        }
                }
        }

        /* initialize this flag to false, set true for special cases below */
        punch.lgPunchToSeparateFiles[punch.npunch] = false;

        /* LAST keyword is an option to punch only on last iteration */
        if( nMatch("LAST",chCard) )
                punch.lgPunLstIter[punch.npunch] = true;
        else
                punch.lgPunLstIter[punch.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
         * last parameter says whether to abort if no quote found        */
        if( GetQuote( chLabel , chCard , true ) )
                /* this can't happen since routine would not return at all if no double quotes found */
                TotalInsanity();

        /* 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",
                        chLabel);
                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");
                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 );

        /* there may be a second file name, and we need to get it off the line
         * before we parse options, last false parameter says not to abort if
         * missing - this is not a problem at this stage */
        if( GetQuote( chSecondFilename , chCard , false ) )
                lgSecondFilename = false;
        else
                lgSecondFilename = true;

        /* CLOBBER clobber keyword is an option to overwrite rather than
         * append to a given file */
        if( nMatch("CLOB",chCard) )
        {
                if( nMatch(" NO ",chCard) )
                {
                        /* do not clobber files */
                        lgNoClobber[punch.npunch] = true;
                }
                else
                {
                        /* clobber files */
                        lgNoClobber[punch.npunch] = false;
                }
        }

        /* 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( !nMatch("O TI",chCard) && nMatch("TITL",chCard))
        {
                sprintf( chHeader, 
                        "# %s %s\n", 
                  t_version::Inst().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 no hash appears the hash marks 
         * will not occur */
        if( nMatch("O HA",chCard) )
                punch.lgHashEndIter[punch.npunch] = false;

        /* punch opacity must come first since elements appear as sub-keywords */
        if( nMatch("OPAC",chCard) )
        {
                /* for grid run, punch results of different models to different files. */
                punch.lgPunchToSeparateFiles[punch.npunch] = true;

                /* 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( nMatch("TOTA",chCard) )
                {
                        /* DoPunch will call punch_opacity to parse the subcommands
                         * punch total opacity */
                        strcpy( punch.chOpcTyp[punch.npunch], "TOTL" );
                        sprintf( chHeader, 
                                "#nu\tTot opac\tAbs opac\tScat opac\tAlbedo\telem\n" );
                }

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

                else if( nMatch("FINE",chCard) )
                {
                        /* punch the fine opacity array */
                        rfield.lgPunchOpacityFine = true;
                        strcpy( punch.chOpcTyp[punch.npunch], "FINE" );
                        /* check for keyword UNITS on line, then scan wavelength or energy units if present,
                         * units are copied into punch.chConPunEnr */
                        ChkUnits(chCard);

                        sprintf( chHeader, 
                                "#nu\topac\n" );
                        ipFFmt = 5;
                        /* range option - important since so much data - usually want to
                         * only give portion of the continuum */
                        if( nMatch("RANGE",chCard) ) 
                        {
                                /* get lower and upper range, must be in Ryd */
                                punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
                                punch.punarg[punch.npunch][1] = (realnum)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");
                                        cdEXIT(EXIT_FAILURE);
                                }
                                if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] )
                                {
                                        fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else
                        {
                                /* these mean full energy range */
                                punch.punarg[punch.npunch][0] = 0.;
                                punch.punarg[punch.npunch][1] = 0.;
                        }
                        /* optional last parameter - how many points to bring together */
                        punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt,
                          INPUT_LINE_LENGTH,&lgEOL);

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

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

                else if( nMatch("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 */
                        sprintf( chHeader, 
                                "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" );
                }

                else if( nMatch("BREM",chCard) )
                {
                        /* punch bremsstrahlung opacity */
                        strcpy( punch.chOpcTyp[punch.npunch], "BREM" );
                        sprintf( chHeader, 
                                "#nu\tbrem opac\n" );
                }

                else if( nMatch("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( punch.chOpcTyp[punch.npunch], "SHEL" );

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

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

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

                        if( lgEOL )
                        {
                                fprintf( ioQQQ, "There must be IO unit, atom weight, ion, shell\nSorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        sprintf( chHeader, 
                                "#sub shell cross section\n" );
                }

                else if( nMatch("ELEM",chCard) )
                {
                        /* punch element 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 name on the card is ignored.
                         * The code stops in punch_opacity 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" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* copy string over */
                        strcpy( punch.chOpcTyp[punch.npunch], elementnames.chElementNameShort[nelem] );
                }
                else
                {
                        fprintf( ioQQQ, " I did not recognize a keyword on this punch opacity command.\n" );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

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

        /* punch grain abundance will be handled later */
        else if( nMatch("ABUN",chCard) && !nMatch("GRAI",chCard) )
        {
                /* punch abundances */
                strcpy( punch.chPunch[punch.npunch], "ABUN" );
                sprintf( chHeader, 
                        "#abund H" );
                for(nelem=ipHELIUM;nelem<LIMELM; ++nelem )
                {
                        sprintf( chTemp, "\t%s",
                                elementnames.chElementNameShort[nelem] );
                        strcat( chHeader, chTemp );
                }
                strcat( chHeader, "\n");
        }

        else if( nMatch(" AGE",chCard) )
        {
                /* punch ages */
                strcpy( punch.chPunch[punch.npunch], "AGES" );
                sprintf( chHeader, 
                        "#ages depth\tt(cool)\tt(H2 dest)\tt(CO dest)\tt(OH dest)\tt(H rec)\n" );
        }

        else if( nMatch(" 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( nMatch("CHAR",chCard) )
                {
                        strcpy( punch.chPunch[punch.npunch], "CHAG" );
                        sprintf( chHeader, 
                                "#charge exchange rate coefficnt\n" );
                }

                else if( nMatch("RECO",chCard) )
                {
                        /* punch recombination rates for AGN3 table */
                        strcpy( punch.chPunch[punch.npunch], "RECA" );
                        sprintf( chHeader, 
                                "#Recom rates for AGN3 table\n" );
                }

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

                else if( nMatch("HECS",chCard) )
                {
                        /* create table for appendix in AGN */
                        strcpy( punch.chPunchArgs[punch.npunch], "HECS" );
                        sprintf( chHeader, 
                                "#AGN3 he cs \n" );
                }

                else if( nMatch("HEMI",chCard) )
                {
                        /* HEMIS - continuum emission needed for chap 4 of AGN3 */
                        strcpy( punch.chPunchArgs[punch.npunch], "HEMI" );

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

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

        else if( nMatch("AVER",chCard) )
        {
                /* punch averages */
                strcpy( punch.chPunch[punch.npunch], "AVER" );
                /* no need to print this standard line of explanation*/
                /*sprintf( chHeader, " asserts\n" );*/

                /* actually get the averages from the input stream, and malloc the 
                 * space in the arrays 
                 * punch io unit not used in read */
                punch_average( punch.npunch, "READ", chHeader );
        }

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

        else if( nMatch("CHEM",chCard) )
        {

                if( nMatch( "RATE" , chCard ) )
                {
                        /* >>chng 06 May 30, NPA.  Punch reaction rates for selected species */

                        if( nMatch( " CO " , chCard ) )
                        {

                                strcpy( punch.chPunch[punch.npunch], "RCO " );
                                sprintf( chHeader, 
                                "#Depth\tH2O_C3HP_CO_C2H3P\tC2H2_HCOP_CO_C2H3P\tC2H_HCOP_CO_C2H2P\tO_C3H_CO_C2H\t"
                                "O_C2H2_CO_CH2\tOH_C2H2_CO_CH3\tHCOP_C3_C3HP_CO\tO2_C3_CO_C2_O\tO_C3_CO_C2\t"
                                "H_COP_CO_HP\tHminus_HCOP_CO_H2\tH2P_CO_HCOP_H\tH2P_CO_COP_H2\tH3P_CO_HCOP_H2\t"
                                "HeP_CO_OP_C_He\tHeP_CO_O_CP_He\tcrnu_CO_O_C\tCRP_CO_COP_e\tnu_CO_O_C\te_HCOP_CO_H\t"
                                "C_COP_CO_CP\tC_HCOP_CO_CHP\tC_O_CO_nu\tC_O2_CO_O\tC_OH_CO_H\tC_SiOP_SiP_CO\tCP_OH_CO_HP\t"
                                "CP_SiO_SiP_CO\tCP_O2_CO_OP\tO_CH_CO_H\tO_CH2_CO_H_H\tO_CH2_CO_H2\tO_COP_CO_OP\tOP_CO_COP_O\t"
                                "CH_COP_CO_CHP\tCH_HCOP_CO_CH2P\tCH_O2_CO_OH\tCH2_COP_CO_CH2P\tCH2_HCOP_CO_CH3P\t"
                                "CH2_O2_CO_H2O\tCOP_O2_O2P_CO\tH2O_COP_CO_H2OP\tH2O_HCOP_CO_H3OP\tH2OP_CO_HCOP_OH\t"
                                "HCOP_SiH_SiH2P_CO\tHCOP_SiO_SiOHP_CO\tOH_COP_CO_OHP\tOH_HCOP_CO_H2OP\tOHP_CO_HCOP_O\t"
                                "COP_CH4_CO_CH4P\tCO_CH4P_HCOP_CH3\tCO_CH5P_HCOP_CH4\tHP_OCS_HSP_CO\tHeP_OCS_SP_CO_He\t"
                                "OCNP_e_CO_N\tOCSP_e_S_CO\tOCS_NU_S_CO\tOCS_CRP_S_CO\tC_NO_CO_N\tC_OCN_CO_CN\t" 
                                "C_SO_S_CO\tO_CN_CO_N\tO_HCN_CO_NH\tO_OCN_NO_CO\tO_CS_S_CO\tO_OCS_SO_CO\tOH_HCN_CO_NH2\t"   
                                "CN_NO_N2_CO\tCN_O2_NO_CO\tCO_HS_OCS_H\tCP_SO_SP_CO\tCP_OCS_CSP_CO\tCHP_OCS_HCSP_CO\t" 
                                "NP_CO_NOP_C\tNP_OCS_SP_CO_N\tNHP_CO_HCOP_N\tNHP_CO_OCNP_H\tNH_HCOP_CO_NH2P\t"    
                                "NH2_HCOP_CO_NH3P\tNH3_HCOP_CO_NH4P\tCNP_O2_NOP_CO\tHCNP_CO_HCOP_CN\tCO_HNOP_NO_HCOP\t"    
                                "N2P_OCS_SP_N2_CO\tHCOP_S_HSP_CO\tHCOP_CS_HCSP_CO\tNH_COP_CO_NHP\tNH2_COP_CO_NH2P\t"    
                                "NH3_COP_CO_NH3P\tCNP_CO_COP_CN\tHCN_COP_CO_HCNP\tCO_N2P_N2_COP\tCOP_NO_NOP_CO\t"      
                                "CO_S_OCS_NU\tNP_CO_COP_N\tCOP_S_SP_CO\tC_CO_C2_O\tO_C2_CO_C\tC2_O2_CO_CO\t"
                                "C2_O2P_COP_CO\tC2_COP_CO_C2P\tC2P_O2_COP_CO\tO_C2H_CO_CH\tO_CCl_Cl_CO\t"
                                "CO_H2ClP_HCl_HCOP\tHNC_HCOP_HCNHP_CO\tHCN_HCOP_HCNHP_CO\tC2H_COP_CO_C2HP\tC2_HCOP_CO_C2HP\n");
                        }

                        else if( nMatch( " OH " , chCard ) )
                        {
                                strcpy( punch.chPunch[punch.npunch], "ROH " );
                                sprintf( chHeader, 
                                "#Depth\tnu_H2O_OH_H\tnu_OH_OHP_e\tnu_OH_O_H\tnu_H2O_OH_H\tnu_OH_OHP_e\tnu_OH_O_H\tcrnu_H2O_OH_H\tcrnu_OH_O_H\tCP_OH_CO_HP\tCP_OH_COP_H\t"
                                "CP_OH_CO_HP\tCP_OH_COP_H\tCH2_OHP_OH_CH2P\tCH2P_O2_HCOP_OH\tH2O_COP_HCOP_OH\tH2OP_H2O_H3OP_OH\tC_H2OP_OH_CHP\tOP_OH_O2P_H\tOP_OH_OHP_O\tSiP_OH_SiOP_H\t"
                                "CH_H2OP_OH_CH2P\tCH_OHP_OH_CHP\tCHP_OH_COP_H2\tCHP_O2_COP_OH\tCH2_H2OP_OH_CH3P\tOH_COP_HCOP_O\tOH_H2OP_H3OP_O\tOHP_H2O_H2OP_OH\tOHP_O2_O2P_OH\tOHP_OH_H2OP_O\t"
                                "O_CH4P_OH_CH3P\tOP_CH4_OH_CH3P\tCH5P_OH_H2OP_CH4\tOHP_C2_C2P_OH\tOH_C2P_C2_OHP\tCH_NO_CN_OH\tNH_O_OH_N\tNH_OH_HNO_H\tO_HNO_NO_OH\tNH2_OH_H2O_NH\t"
                                "NH2_NO_N2_OH_H\tOH_CN_OCN_H\tOH_S_HS_O\tOH_S_SO_H\tNHP_OH_H2OP_N\tNHP_H2O_OH_NH2P\tNHP_O2_NOP_OH\tO_HSP_SP_OH\tNH2P_OH_H2OP_NH\tNH2P_H2O_NH3P_OH\t"
                                "NH2_H2OP_NH3P_OH\tNH2P_O2_HNOP_OH\tOH_NH3P_NH4P_O\tOH_HCNP_CN_H2OP\tOH_HNOP_NO_H2OP\tOH_SP_SOP_H\tNH3P_H2O_NH4P_OH\tNH3_H2OP_NH4P_OH\tH2O_CNP_HCNP_OH\tH2OP_S_HSP_OH\t"
                                "NH_OHP_OH_NHP\tNH2_OHP_OH_NH2P\tOHP_NH3_NH3P_OH\tOH_CNP_CN_OHP\tOH_N2P_N2_OHP\tOHP_NO_NOP_OH\tNP_OH_OHP_N\tOHP_S_SP_OH\tH2OP_HNC_HCNHP_OH\tH2OP_HCN_HCNHP_OH\t"
                                "H2O_C2P_C2HP_OH\tH2OP_C2_C2HP_OH\tH2OP_C2H_C2H2P_OH\tOHP_C2H_C2HP_OH\tHminus_OH_H2O_e\tHminus_O_OH_e\tHP_OH_OHP_H\tH2s_OH_O_H2_H\tH2s_H2O_OH_H2_H\tH2P_OH_H2OP_H\t"
                                "H2P_OH_OHP_H2\tH3P_OH_H2OP_H2\tHeP_OH_OP_He_H\tHeP_H2O_OH_He_HP\tH3P_NO2_NOP_OH_H2\tCH_O2_CO_OH\tH2OP_CO_HCOP_OH\tOH_COP_CO_OHP\tOH_HCOP_CO_H2OP\tCH2_OH_H2O_CH\t"
                                "C_OH_O_CH\tO_CH_OH_C\tO_CH2_OH_CH\tO_H2O_OH_OH\tO_OH_O2_H\tSi_OH_SiO_H\tOH_OH_H2O_O\tO_CH4_OH_CH3\tOH_CH2_O_CH3\tOH_CH3_CH4_O\t"
                                "OH_CH3_H2O_CH2\tH2O_CH3_OH_CH4\tOH_CH4_H2O_CH3\tN_OH_O_NH\tN_OH_NO_H\tCH2_NO_HCN_OH\tNH_OH_NH2_O\tNH_OH_H2O_N\tNH_H2O_OH_NH2\tNH_NO_N2_OH\t"
                                "NH_NO2_N2O_OH\tO_NH2_OH_NH\tO_NH3_OH_NH2\tO_HCN_CN_OH\tO_HS_S_OH\tNH2_OH_NH3_O\tOH_NH3_H2O_NH2\tOH_CN_HCN_O\tOH_HCN_CN_H2O\tOH_NO_NO2_H\t"
                                "OH_N2O_HNO_NO\tOH_CS_OCS_H\tNO_HNO_N2O_OH\tO_C2H2_C2H_OH\tOH_C2H2_C2H_H2O\te_H2OP_OH_H\te_H3OP_OH_H_H\te_H3OP_OH_H2\te_SiOHP_Si_OH\tH_OH_O_H_H\t"
                                "H_H2O_OH_H_H\tH_OH_O_H2\tH_H2O_OH_H2\tH_OH_H2O_nu\tHminus_H3OP_OH_H2_H\tH_O_OH_nu\tH2_OH_H2O_H\tH2_O_OH_H\tH2_OH_O_H2_H\tH2_H2O_OH_H2_H\t"
                                "H2_O2_OH_OH\tH_O2_OH_O\tC_OH_CO_H\tOH_HCN_CO_NH2\tOH_C2H2_CO_CH3\tOH_C2H2_CO_CH3\n");
                        }

                        else
                        {
                                fprintf(ioQQQ," The keyword CO or OH must appear on the PUNCH CHEMISTRY RATES command.\n");
                                fprintf( ioQQQ, " Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                }
        }

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

        else if( nMatch("COOL",chCard) )
        {
                /* punch cooling, actually done by routine cool_punch */
                strcpy( punch.chPunch[punch.npunch], "COOL" );
                /*>>chng 06 jun 06, revise to be same as punch cooling */
                sprintf( chHeader, 
                        "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\tcool fracs\n" );
        }

        else if( nMatch("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( nMatch( "ADVE" , chCard ) )
                {
                        /* punch information relating to advection */
                        strcpy( punch.chPunch[punch.npunch], "DYNa");
                        sprintf( chHeader, 
                                "#advection depth\tHtot\tadCool\tadHeat\tdCoolHeatdT\t"
                                "Source[hyd][hyd]\tRate\tEnthalph\tadSpecEnthal\n" );
                }
                else
                {
                        fprintf( ioQQQ, " I did not recognize a sub keyword on this PUNCH DYNAMICS command.\n" );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

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

        else if( nMatch("EXEC",chCard) && nMatch("TIME",chCard) )
        {
                /* output the execution time per zone */
                strcpy( punch.chPunch[punch.npunch], "XTIM" );
                sprintf( chHeader, 
                        "#zone\tdTime\tElapsed t\n" );
        }

        else if( nMatch("FEII",chCard) || nMatch("FE II",chCard) )
        {
                /* something to do with FeII atom - several options
                 * FeII column densities */
                if( nMatch("COLU",chCard) && nMatch("DENS",chCard) )
                {
                        /* punch FeII column density */
                        strcpy( punch.chPunch[punch.npunch], "FEcl" );

                        /* file will give energy, statistical weight, and column density [cm-2] */
                        sprintf( chHeader, 
                                "#FeII: energy\tstat wght\tcol den\n" );
                }

                /* FeII continuum, only valid if large atom is set */
                else if( nMatch("CONT",chCard) )
                {
                        /* punch FeII continuum */
                        strcpy( punch.chPunch[punch.npunch], "FEco" );

                        sprintf( chHeader, 
                                "#FeII: Wl(A)\tInt[erg cm-2 s-1]\n" );
                }

                else if( nMatch("DEPA",chCard) )
                {
                        /* punch out departure coefficient for the large FeII atom */
                        sprintf( chHeader, 
                                "#FeII departure coefficient \n" );
                        /* optional keyword all means do all levels, if not then just do subset */
                        if( nMatch(" 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( nMatch("LEVE",chCard) )
                {
                        /* punch level energies and statistical weights for large FeII atom */
                        sprintf( chHeader, 
                                "#FeII energy(wn)\tstat weight\n" );
                        strcpy( punch.chPunch[punch.npunch], "FE2l" );
                }

                else if( nMatch("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( nMatch("SHOR",chCard) )
                        {
                                FeII.lgShortFe2 = true;
                        }
                        else
                        {
                                FeII.lgShortFe2 = false;
                        }

                        /* first optional number changes the threshold 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 = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                FeII.fe2thresh = 0.;
                        }

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

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

                        FeII.fe2ener[1] = (realnum)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] = (realnum)pow((realnum)10.f,FeII.fe2ener[0]);
                                FeII.fe2ener[1] = (realnum)pow((realnum)10.f,FeII.fe2ener[1]);
                        }

                        /* entered in Ryd in above, convert to wavenumbers */
                        FeII.fe2ener[0] /= (realnum)WAVNRYD;
                        FeII.fe2ener[1] /= (realnum)WAVNRYD;

                        /* these results are actually created by the FeIIPunchLines routine
                        * that lives in the FeIILevelPops file */
                        sprintf( chHeader, 
                                "#FeII ipHi\tipLo\tWL(A)\tlog I\tI/Inorm\t\tTau\n" );
                }

                else if( nMatch("OPTI",chCard) && nMatch("DEPT",chCard) )
                {
                        /* punch optical depths for large FeII atom */
                        sprintf( chHeader, 
                                "#FeII hi\tlow\twl(A)\ttau\n" );
                        strcpy( punch.chPunch[punch.npunch], "FE2o" );
                }

                else if( nMatch("POPU",chCard) )
                {
                        /* punch out level populations for the large FeII atom */
                        sprintf( chHeader, 
                                "#FeII level populations [cm^-3]\n" );

                        /* this is keyword RELATIVE that says to punch relative to total Fe+, 
                         * default is actual density (cm-3) */
                        if( nMatch("RELA",chCard) )
                        {
                                punch.punarg[punch.npunch][2] = 0.;
                        }
                        else
                        {
                                /* default is to punch density (cm-3) */
                                punch.punarg[punch.npunch][2] = 1.;
                        }

                        /* optional keyword all means do all levels, if not then just do subset */
                        if( nMatch(" ALL",chCard) )
                        {
                                /* punch all levels, calls routine FeIIPunPop */
                                strcpy( punch.chPunch[punch.npunch], "FE2P" );
                                punch.punarg[punch.npunch][0] = 0.;
                                punch.punarg[punch.npunch][1] = (realnum)FeII.nFeIILevel;
                        }

                        /* optional keyword range means read lower and upper bound, do these */
                        else if( nMatch("RANG",chCard) )
                        {
                                /* punch range of levels, calls routine FeIIPunPop */
                                strcpy( punch.chPunch[punch.npunch], "FE2P" );
                                ipFFmt = 5;
                                punch.punarg[punch.npunch][0] = 
                                        (realnum)FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL);
                                punch.punarg[punch.npunch][1] = 
                                        (realnum)FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL);
                                if( lgEOL || punch.punarg[punch.npunch][0] <0 ||
                                        punch.punarg[punch.npunch][0]> punch.punarg[punch.npunch][1] ||
                                        punch.punarg[punch.npunch][1]> FeII.nFeIILevel)
                                {
                                        fprintf( ioQQQ, "There must be two numbers on this punch FeII populations range command.\n" );
                                        fprintf( ioQQQ, "These give the lower and upper levels for the range of FeII levels.\n" );
                                        fprintf( ioQQQ, "The first must be less than the second and the second must be <= the number of FeII levels.\n" );
                                        fprintf( ioQQQ, "Sorry.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }

                        else
                        {
                                /* punch a very few selected levels, calls routine FeIIPunPop */
                                strcpy( punch.chPunch[punch.npunch], "FE2p" );
                        }
                }
                else
                {
                        fprintf( ioQQQ, "There must be a second keyword on this PUNCH FEII command.\n" );
                        fprintf( ioQQQ, "The ones I know about are COLUmn, CONTinuum, "
                                "DEPArture, LEVEls, LINE, OPTIcal DEPTh, and POPUlations.\n" );
                        fprintf( ioQQQ, "Sorry.\n" );
                        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( nMatch("CONT",chCard) && !nMatch("XSPE",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( nMatch("BINS",chCard) )
                {
                        /* continuum binning */
                        strcpy( punch.chPunch[punch.npunch], "CONB" );

                        sprintf( chHeader, 
                                "#Continuum binning, enrOrg, Energy, width of cells\n" );
                }

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

                        sprintf( chHeader, 
                                "#energy\t ConEmitLocal \n" );
                }

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

                        sprintf( chHeader, 
                                "#Energy\treflec\toutward\ttotal\tline\tcont\n" );
                }

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

                        sprintf( chHeader, 
                                "#Energy\tTransmitted\n" );

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

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

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

                else if( nMatch("GRAI",chCard) )
                {
                        /* punch grain continuum in optically thin limit */
                        strcpy( punch.chPunch[punch.npunch], "CONG" );

                        sprintf( chHeader, 
                                "#energy\tgraphite\trest\ttotal\n" );
                }

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

                        sprintf( chHeader, 
                                "#Incident Continuum, Enr\tnFn \n" );
                }

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

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

                else if( nMatch("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[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,
                          INPUT_LINE_LENGTH,&lgEOL);

                        /* this is option for smallest interaction to punch, def is 1 percent */
                        punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                punch.punarg[punch.npunch][1] = 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( nMatch("EVER", chCard ) )
                        {
                                /* punch every zone */
                                punch.lgPunchEveryZone[punch.npunch] = true;
                                punch.nPunchEveryZone[punch.npunch] = 1;
                        }
                        else
                        {
                                /* only punch last zone */
                                punch.lgPunchEveryZone[punch.npunch] = false;
                                punch.nPunchEveryZone[punch.npunch] = 1;
                        }

                        /* put the header at the top of the file */
                        sprintf( chHeader, 
                                "#cell(on C scale)\tnu\tflux\tflx*cs\tFinc\totsl\totsc\toutlin\toutcon\trate/tot\tintegral\tline\tcont\n" );
                }

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

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

                        sprintf( chHeader, 
                                "#ener\tTran Contin\ttrn coef\n" );
                }

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

                        sprintf( chHeader, 
                                "#energy\t n_nu\tnuF_nu \n" );
                }

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

                        sprintf( chHeader, 
                                "#Raw Con anu\tflux\totslin\totscon\tConRefIncid\tConEmitReflec\tConInterOut\toutlin\tConEmitOut\tline\tcont\tnLines\n" );
                }

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

                        sprintf( chHeader, 
                                "#Reflected\tcont\tline\ttotal\talbedo\tConID\n" );
                }

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

                        /* >>chng 06 apr 03, add "every" option to punch this on every zone -
                         * if every is not present then only last zone is punched */
                        if( nMatch("EVER", chCard ) )
                        {
                                /* punch every zone */
                                punch.lgPunchEveryZone[punch.npunch] = true;
                                /* option to say how many to skip */
                                ipFFmt = 5;
                                punch.nPunchEveryZone[punch.npunch] = (long)FFmtRead(chCard,&ipFFmt,
                                  INPUT_LINE_LENGTH,&lgEOL);
                                if( lgEOL )
                                        punch.nPunchEveryZone[punch.npunch] = 1;
                        }
                        else
                        {
                                /* only punch last zone */
                                punch.lgPunchEveryZone[punch.npunch] = false;
                                punch.nPunchEveryZone[punch.npunch] = 1;
                        }
                }
        }

        /* 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( nMatch("CONV",chCard) )
        {
                if( nMatch("REAS",chCard) )
                {
                        /* this does not count as a punch option (really) */
                        punch.lgPunConv = true;
                        /* this is done below */
                        strcpy( punch.chPunch[punch.npunch], "" );
                        punch.lgRealPunch[punch.npunch] = false;
                }
                else if( nMatch("ERRO",chCard) )
                {
                        /* punch zone by zone errors in pressure, electron density, and heating-cooling */
                        /* convergence error */
                        strcpy( punch.chPunch[punch.npunch], "CNVE" );
                        sprintf( chHeader, 
                                "#depth\tnPres2Ioniz\tP(cor)\tP(cur)\tP%%error\tNE(cor)\tNE(cur)\tNE%%error\tHeat\tCool\tHC%%error\n" );
                }
                else if( nMatch("BASE",chCard) )
                {
                        /* punch converged quantities in Converge base for each pass through
                         * solvers - not one pass per zone */
                        strcpy( punch.chPunch[punch.npunch], "CNVB" );
                        strcpy( punch.chPunch[punch.npunch], "" );
                        punch.lgRealPunch[punch.npunch] = false;
                }
                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" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch(" DR ",chCard) )
        {
                /* first occurrence of punch dr to follow choice in change of zoning */
                punch.lgDROn = true;
                strcpy( punch.chPunch[punch.npunch], "" );
                punch.lgRealPunch[punch.npunch] = false;
        }

        else if( nMatch("ELEM",chCard) && !nMatch("GAMMA" , chCard) )
        {
                /* option to punch ionization structure of some element
                 * will give each stage of ionization, vs depth */
                strcpy( punch.chPunch[punch.npunch], "ELEM" );

                /* this returns element number on c scale */
                /* >>chng 04 nov 23, had converted to f scale, leave on c */
                nelem = GetElem(chCard);

                /* this is the atomic number on the c scale */
                punch.punarg[punch.npunch][0] = (realnum)nelem;
                ASSERT( nelem>=ipHYDROGEN);

                /* >>chng 04 nov 24, add DENSE option to print density rather than fraction */
                punch.punarg[punch.npunch][1] = 0;
                if( nMatch("DENS",chCard)  )
                        punch.punarg[punch.npunch][1] = 1.;

                /* start printing header line - first will be the depth in cm */
                sprintf( chHeader, "#depth");

                /* next come the nelem+1 ion stages */
                for(i=0; i<=nelem+1;++i )
                {
                        sprintf( chTemp, 
                                "\t%2s%2li", elementnames.chElementSym[nelem],i+1);
                        strcat( chHeader, chTemp );
                }

                /* finally some fine structure or molecular populations */
                /* >>chng 04 nov 23, add fs pops of C, O 
                 * >>chng 04 nov 25, add molecules */
                if( nelem==ipHYDROGEN )
                {
                        sprintf( chTemp, "\tH2");
                        strcat( chHeader, chTemp );
                }
                if( nelem==ipCARBON )
                {
                        sprintf( chTemp, "\tC1\tC1*\tC1**\tC2\tC2*\tCO");
                        strcat( chHeader, chTemp );
                }
                else if( nelem==ipOXYGEN )
                {
                        sprintf( chTemp, "\tO1\tO1*\tO1**");
                        strcat( chHeader, chTemp );
                }

                /* finally the new line */
                strcat( chHeader, "\n");
        }

        else if( nMatch("FITS",chCard) )
        {

#ifdef FLT_IS_DBL
                fprintf( ioQQQ, "Punching FITS files is not currently supported in double precision.\n" );
                fprintf( ioQQQ, "Please recompile without the FLT_IS_DBL option.\n" );
                fprintf( ioQQQ, "Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
#else
                /* say that this is a FITS file output */
                punch.lgFITS[punch.npunch] = true;

                strcpy( punch.chPunch[punch.npunch], "FITS" );
#endif

        }

        else if( nMatch("FRED",chCard) )
        {
                /* punch out some stuff for Fred's dynamics project */
                sprintf( chHeader, 
                        "#Radius\tDepth\tVelocity\thden\teden\tTemperature\tRadAccel line\tRadAccel con\t"
                        "Force multiplier\t"
                        "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t"
                        "O2\tO3\tO4\tO5\tO6\tO7\tO8\t" 
                        "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t"
                        "O2\tO3\tO4\tO5\tO6\tO7\tO8\tMg2\tMg2\n");
                strcpy( punch.chPunch[punch.npunch], "FRED" );
        }

        else if( nMatch("GAMM",chCard) )
        {
                /* punch all photoionization rates for all subshells */
                sprintf( chHeader, 
                        "#Photoionization rates \n" );
                if( nMatch("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 */
                        nelem = GetElem(chCard);
                        /* this is the atomic number on the C scale */
                        punch.punarg[punch.npunch][0] = (realnum)nelem;

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

        }
        else if( nMatch("GRAI",chCard) )
        {
                /* punch grain ... options */
                if( nMatch("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 */
                        sprintf( chHeader, 
                                "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" );
                }
                else if( nMatch("ABUN",chCard) )
                {
                        /* punch grain abundance */
                        strcpy( punch.chPunch[punch.npunch], "DUSA" );
                        sprintf( chHeader, 
                                 "#grain\tdepth\tabundance (g/cm^3)\n" );
                }
                else if( nMatch("D/G ",chCard) )
                {
                        /* punch grain dust/gas mass ratio */
                        strcpy( punch.chPunch[punch.npunch], "DUSD" );
                        sprintf( chHeader, 
                                 "#grain\tdepth\tdust/gas mass ratio\n" );
                }
                else if( nMatch("PHYS",chCard) )
                {
                        /* punch grain physical conditions */
                        strcpy( punch.chPunch[punch.npunch], "DUSP" );
                        sprintf( chHeader, 
                                "#grain\tdepth\tpotential\n" );
                }
                else if( nMatch(" QS ",chCard) )
                {
                        strcpy( punch.chPunch[punch.npunch], "DUSQ" );
                        sprintf( chHeader, 
                                "#grain\tnu\tQ_abs\tQ_scat\n" );
                }
                else if( nMatch("TEMP",chCard) )
                {
                        /* punch temperatures of each grain species */
                        strcpy( punch.chPunch[punch.npunch], "DUST" );
                        /* cannot punch grain labels since they are not known yet */
                        sprintf( chHeader, 
                                "#grain temperature\n" );
                }
                else if( nMatch("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 */
                        sprintf( chHeader, 
                                "#grain drift velocity\n" );
                }
                else if( nMatch("EXTI",chCard) )
                {
                        /* punch grain extinction */
                        strcpy( punch.chPunch[punch.npunch], "DUSE" );
                        /* cannot punch grain labels since they are not known yet */
                        sprintf( chHeader, 
                                "#depth\tA_V(extended)\tA_V(point)\n" );
                }
                else if( nMatch("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 */
                        sprintf( chHeader, 
                                "#grain charge\n" );
                }
                else if( nMatch("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 */
                        sprintf( chHeader, 
                                "#grain heating\n" );
                }
                else if( nMatch("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 */
                        sprintf( chHeader, 
                                "#grain\tdepth\tpotential\n" );
                }
                else if( nMatch("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 */
                        sprintf( chHeader, 
                                "#grain H2 formation rates\n" );
                }
                else
                {
                        fprintf( ioQQQ, "There must be a second key on this GRAIN command; The options I know about follow (required key in CAPS):\n");
                        fprintf( ioQQQ, "OPACity, ABUNdance, D/G mass ratio, PHYSical conditions,  QS , TEMPerature, DRIFt velocity, EXTInction, CHARge, HEATing, POTEntial, H2RAtes\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("GAUN",chCard) )
        {
                strcpy( punch.chPunch[punch.npunch], "GAUN" );
                sprintf( chHeader, 
                        "#Gaunt factors.\n" );
        }
        else if( nMatch("GRID",chCard) )
        {
                strcpy( punch.chPunch[punch.npunch], "GRID" );
                /* automatically generate no hash option */
                punch.lgHashEndIter[punch.npunch] = false;
        }
        else if( nMatch( "HIST" , chCard ) )
        {
                /* punch pressure history of current zone */
                if( nMatch( "PRES" , chCard ) )
                {
                        /* punch pressure history - density - pressure for this zone */
                        strcpy( punch.chPunch[punch.npunch], "HISp" );
                        sprintf( chHeader, 
                                "#iter zon\tdensity\tpres cur\tpres corr\n" );
                }
                /* punch temperature history of current zone */
                else if( nMatch( "TEMP" , chCard ) )
                {
                        /* punch pressure history - density - pressure for this zone */
                        strcpy( punch.chPunch[punch.npunch], "HISt" );
                        sprintf( chHeader, 
                                "#iter zon\ttemperature\theating\tcooling\n" );
                }
        }

        else if( nMatch("HTWO",chCard) )
        {
                fprintf(ioQQQ," Sorry, this command has been replaced with the "
                        "PUNCH H2 CREATION and PUNCH H2 DESTRUCTION commands.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* QHEAT has to come before HEAT... */
        else if( nMatch("QHEA",chCard) )
        {
                /* this is just a dummy clause, do the work below after parsing is over. 
                 * this is a no-nothing, picked up to stop optimizer */
                ((void)0);
        }

        else if( nMatch("HEAT",chCard) )
        {
                /* punch heating */
                strcpy( punch.chPunch[punch.npunch], "HEAT" );
                /*>>chng 06 jun 06, revise to be same as punch cooling */
                sprintf( chHeader, 
                        "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\theat fracs\n" );
        }

        else if( nMatch("HELI",chCard) &&!( nMatch("IONI",chCard)))
        {
                /* punch helium & helium-like iso sequence, but not punch helium ionization rate
                 * punch helium line wavelengths */
                if( nMatch("LINE",chCard) && nMatch("WAVE",chCard) )
                {
                        strcpy( punch.chPunch[punch.npunch], "HELW" );
                        sprintf( chHeader, 
                                "#wavelengths of lines from He-like ions\n" );
                }
                else
                {
                        fprintf( ioQQQ, "punch helium has options: LINE WAVElength.\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                        /* no key */
                }
        }

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

        else if( nMatch("HYDR",chCard) )
        {
                /* punch hydrogen physical conditions */
                if( nMatch("COND",chCard) )
                {
                        strcpy( punch.chPunch[punch.npunch], "HYDc" );
                        sprintf( chHeader, 
                                "#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 - accept either keyword 21cm or 21 cm */
                else if( nMatch("21 CM",chCard) ||nMatch("21CM",chCard))
                {
                        /* punch information about 21 cm line */
                        strcpy( punch.chPunch[punch.npunch], "21CM" );
                        sprintf( chHeader, 
                                "#depth\tT(spin)\tT(kin)\tT(Lya/21cm)\tnLo\tnHi\tOccLya\ttau(21cm)"
                                "\ttau(Lya)\topac(21 cm)\tn/Ts\ttau(21)\tTex(Lya)\tN(H0)/Tspin"
                                "\tSum_F0\tSum_F1\tSum_T21\n" );
                }

                else if( nMatch("IONI",chCard) )
                {
                        /* punch hydrogen ionization */
                        strcpy( punch.chPunch[punch.npunch], "HYDi" );
                        sprintf( chHeader, 
                                "#hion\tzn\tgam1\tcoll ion1\tRecTot\tHRecCaB\thii/hi\tSim hii/hi"
                                "\time_Hrecom_long(esc)\tdec2grd\texc pht\texc col\trec eff\tsec ion\n" );
                }
                else if( nMatch("POPU",chCard) )
                {
                        /* punch hydrogen populations */
                        strcpy( punch.chPunch[punch.npunch], "HYDp" );
                        sprintf( chHeader, 
                                "#depth\tn(H0)\tn(H+)\tn(1s)\tn(2s)\tn(2p)\tetc\n" );
                }
                else if( nMatch("LINE",chCard) )
                {
                        /* punch hydrogen lines
                         * hydrogen line intensities and optical depths  */
                        strcpy( punch.chPunch[punch.npunch], "HYDl" );
                        sprintf( chHeader, 
                                "#nup\tnlo\tE(ryd)\ttau\n" );
                }
                else if( nMatch(" LYA",chCard) )
                {
                        /* punch hydrogen Lya some details about Lyman alpha  */
                        strcpy( punch.chPunch[punch.npunch], "HYDL" );
                        sprintf( chHeader, 
                                "#depth\tTauIn\tTauTot\tn(2p)/n(1s)\tTexc\tTe\tTex/T\tPesc\tPdes\tpump\topacity\talbedo\n" );
                }
                else
                {
                        fprintf( ioQQQ, "Punch hydrogen has options: CONDitions, 21 CM, LINE, POPUlations, and IONIzation.\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("IONI",chCard) )
        {
                if( nMatch("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" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        punch.punarg[punch.npunch][0] = (realnum)nelem;
                        strcpy( punch.chPunch[punch.npunch], "IONR" );
                        sprintf( chHeader, 
                                "#%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" );
                        sprintf( chHeader, 
                                "#Mean ionization distribution\n" );
                }
        }

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

        else if( nMatch("LEID",chCard) )
        {
                if( nMatch( "LINE" , chCard ) )
                {
                        /* punch Leiden lines
                         * final intensities of the Leiden PDR models */
                        strcpy( punch.chPunch[punch.npunch], "LEIL" );
                        sprintf( chHeader, "#ion\twl\tInt\trel int\n");
                }
                else
                {
                        /* punch Leiden structure
                        * structure of the Leiden PDR models */
                        strcpy( punch.chPunch[punch.npunch], "LEIS" );
                        sprintf( chHeader, 
                                /* 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( nMatch("LINE",chCard) && nMatch("LIST",chCard) )
        {
                /* punch line list "output file" "Line List file" */
                long int j;
                /* 
                 * we parsed off the second file name at start of this routine
                 * check if file was found, use it if it was, else abort
                 */
                if( !lgSecondFilename )
                {
                        fprintf(ioQQQ , "There must be a second file name between "
                                "double quotes on the PUNCH LINE LIST command.  This second"
                                " file contains the input line list.  I did not find it.\nSorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }

                /* actually get the lines, and malloc the space in the arrays 
                 * cdGetLineList will look on path */
                if( punch.ipPnunit[punch.npunch] == NULL  &&
                        ( punch.nLineList[punch.npunch] = cdGetLineList( chSecondFilename , 
                        &punch.chLineListLabel[punch.npunch] , 
                        &punch.wlLineList[punch.npunch] ) ) < 0 )
                {
                        fprintf(ioQQQ,"DISASTER could not open PUNCH LINE LIST file %s \n",
                                chSecondFilename );
                        cdEXIT(EXIT_FAILURE);
                }

                /* ratio option, in which pairs of lines form ratios, first over
                 * second */
                if( nMatch("RATI",chCard) )
                {
                        punch.lgLineListRatio[punch.npunch] = true;
                        if( punch.nLineList[punch.npunch]%2 )
                        {
                                /* odd number of lines - cannot take ratio */
                                fprintf(ioQQQ , "There must be an even number of lines to"
                                        " take ratios of lines.  There were %li, an odd number."
                                        "\nSorry.\n", punch.nLineList[punch.npunch]);
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        /* no ratio */
                        punch.lgLineListRatio[punch.npunch] = false;
                }

                /* do special string saying our job, and then print labels */
                strcpy( punch.chPunch[punch.npunch], "LLST" );

                /* keyword absolute says to do absolute rather than relative intensities 
                 * relative intensities are the default */
                if( nMatch("ABSO",chCard) )
                {
                        punch.punarg[punch.npunch][0] = 1;
                }
                else
                {
                        punch.punarg[punch.npunch][0] = 0;
                }

                /* give header line */
                sprintf( chHeader, "#lineslist" );
                for( j=0; j<punch.nLineList[punch.npunch]; ++j )
                {
                        /* if taking ratio then put div sign between pairs */
                        if( punch.lgLineListRatio[punch.npunch] && is_odd(j) )
                                strcat( chHeader , "/" );
                        else
                                strcat( chHeader , "\t" );
                        sprintf( chTemp, "%s ", punch.chLineListLabel[punch.npunch][j] ); 
                        strcat( chHeader, chTemp );
                        sprt_wl( chTemp, punch.wlLineList[punch.npunch][j] );
                        strcat( chHeader, chTemp );
                }
                strcat( chHeader, "\n" );
        }

        else if( nMatch("LINE",chCard) && !nMatch("XSPE",chCard) && !nMatch("NEAR",chCard))
        {
                /* punch line emissivity -
                 * this is not punch xspec lines and not linear option
                 * check for keyword UNITS on line, then scan wavelength or energy units, 
                 * sets punch.chConPunEnr*/
                ChkUnits(chCard);

                /* punch line emissivity, line intensity, line array,
                 * and line data */
                if( nMatch("STRU",chCard) )
                {
                        fprintf(ioQQQ," The     PUNCH LINES STRUCTURE command is now PUNCH LINES "
                                "EMISSIVITY.\n Sorry.\n\n");
                        cdEXIT(EXIT_FAILURE);
                }

                else if( nMatch("EMIS",chCard) )
                {
                        /* this used to be the punch lines structure command, is now
                         * the punch lines emissivity command 
                         * give line emissivity vs depth */
                        strcpy( punch.chPunch[punch.npunch], "LINS" );
                        sprintf( chHeader, 
                                "#Emission line structure:");
                        /* read in the list of lines to examine */
                        punch_line(punch.ipPnunit[punch.npunch],"READ",false, chTemp);
                        strcat( chHeader, chTemp );
                }

                else if( nMatch(" RT ", chCard ) )
                {
                        /* punch line RT */
                        strcpy( punch.chPunch[punch.npunch], "LINR" );
                        /* punch some details needed for line radiative transfer 
                         * routine in punch_line.c */
                        Punch_Line_RT(punch.ipPnunit[punch.npunch],"READ");
                }

                else if( nMatch("CUMU",chCard) )
                {
                        /* punch lines cumulative 
                         * this will be integrated line intensity, function of depth */
                        strcpy( punch.chPunch[punch.npunch], "LINC" );
                        /* option for either relative intensity or abs luminosity */
                        if( nMatch("RELA",chCard) )
                        {
                                lgEOL = true;
                                sprintf( chHeader, 
                                        "#Cumulative emission line relative intensity.\n" );
                        }
                        else
                        {
                                sprintf( chHeader, 
                                        "#Cumulative emission line absolute intensity.\n" );
                                lgEOL = false;
                        }
                        /* read in the list of lines to examine */
                        punch_line(punch.ipPnunit[punch.npunch],"READ",lgEOL,chTemp);
                        strcat( chHeader, chTemp );
                }

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

                else if( nMatch("ARRA",chCard) )
                {
                        /* punch line array -
                         * output energies and luminosities of predicted lines */
                        strcpy( punch.chPunch[punch.npunch], "LINA" );
                        sprintf( chHeader, 
                                "#enr\tID\tI(intrinsic)\tI(emergent)\ttype\n" );
                }

                else if( nMatch("LABE",chCard) )
                {
                        /* punch line labels */
                        strcpy( punch.chPunch[punch.npunch], "LINL" );
                        sprintf( chHeader, 
                                "#index\tlabel\twavelength\tcomment\n" );
                        /* this controls whether we will print lots of redundant 
                         * info labels for transferred lines - if keyword LONG appears
                         * then do so, if does not appear then do not - this is default */
                        if( nMatch("LONG",chCard) )
                                punch.punarg[punch.npunch][0] = 1;
                        else
                                punch.punarg[punch.npunch][0] = 0;
                }

                else if( nMatch("OPTI",chCard) )
                {
                        /* punch line optical depths, done in PunchLineStuff, located in punchdo.c */
                        strcpy( punch.chPunch[punch.npunch], "LINO" );
                        sprintf( chHeader, 
                                "#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( nMatch("UNIT",chCard) )
                                ChkUnits(chCard);

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

                else if( nMatch("POPU",chCard) )
                {
                        /* punch line populations command - first give index and inforamtion
                         * for all lines, then populations for lines as a function of
                         * depth, using this index */
                        strcpy( punch.chPunch[punch.npunch], "LINP" );
                        sprintf( chHeader, 
                                "#population information\n" );
                        /* this is optional limit to smallest population to punch - always
                         * interpreted as a log */
                        ipFFmt = 5;
                        punch.punarg[punch.npunch][0] = (realnum)pow(10.,FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL));

                        /* this is default - all positive populations */
                        if( lgEOL )
                                punch.punarg[punch.npunch][0] = 0.f;

                        if( nMatch(" OFF",chCard) )
                        {
                                /* no lower limit - print all lines */
                                punch.punarg[punch.npunch][0] = -1.f;
                        }
                }

                else if( nMatch("INTE",chCard) )
                {
                        /* this will be full set of line intensities */
                        strcpy( punch.chPunch[punch.npunch], "LINI" );
                        sprintf( chHeader, 
                                "#Emission line intensities per unit inner area\n" );
                        if( nMatch("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( nMatch("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" );
                                        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" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch(" MAP",chCard) )
        {
                strcpy( punch.chPunch[punch.npunch], "MAP " );
                sprintf( chHeader, 
                        "#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;
                hcmap.MapZone = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        hcmap.MapZone = 1;
                }

                if( nMatch("RANG",chCard) )
                {
                        bool lgLogOn;
                        hcmap.RangeMap[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( hcmap.RangeMap[0] <= 10. && !nMatch("LINE",chCard) )
                        {
                                hcmap.RangeMap[0] = (realnum)pow((realnum)10.f,hcmap.RangeMap[0]);
                                lgLogOn = true;
                        }
                        else
                        {
                                lgLogOn = false;
                        }

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

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

        else if( nMatch("MOLE",chCard) )
        {
                /* molecules, especially for PDR calculations */
                strcpy( punch.chPunch[punch.npunch], "MOLE" );
                sprintf( chHeader, 
                        "#molecular species will follow:\n");
        }

        else if( nMatch("OPTICAL",chCard) && nMatch("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( nMatch("FINE",chCard) )
                {
                        /* punch fine continuum optical depths */
                        rfield.lgPunchOpacityFine = true;
                        strcpy( punch.chPunch[punch.npunch], "OPTf" );
                        sprintf( chHeader, "#energy\tTau tot\topacity\n" );
                        ipFFmt = 5;
                        /* range option - important since so much data */
                        if( nMatch("RANGE",chCard) ) 
                        {
                                /* get lower and upper range, must be in Ryd */
                                punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
                                punch.punarg[punch.npunch][1] = (realnum)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");
                                        cdEXIT(EXIT_FAILURE);
                                }
                                if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] )
                                {
                                        fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else
                        {
                                /* these mean full energy range */
                                punch.punarg[punch.npunch][0] = 0.;
                                punch.punarg[punch.npunch][1] = 0.;
                        }
                        /* optional last parameter - how many points to bring together */
                        punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt,
                          INPUT_LINE_LENGTH,&lgEOL);
                        /* default is to bring together ten */
                        if( lgEOL )
                                punch.punarg[punch.npunch][2] = 10;
                        if( punch.punarg[punch.npunch][2] < 1 )
                        {
                                fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        /* punch coarse continuum optical depths */
                        strcpy( punch.chPunch[punch.npunch], "OPTc" );
                        sprintf( chHeader, 
                                "#energy\ttotal\tabsorp\tscat\n" );
                }

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

        else if( nMatch("OVER",chCard) && nMatch(" OVE",chCard) )
        {
                /* punch overview of model results */
                strcpy( punch.chPunch[punch.npunch], "OVER" );
                sprintf( chHeader, 
                        "#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( nMatch(" PDR",chCard) )
        {
                strcpy( punch.chPunch[punch.npunch], " PDR" );
                sprintf( chHeader, 
                        "#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( nMatch("PHYS",chCard) )
        {
                /* punch physical conditions */
                strcpy( punch.chPunch[punch.npunch], "PHYS" );
                sprintf( chHeader, 
                        "#PhyC depth\tTe\tn(H)\tn(e)\tHtot\taccel\tfillfac\n" );
        }

        else if( nMatch("POIN",chCard) )
        {
                /* punch out the pointers */
                punch.lgPunPoint = true;
                /* this does not count as a punch option (really) */
                strcpy( punch.chPunch[punch.npunch], "" );
                punch.lgRealPunch[punch.npunch] = false;
        }

        else if( nMatch("PRES",chCard) )
        {
                /* the punch pressure command */
                strcpy( punch.chPunch[punch.npunch], "PRES" );
                sprintf( chHeader, 
                        "#P depth\tPcorrect\tPcurrent\tPIn+Pinteg\tPgas(r0)\tPgas\tPram\tPrad(line)\tPinteg\tV(wind km/s)\tcad(wind km/s)\tP(mag)\tV(turb km/s)\tP(turb)\tconv?\n" );
        }

        else if( nMatch("RADI",chCard) )
        {
                /* the punch radius command */
                sprintf( chHeader, "#NZONE\tradius\tdepth\tdr\n" );
                /* option to only punch the outer radius */
                if( nMatch( "OUTE" , chCard ) )
                {
                        /* only outer radius */
                        strcpy( punch.chPunch[punch.npunch], "RADO" );
                }
                else
                {
                        /* all radii */
                        strcpy( punch.chPunch[punch.npunch], "RADI" );
                }
        }

        else if( nMatch("RECO",chCard) )
        {
                if( nMatch("COEF",chCard) )
                {
                        /* recombination coefficients for everything */

                        /* this is logical flag used in routine ion_recom to create the punch output */
                        punch.lgioRecom = true;
                        /* this does not count as a punch option (really) */
                        strcpy( punch.chPunch[punch.npunch], "" );
                        punch.lgRealPunch[punch.npunch] = false;
                }

                else if( nMatch("EFFI",chCard) )
                {
                        /* punch recombination efficiency */
                        strcpy( punch.chPunch[punch.npunch], "RECE" );
                        sprintf( chHeader, 
                                "#Recom 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" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* punch results command, either as single column or wide array */
        else if( nMatch("RESU",chCard) )
        {
                strcpy( punch.chPunch[punch.npunch], "RESU" );
                if( nMatch("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 */
                sprintf( chHeader, 
                        "#results of calculation\n" );
        }

        else if( nMatch("SECO",chCard) )
        {
                /* punch secondary ionization rate */
                strcpy( punch.chPunch[punch.npunch], "SECO" );
                sprintf( chHeader, 
                        "#depth\tIon(H^0)\tDiss(H_2)\tExcit(Lya)\n" );
        }

        else if( nMatch("SOUR",chCard) )
        {
                if( nMatch("DEPT",chCard) )
                {
                        /* print continuum source function as function of depth */
                        strcpy( punch.chPunch[punch.npunch], "SOUD" );
                        sprintf( chHeader, 
                                "#continuum source function vs depth\n" );
                }
                else if( nMatch("SPEC",chCard) )
                {
                        /* print spectrum continuum source function at 1 depth */
                        strcpy( punch.chPunch[punch.npunch], "SOUS" );
                        sprintf( chHeader, 
                                "#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" );
                        cdEXIT(EXIT_FAILURE);
                }
        }


        /* punch spectrum the new form of the punch continuum, will eventually replace the standard
         * punch continuum command */
        else if( nMatch("SPEC",chCard) && nMatch("TRUM",chCard) && !nMatch("XSPE",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" );

                sprintf( chHeader, 
                        "#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 spectrum commands */
                punch.punarg[punch.npunch][0] = (realnum)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 );
                        cdEXIT(EXIT_FAILURE);
                }

        }

        else if( nMatch("SPEC",chCard) && nMatch("CIAL",chCard) )
        {
                /* punch special, will call routine PunchSpecial */
                strcpy( punch.chPunch[punch.npunch], "SPEC" );
                sprintf( chHeader, "#Special.\n" );
        }

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

        else if( nMatch("TEMP",chCard) )
        {
                /* punch temperature command */
                strcpy( punch.chPunch[punch.npunch], "TEMP" );
                sprintf( chHeader, 
                        "#depth\tTe\tcC/dT\tdt/dr\td^2T/dr^2\n" );
        }

        else if( nMatch("TIME",chCard) && nMatch("DEPE",chCard) )
        {
                /* information about time dependent solutions */
                strcpy( punch.chPunch[punch.npunch], "TIMD" );
                /* do not want to separate iterations with special character */
                punch.lg_separate_iterations[punch.npunch] = false;
                /* write header */
                sprintf( chHeader , 
                        "#elapsed time\ttime step \tscale cont\tn(H)\t<T>\t<H+/H rad>\t<H0/H rad>\t<H2/H rad>\t<He+/H rad>\t<CO/H>\n" );
        }

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

        else if( nMatch("WIND",chCard) )
        {
                strcpy( punch.chPunch[punch.npunch], "WIND" );
                sprintf( chHeader, 
                        "#radius\tdepth\tvel [cm/s]\tTot accel [cm s-2]\tLin accel [cm s-2]"
                        "\tCon accel [cm s-2]\tforce multiplier\ta_gravity\n" );
                if( nMatch( "TERM" , chCard ) )
                {
                        /* only punch for last zone, the terminal velocity, for grids */
                        punch.punarg[punch.npunch][0] = 0.;
                }
                else
                {
                        /* one means punch every zone */
                        punch.punarg[punch.npunch][0] = 1.;
                }
        }

        else if( nMatch("XSPE",chCard) )
        {
                /* say that this is a FITS file output */
                punch.lgFITS[punch.npunch] = true;

                /* the punch xspec commands */
                if( !punch.lgPunLstIter[punch.npunch] )
                {
                        punch.lgPunLstIter[punch.npunch] = true;
                }

                if( nMatch("ATAB",chCard) )
                {
                        /* punch xspec atable command */
                        if( nMatch("INCI",chCard) )
                        {
                                if( nMatch("ATTE",chCard) )
                                {
                                        /* attenuated incident continuum */
                                        strcpy( punch.chPunch[punch.npunch], "XATT" );
                                        grid.lgOutputTypeOn[2] = true;
                                }
                                else if( nMatch("REFL",chCard) )
                                {
                                        /* reflected incident continuum */
                                        strcpy( punch.chPunch[punch.npunch], "XRFI" );
                                        grid.lgOutputTypeOn[3] = true;
                                }
                                else
                                {
                                        /* incident continuum */
                                        strcpy( punch.chPunch[punch.npunch], "XINC" );
                                        grid.lgOutputTypeOn[1] = true;
                                }
                        }
                        else if( nMatch("DIFF",chCard) )
                        {
                                if( nMatch("REFL",chCard) )
                                {
                                        /* reflected diffuse continuous emission */
                                        strcpy( punch.chPunch[punch.npunch], "XDFR" );
                                        grid.lgOutputTypeOn[5] = true;
                                }
                                else
                                {
                                        /* diffuse continuous emission outward */
                                        strcpy( punch.chPunch[punch.npunch], "XDFO" );
                                        grid.lgOutputTypeOn[4] = true;
                                }
                        }
                        else if( nMatch("LINE",chCard) )
                        {
                                if( nMatch("REFL",chCard) )
                                {
                                        /* reflected lines */
                                        strcpy( punch.chPunch[punch.npunch], "XLNR" );
                                        grid.lgOutputTypeOn[7] = true;
                                }
                                else
                                {
                                        /* outward lines */
                                        strcpy( punch.chPunch[punch.npunch], "XLNO" );
                                        grid.lgOutputTypeOn[6] = true;
                                }
                        }
                        else if( nMatch("SPEC",chCard) )
                        {
                                if( nMatch("REFL",chCard) )
                                {
                                        /* reflected spectrum */
                                        strcpy( punch.chPunch[punch.npunch], "XREF" );
                                        grid.lgOutputTypeOn[9] = true;
                                }
                                else
                                {
                                        /* transmitted spectrum */
                                        strcpy( punch.chPunch[punch.npunch], "XTRN" );
                                        grid.lgOutputTypeOn[8] = true;
                                }
                        }
                        else
                        {
                                /* transmitted spectrum */
                                strcpy( punch.chPunch[punch.npunch], "XTRN" );
                                grid.lgOutputTypeOn[8] = true;
                        }
                }
                else if( nMatch("MTAB",chCard) )
                {
                        /* punch xspec mtable */
                        strcpy( punch.chPunch[punch.npunch], "XSPM" );
                        grid.lgOutputTypeOn[10] = true;
                }
                else
                {
                        fprintf( ioQQQ, "Support only for xspec atable and xspec mtable.\n" );
                        cdEXIT( EXIT_FAILURE );
                }
        }

        /* punch column density has to come last so do not trigger specific column 
         * densities, H2, FeII, etc.
         * Need both keywords since column is also the keyword for one line per line */
        else if( nMatch("COLU",chCard) && nMatch("DENS",chCard) )
        {
                if( nMatch("SOME" , chCard ))
                {
                        /* flag saying punch some column densities */
                        strcpy( punch.chPunch[punch.npunch], "COLS" );
                        punch_colden( chHeader , punch.ipPnunit[punch.npunch] , "READ" );
                }
                else
                {
                        /* punch them all as in table */
                        /* punch column density */
                        strcpy( punch.chPunch[punch.npunch], "COLU" );
                        sprintf( chHeader, 
                                "#column densities\n" );
                }
        }
        else
        {
                fprintf( ioQQQ, 
                        "ParsePunch cannot find a recognized keyword on this PUNCH command line.\nSorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* only open if file has not already been opened during a previous call */
        if( punch.ipPnunit[punch.npunch] == NULL )
        {
                if( nMatch("XSPE",chCard) || nMatch("FITS",chCard) )
                {
                        /* open the file with the name generated above */
                        punch.ipPnunit[punch.npunch] = open_data( chFilename, "wb", AS_LOCAL_ONLY );
                }
                else if( punch.lgPunchToSeparateFiles[punch.npunch] && optimize.lgOptimr )
                {
                        char chIndex[4];

                        /* take max because optimize.nOptimiz has not yet been incremented on first pass. */
                        sprintf( chIndex, "%li", MAX2( 1, optimize.nOptimiz ) );

                        strcat( chFilename, chIndex );

                        fixit(); 
                        /* must create unique name here.  considerations include:
                         *   1.  must zero fill to the left 
                         *   2.  must make sure that we have enough digits  (don't hard-wire 3 or 4 digits, 
                                         count the number of terms.             
                         *   3.  should we look for a file ending (i.e., ".txt" and insert number before that? */

                        /* open the file with the name generated above */
                        punch.ipPnunit[punch.npunch] = open_data( chFilename, "w", AS_LOCAL_ONLY );
                }
                else
                {
                        /* open the file with the name generated above */
                        punch.ipPnunit[punch.npunch] = open_data( chFilename, "w", AS_LOCAL_ONLY );
                }

                /* option to set no buffering for this file.  The setbuf command may 
                 * ONLY be issued right after the open of the file.  Giving it after 
                 * i/o has been done may result in loss of the contents of the buffer, PvH */
                if( nMatch("NO BUFFER",chCard) )
                        setbuf( punch.ipPnunit[punch.npunch] , NULL );
        }

        /***************************************************************/
        /*                                                             */
        /*  The following are special punch options and must be done   */
        /*  after the parsing and file opening above.                  */
        /*                                                             */
        /*  NB:  these are ALSO parsed above.  Here we DO something.   */
        /*                                                             */
        /***************************************************************/

        if( nMatch("CONV",chCard) && nMatch("REAS",chCard) )
        {
                /* punch reason model declared not converged
                 * not a true punch command, since done elsewhere */
                punch.ipPunConv = punch.ipPnunit[punch.npunch];
                lgPunConv_noclobber = lgNoClobber[punch.npunch];
                punch.lgPunConv = true;
                fprintf( punch.ipPunConv, 
                        "# reason for continued iterations\n" );
                strcpy( punch.chPunch[punch.npunch], "" );
                punch.lgRealPunch[punch.npunch] = false;
        }

        else if( nMatch("CONV",chCard) && nMatch("BASE",chCard) )
        {
                /* punch some quantities we are converging */
                punch.lgTraceConvergeBase = true;
                /* the second punch occurrence - file has been opened,
                * copy handle, also pass on special no hash option */
                if( nMatch("O HA",chCard) )
                        punch.lgTraceConvergeBaseHash = false;
                punch.ipTraceConvergeBase = punch.ipPnunit[punch.npunch];
                /* set punch last flag to whatever it was above */
                lgTraceConvergeBase_noclobber = lgNoClobber[punch.npunch];
                static bool lgPrtHeader = true;
                if( lgPrtHeader )
                        fprintf( punch.ipTraceConvergeBase, 
                        "#zone\theat\tcool\teden\n" );
                lgPrtHeader = false;
        }

        else if( nMatch(" DR ",chCard) )
        {
                static bool lgPrtHeader = true;
                /* the second punch dr occurrence - file has been opened,
                 * copy handle to ipDRout, also pass on special no hash option */
                if( nMatch("O HA",chCard) )
                        punch.lgDRHash = false;
                punch.ipDRout = punch.ipPnunit[punch.npunch];
                /* set punch last flag to whatever it was above */
                punch.lgDRPLst = punch.lgPunLstIter[punch.npunch];
                lgDROn_noclobber = lgNoClobber[punch.npunch];
                if( lgPrtHeader )
                        fprintf( punch.ipDRout, 
                        "#zone\tdepth\tdr\tdr 2 go\treason \n" );
                lgPrtHeader = false;
                strcpy( punch.chPunch[punch.npunch], "" );
                punch.lgRealPunch[punch.npunch] = false;
        }

        else if( nMatch("QHEA",chCard) )
        {
                gv.QHPunchFile = punch.ipPnunit[punch.npunch];
                gv.lgQHPunLast = punch.lgPunLstIter[punch.npunch];
                lgQHPunchFile_noclobber = lgNoClobber[punch.npunch];
                fprintf( gv.QHPunchFile, 
                        "#Probability distributions from quantum heating routine.\n" );
        }

        else if( nMatch("POIN",chCard) )
        {
                /* punch out the pointers */
                punch.ipPoint = punch.ipPnunit[punch.npunch];
                lgPunPoint_noclobber = lgNoClobber[punch.npunch];
                punch.lgPunPoint = true;
                fprintf( punch.ipPoint, 
                        "#pointers. \n" );
                strcpy( punch.chPunch[punch.npunch], "" );
                punch.lgRealPunch[punch.npunch] = false;
        }

        else if( nMatch("RECO",chCard) && nMatch("COEF",chCard) )
        {
                /* recombination coefficients for everything
                 * punch.lgioRecom set to false 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 
                 * diel supres, helium, hydrorecom, iibod, and makerecomb*/
                punch.ioRecom = punch.ipPnunit[punch.npunch];
                lgioRecom_noclobber = lgNoClobber[punch.npunch];
                /* this is logical flag used in routine ion_recom to create the punch output */
                punch.lgioRecom = true;
                fprintf( punch.ioRecom, 
                        "#recombination coefficients cm3 s-1 for current density and temperature\n" );
                strcpy( punch.chPunch[punch.npunch], "" );
                punch.lgRealPunch[punch.npunch] = false;
        }

        else if( nMatch(" MAP",chCard) )
        {
                /* say output goes to special punch */
                ioMAP = punch.ipPnunit[punch.npunch];
        }

        /* check that string written into chHeader can actually fit there
         * we may have overrun this buffer, an internal error */
        /* check that there are less than nChar characters in the line */
        char *chEOL = strchr(chHeader , '\0' );

        /* return null if input string longer than nChar, the longest we can read.
        * Print and return null but chLine still has as much of the line as 
        * could be placed in cdLine */
        if( (chEOL==NULL) || (chEOL - chHeader)>=MAX_HEADER_SIZE-1 )
        {
                fprintf( ioQQQ, "DISASTER chHeader has been overwritten "
                        "with a line too long to be read.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* if flag set and cdHeader not still equal to initialized string, print header
         * punch.lgPunHeader normally true, is false during grid calculation */
        if( nSimThisCoreload > 1 )
                punch.lgPunHeader = false;
        if( punch.lgPunHeader && !nMatch("ArNdY38dZ9us4N4e12SEcuQ",chHeader) )
        {
                fprintf( punch.ipPnunit[punch.npunch], "%s", chHeader );
        }

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

/*PunchFilesInit initialize punch file pointers, called from InitCoreload
 * called one time per core load 
 * NB KEEP THIS ROUTINE SYNCHED UP WITH THE NEXT ONE, ClosePunchFiles */
void PunchFilesInit(void)
{
        long int i;
        static bool lgFIRST = true;

        DEBUG_ENTRY( "PunchFilesInit()" );

        ASSERT( lgFIRST );
        lgFIRST = false;

        /* set lgNoClobber to not overwrite files, reset with clobber on punch line
         * if we are running a grid (grid command entered in cdRead) grid.lgGrid 
         * true, is false if single sim.  For grid we want to not clobber files 
         * by default, do clobber for optimizer since this was behavior before */
        bool lgNoClobberDefault = false;
        if( grid.lgGrid )
        {
                /* cdRead encountered grid command - do not want to clobber files */
                lgNoClobberDefault = true;
        }

        for( i=0; i < LIMPUN; i++ )
        {
                lgNoClobber[i] = lgNoClobberDefault;
        }
        lgPunConv_noclobber = lgNoClobberDefault;
        lgDROn_noclobber = lgNoClobberDefault;
        lgTraceConvergeBase_noclobber = lgNoClobberDefault;
        lgPunPoint_noclobber = lgNoClobberDefault;
        lgioRecom_noclobber = lgNoClobberDefault;
        lgQHPunchFile_noclobber = lgNoClobberDefault;

        /* done parsing, now turn flag to punch headers, the default behavior
         * this is set false during a grid calculation to do only one header per grid*/
        punch.lgPunHeader = true;

        for( i=0; i < LIMPUN; i++ )
        {
                if( !lgNoClobber[i] )
                {
                        punch.ipPnunit[i] = NULL;
                }
                /* set false with the dummy punch commands like punch dr */
                punch.lgRealPunch[i] = true;
                /* 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 punch resolution command */
                punch.cp_resolving_power[i] = 0.;
        }

        punch.lgTraceConvergeBase = false;

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

        if( !lgTraceConvergeBase_noclobber )
        {
                punch.ipTraceConvergeBase = NULL;
                punch.lgTraceConvergeBase = false;
        }

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

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

        if( !lgQHPunchFile_noclobber )
                gv.QHPunchFile = NULL;

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

        ioMAP = NULL;
        return;
}

/*ClosePunchFiles close punch files called from cdEXIT upon termination,
 * from cloudy before returning 
 * NB - KEEP THIS ROUTINE SYNCHED UP WITH THE PREVIOUS ONE, PunchFilesInit */
void ClosePunchFiles( bool lgFinal )
{
        long int i;

        DEBUG_ENTRY( "ClosePunchFiles()" );

        /* close all punch units cloudy opened with punch command,
         * lgNoClobber is set false with CLOBBER option on punch, says to
         * overwrite the files */
        for( i=0; i < punch.npunch; i++ )
        {
                /* if lgFinal is true, we close everything, no matter what. 
                 * this means ignoring "no clobber" options */
                if( punch.ipPnunit[i] != NULL && ( !lgNoClobber[i] || lgFinal ) )
                {
                        /* Test that any FITS files are the right size! */ 
                        if( punch.lgFITS[i] )
                        {
                                /* \todo 2 This overflows for file sizes larger (in bytes) than
                                 * a long int can represent (about 2GB on most 2007 systems)  */
                                fseek(punch.ipPnunit[i], 0, SEEK_END);
                                long file_size = ftell(punch.ipPnunit[i]);
                                if( file_size%2880 )
                                {
                                        fprintf( ioQQQ, " PROBLEM  FITS file is wrong size!\n" );
                                }
                        }

                        fclose( punch.ipPnunit[i] );
                        punch.ipPnunit[i] = NULL;
                }
        }

        /* following file handles are aliased to ipPnunit which
         * was closed above */
        if( punch.ipDRout != NULL && !lgDROn_noclobber  )
        {
                /*fclose( punch.ipDRout );*/
                punch.ipDRout = NULL;
                punch.lgDROn = false;
        }

        if( punch.ipTraceConvergeBase != NULL && !lgTraceConvergeBase_noclobber  )
        {
                /*fclose( punch.ipDRout );*/
                punch.ipTraceConvergeBase = NULL;
                punch.lgTraceConvergeBase = false;
        }

        if( punch.ipPunConv != NULL && !lgPunConv_noclobber )
        {
                /*fclose( punch.ipPunConv );*/
                punch.ipPunConv = NULL;
                punch.lgPunConv = false;
        }
        if( punch.ipPoint != NULL && !lgPunPoint_noclobber  )
        {
                /*fclose( punch.ipPoint );*/
                punch.ipPoint = NULL;
                punch.lgPunPoint = false;
        }
        if( gv.QHPunchFile != NULL  && !lgQHPunchFile_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;
        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)
{

        DEBUG_ENTRY( "ChkUnits()" );

        /* option to set units for continuum energy in punch output */
        if( nMatch("NITS",chCard) )
        {
                if( nMatch("MICR",chCard) )
                {
                        /* microns */
                        strcpy( punch.chConPunEnr[punch.npunch], "micr" );
                }
                else if( nMatch(" KEV",chCard) )
                {
                        /* keV */
                        strcpy( punch.chConPunEnr[punch.npunch], " kev" );
                }
                else if( nMatch("CENT",chCard) || nMatch(" CM ",chCard) )
                {
                        /* centimeters*/
                        strcpy( punch.chConPunEnr[punch.npunch], "cent" );
                }
                else if( nMatch(" EV ",chCard) )
                {
                        /* eV */
                        strcpy( punch.chConPunEnr[punch.npunch], " ev " );
                }
                else if( nMatch("ANGS",chCard) )
                {
                        /* angstroms */
                        strcpy( punch.chConPunEnr[punch.npunch], "angs" );
                }
                else if( nMatch("WAVE",chCard) )
                {
                        /* wavenumbers */
                        strcpy( punch.chConPunEnr[punch.npunch], "wave" );
                }
                else if( nMatch(" MHZ",chCard) )
                {
                        /* megaherz */
                        strcpy( punch.chConPunEnr[punch.npunch], " mhz" );
                }
                else if( nMatch(" 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" );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else
        {
                strcpy( punch.chConPunEnr[punch.npunch], "ryd " );
        }
        return;
}

---