/home66/gary/public_html/cloudy/c13_branch/source/parse_save.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ParseSave parse the save command */
/*SaveFilesInit initialize save file pointers, called from cdInit */
/*CloseSaveFiles close all save 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 "save.h"
#include "mpi_utilities.h"
#include "parser.h"

/* check for keyword UNITS on line, then scan wavelength or energy units if present */
STATIC void ChkUnits(Parser &p);

/* NB NB NB NB NB NB NB NB NB NB
 *
 * if any "special" save commands are added (commands that copy the file pointer
 * into a separate variable, e.g. like SAVE _DR_), be sure to add that file pointer
 * to SaveFilesInit and CloseSaveFiles !!!
 *
 * SAVE FILE POINTERS SHOULD NEVER BE ALTERED BY ROUTINES OUTSIDE THIS MODULE !!!
 *
 * hence initializations of save 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 */

/* save file header headers - these are written into the string save.chHeader[save.nsave] when
 * the command is parsed
 * save.lgPunHeader[] determines whether header is saved 
 */


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

        char chTemp[MAX_HEADER_SIZE];

        DEBUG_ENTRY( "ParseSave()" );

        /* check that limit not exceeded */
        if( save.nsave >= LIMPUN )
        {
                fprintf( ioQQQ, 
                        "The limit to the number of SAVE options is %ld.  Increase "
                        "LIMPUN in save.h if more are needed.\nSorry.\n", 
                  LIMPUN );
                cdEXIT(EXIT_FAILURE);
        }

        /* initialize this flag, forced true for special cases below (e.g. for FITS files) */
        save.lgSaveToSeparateFiles[save.nsave] = p.nMatch("SEPA");

        /* LAST keyword is an option to save only on last iteration */
        save.lgPunLstIter[save.nsave] = p.nMatch("LAST");

        /* get file name for this save 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( p.GetQuote( chLabel , 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, "ParseSave will not allow save file name %s, please choose another.\nSorry.\n",
                        chLabel);
                cdEXIT(EXIT_FAILURE);
        }
        else if( chLabel[0]=='\0' )
        {
                fprintf( ioQQQ, "ParseSave found a null file name between double quotes, please check command line.\nSorry.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* now copy to chFilename, with optional grid prefix first */
        strcpy( chFilename , save.chGridPrefix.c_str() );
        /* this is optional prefix, normally a null string, set with set save prefix command */
        strcat( chFilename , save.chFilenamePrefix.c_str() );
        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( p.GetQuote( chSecondFilename , false ) )
                lgSecondFilename = false;
        else
                lgSecondFilename = true;

        /* CLOBBER clobber keyword is an option to overwrite rather than
         * append to a given file */
        if( p.nMatch("CLOB") )
        {
                if( p.nMatch(" NO ") )
                {
                        /* do not clobber files */
                        save.lgNoClobber[save.nsave] = true;
                }
                else
                {
                        /* clobber files */
                        save.lgNoClobber[save.nsave] = 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( !p.nMatch("NO TI") && p.nMatch("TITL"))
        {
                sprintf( save.chHeader[save.nsave], 
                        "# %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( p.nMatch("NO HA") )
                save.lgHashEndIter[save.nsave] = false;

        /* save opacity must come first since elements appear as sub-keywords */
        if( p.nMatch("OPAC") )
        {
                /* check for keyword UNITS on line, then scan wavelength or energy units if present,
                 * units are copied into save.chConPunEnr */
                ChkUnits(p);

                strcpy( save.chSave[save.nsave], "OPAC" );

                /* "every" option to save this on every zone -
                 * not present then only last zone is saved */
                if( p.nMatch("EVER" ) )
                {
                        /* save every zone */
                        save.lgSaveEveryZone[save.nsave] = true;
                        save.nSaveEveryZone[save.nsave] = 1;
                }
                else
                {
                        /* only save last zone */
                        save.lgSaveEveryZone[save.nsave] = false;
                        save.nSaveEveryZone[save.nsave] = 1;
                }

                if( p.nMatch("TOTA") )
                {
                        /* DoPunch will call save_opacity to parse the subcommands
                         * save total opacity */
                        strcpy( save.chOpcTyp[save.nsave], "TOTL" );
                        sprintf( save.chHeader[save.nsave], 
                                "#nu/%s\tTot opac\tAbs opac\tScat opac\tAlbedo\telem\n",
                                save.chConPunEnr[save.nsave]);
                }

                else if( p.nMatch("FIGU") )
                {
                        /* do figure for hazy */
                        strcpy( save.chOpcTyp[save.nsave], "FIGU" );
                        sprintf( save.chHeader[save.nsave], 
                                "#nu/%s\tH\tHe\ttot opac\n",
                                save.chConPunEnr[save.nsave] );
                }

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

                        sprintf( save.chHeader[save.nsave], 
                                "#nu/%s\topac\n",
                                save.chConPunEnr[save.nsave] );

                        /* range option - important since so much data - usually want to
                         * only give portion of the continuum */
                        if( p.nMatch("RANGE") ) 
                        {
                                /* get lower and upper range, eventually must be in Ryd */
                                double Energy1 = p.FFmtRead();
                                double Energy2 = p.FFmtRead();
                                if( p.lgEOL() )
                                {
                                        fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                                if( p.nMatch("UNIT" ) )
                                {
                                        // apply units to range option
                                        const char *energyUnits = p.StandardEnergyUnit();
                                        Energy unitChange;
                                        unitChange.set(Energy1, energyUnits );
                                        Energy1 = unitChange.Ryd();
                                        unitChange.set(Energy2, energyUnits );
                                        Energy2 = unitChange.Ryd();
                                }
                                /* get lower and upper rang in Ryd */
                                save.punarg[save.nsave][0] = (realnum)MIN2( Energy1 , Energy2 );
                                save.punarg[save.nsave][1] = (realnum)MAX2( Energy1 , Energy2 );
                                //fprintf(ioQQQ , "DEBUG units change fine %.3e %.3e\n" , save.punarg[save.nsave][0] ,
                                //              save.punarg[save.nsave][1] );
                                //cdEXIT(EXIT_FAILURE);
                        }
                        else
                        {
                                /* these mean full energy range */
                                save.punarg[save.nsave][0] = 0.;
                                save.punarg[save.nsave][1] = 0.;
                        }
                        /* optional last parameter - how many points to bring together */
                        save.punarg[save.nsave][2] = (realnum)p.FFmtRead();

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

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

                else if( p.nMatch("GRAI") )
                {
                        /* save grain opacity command, give optical properties of gains in calculation */
                        strcpy( save.chSave[save.nsave], "DUSO" );
                        /* save grain opacity command in twice, here and above in opacity */
                        sprintf( save.chHeader[save.nsave], 
                                "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" );
                }

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

                else if( p.nMatch("SHEL") )
                {
                        /* save shells, a form of the save opacity command for showing subshell crossections*/
                        strcpy( save.chSave[save.nsave], "OPAC" );

                        /* save subshell cross sections */
                        strcpy( save.chOpcTyp[save.nsave], "SHEL" );

                        /* this is element */
                        save.punarg[save.nsave][0] = (realnum)p.FFmtRead();

                        /* this is ion */
                        save.punarg[save.nsave][1] = (realnum)p.FFmtRead();

                        /* this is shell */
                        save.punarg[save.nsave][2] = (realnum)p.FFmtRead();

                        if( p.lgEOL() )
                        {
                                fprintf( ioQQQ, "There must be atom number, ion, shell\nSorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        sprintf( save.chHeader[save.nsave], 
                                "#sub shell cross section\n" );
                }

                else if( p.nMatch("ELEM") )
                {
                        /* save 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 save_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 = p.GetElem() ) < 0 )
                        {
                                fprintf( ioQQQ, "I did not find an element name on the opacity element command.  Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

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

        /* save H2 has to come early since it has many suboptions */
        else if( p.nMatchErase(" H2 ") )
        {
                /* this is in mole_h2_io.c */
                h2.H2_ParseSave( p , save.chHeader[save.nsave] );
        }

        /* save HD has to come early since it has many suboptions */
        else if( p.nMatchErase(" HD ") )
        {
                /* this is in mole_h2_io.c */
                hd.H2_ParseSave( p , save.chHeader[save.nsave] );
        }

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

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

        else if( p.nMatch(" AGN") )
        {
                /* save tables needed for AGN3 */
                strcpy( save.chSave[save.nsave], " AGN" );
                /* this is the AGN option, to produce a table for AGN */

                /* charge exchange rate coefficients */
                if( p.nMatch("CHAR") )
                {
                        strcpy( save.chSave[save.nsave], "CHAG" );
                        sprintf( save.chHeader[save.nsave], 
                                "#charge exchange rate coefficnt\n" );
                }

                else if( p.nMatch("RECO") )
                {
                        /* save recombination rates for AGN3 table */
                        strcpy( save.chSave[save.nsave], "RECA" );
                        sprintf( save.chHeader[save.nsave], 
                                "#Recom rates for AGN3 table\n" );
                }

                else if( p.nMatch("OPAC") )
                {
                        /* create table for appendix in AGN */
                        strcpy( save.chOpcTyp[save.nsave], " AGN" );
                        strcpy( save.chSave[save.nsave], "OPAC" );
                }

                else if( p.nMatch("HECS") )
                {
                        /* create table for appendix in AGN */
                        strcpy( save.chSaveArgs[save.nsave], "HECS" );
                        sprintf( save.chHeader[save.nsave], 
                                "#AGN3 he cs \n" );
                }

                else if( p.nMatch("HEMI") )
                {
                        /* HEMIS - continuum emission needed for chap 4 of AGN3 */
                        strcpy( save.chSaveArgs[save.nsave], "HEMI" );

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

        else if( p.nMatch("AVER") )
        {
                /* save averages */
                strcpy( save.chSave[save.nsave], "AVER" );
                /* no need to print this standard line of explanation*/
                /*sprintf( save.chHeader[save.nsave], " asserts\n" );*/

                /* actually get the averages from the input stream, and malloc the 
                 * space in the arrays 
                 * save io unit not used in read */
                parse_save_average( p, save.nsave, save.chHeader[save.nsave] );
        }

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

        // save chianti collision strengths in physical units
        else if( p.nMatch("CHIA"))
        {
                strcpy( save.chSave[save.nsave], "CHIA" );

        }

        else if( p.nMatch("CHEM") )
        {

                if( p.nMatch( "RATE" ) )
                {
                        /* >>chng 06 May 30, NPA.  Save reaction rates for selected species */
                        if( lgSecondFilename )
                        {
                                if( p.nMatch( "DEST" ) )
                                        strcpy( save.chSaveArgs[save.nsave], "DEST" );
                                else if( p.nMatch( "CREA" ) )
                                        strcpy( save.chSaveArgs[save.nsave], "CREA" );
                                else if( p.nMatch( "CATA" ) )   
                                        strcpy( save.chSaveArgs[save.nsave], "CATA" );
                                else if( p.nMatch( "ALL" ) )
                                        strcpy( save.chSaveArgs[save.nsave], "ALL " );
                                else
                                        strcpy( save.chSaveArgs[save.nsave], "DFLT" );
                                        
                                strcpy( save.chSave[save.nsave], "CHRT" );
                                save.optname[save.nsave] = chSecondFilename;
                                // Haven't read chemistry database yet, so put off setting up header
                                //sprintf( save.chHeader[save.nsave], "#");  
                        } 

                        else
                        {
                                fprintf(ioQQQ," A species label must appear within a second set of quotes (following the output filename).\n" );
                                fprintf( ioQQQ, " Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                }
        }

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

        else if( p.nMatch("COOL") )
        {
                /* save cooling, actually done by routine cool_save */
                if( p.nMatch("EACH") )
                {
                        strcpy( save.chSave[save.nsave], "EACH");
                        sprintf( save.chHeader[save.nsave], 
                                         "#depth(cm)\tTemp(K)\tCtot(erg/cm3/s)\t" );
                        for( int i = 0 ; i < LIMELM ; i++ )
                        {
                                strcat(save.chHeader[save.nsave], elementnames.chElementSym[i] );
                                strcat(save.chHeader[save.nsave], "\t" );
                        }
                        strcat(save.chHeader[save.nsave], "molecule\tdust\tH2cX\tCT C\tH-fb\tH2ln\tHDro\tH2+ \tFFcm\thvFB\teeff\tComp\tExtr\tExpn\tCycl\tHvin\tdima\n" );
                }
                else
                {
                        strcpy( save.chSave[save.nsave], "COOL");
                        /*>>chng 06 jun 06, revise to be same as save cooling */
                        sprintf( save.chHeader[save.nsave],
                                         "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\tcool fracs\n" );
                }
        }

        // punch the dominant rates for a given species
        else if( p.nMatch("DOMI") && p.nMatch("RATE"))
        {                       
                if( !lgSecondFilename )
                {
                        fprintf( ioQQQ,"This command requires two items in quotes (a filename and a species label).  Only one set of quotes was found.\nSorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
                /* in this case the "second filename" is really the species label. */
                strncpy( save.chSpeciesDominantRates[save.nsave], chSecondFilename, CHARS_SPECIES );

                /* save dominant rates "species" */
                strcpy( save.chSave[save.nsave], "DOMI" );
                sprintf( save.chHeader[save.nsave], 
                        "#depth cm\t%s col cm-2\tsrc s-1\tsnk s-1\n", 
                        save.chSpeciesDominantRates[save.nsave] );
        }

        else if( p.nMatch("DYNA") )
        {
                /* save something dealing with dynamics 
                 * in SaveDo the DYN part of key is used to call DynaPunch,
                 * with the 4th char as its second argument.  DynaSave uses that
                 * 4th letter to decide the job */
                if( p.nMatch( "ADVE" ) )
                {
                        /* save information relating to advection */
                        strcpy( save.chSave[save.nsave], "DYNa");
                        sprintf( save.chHeader[save.nsave], 
                                "#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 SAVE DYNAMICS command.\n" );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

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

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

        else if( p.nMatch("FEII") || p.nMatch("FE II") )
        {
                /* something to do with FeII atom - several options
                 * FeII column densities */
                if( p.nMatch("COLU") && p.nMatch("DENS") )
                {
                        /* save FeII column density */
                        strcpy( save.chSave[save.nsave], "FENl" );

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

                /* FeII continuum, only valid if large atom is set */
                else if( p.nMatch("CONT") )
                {
                        // save FeII continuum, options are total (default), inward,
                        // and outward
                        if( p.nMatch("INWA") )
                        {
                                // inward continuum
                                strcpy( save.chSave[save.nsave], "FEcI" );
                                //sprintf( save.chHeader[save.nsave], 
                                //      "#FeII inward: Wl(A)\tInt[erg cm-2 s-1]\n" );
                        }
                        else if( p.nMatch(" OUT") )
                        {
                                // outward continuum
                                strcpy( save.chSave[save.nsave], "FEcO" );
                                //sprintf( save.chHeader[save.nsave], 
                                //      "#FeII outward: Wl(A)\tInt[erg cm-2 s-1]\n" );
                        }
                        else
                        {
                                // total continuum
                                strcpy( save.chSave[save.nsave], "FEcT" );
                                //sprintf( save.chHeader[save.nsave], 
                                //      "#FeII total: Wl(A)\tInt[erg cm-2 s-1]\n" );
                        }

                        // default units of spectral energy are Ryd, this can change to
                        // many other units
                        ChkUnits(p);

                        // by default give numbers in two columns, row keyword says to
                        // write the numbers across as one long row 
                        if( p.nMatch(" ROW") )
                                save.punarg[save.nsave][0] = 1;
                        else
                                // the default, two columns
                                save.punarg[save.nsave][0] = 2;
                }

                else if( p.nMatch("DEPA") )
                {
                        /* save out departure coefficient for the large FeII atom */
                        sprintf( save.chHeader[save.nsave], 
                                "#FeII departure coefficient \n" );
                        /* optional keyword all means do all levels, if not then just do subset */
                        if( p.nMatch(" ALL") )
                        {
                                /* save all levels, calls routine FeIIPunDepart */
                                strcpy( save.chSave[save.nsave], "FE2D" );
                        }
                        else
                        {
                                /* save a very few selected levels, calls routine FeIIPunDepart */
                                strcpy( save.chSave[save.nsave], "FE2d" );
                        }
                }

                else if( p.nMatch("LEVE") )
                {
                        /* save level energies and statistical weights for large FeII atom */
                        sprintf( save.chHeader[save.nsave], 
                                "#FeII energy(wn)\tstat weight\n" );
                        strcpy( save.chSave[save.nsave], "FE2l" );
                }

                else if( p.nMatch("LINE") )
                {
                        /* save FeII lines command
                         * three optional parameters, threshold for faintest
                         * line to print, lower and upper energy bounds */

                        /* save intensities from large FeII atom */
                        strcpy( save.chSave[save.nsave], "FEli" );

                        /* short keyword makes save half as big */
                        if( p.nMatch("SHOR") )
                        {
                                FeII.lgShortFe2 = true;
                        }
                        else
                        {
                                FeII.lgShortFe2 = false;
                        }

                        /* first optional number changes the threshold of weakest line to print*/
                        /* fe2thresh is intensity relative to normalization line,
                        * normally Hbeta, and is set to zero in zero.c */
                        FeII.fe2thresh = (realnum)p.FFmtRead();
                        if( p.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 saved,
                         * this is to limit size of output file */
                        FeII.fe2ener[0] = (realnum)p.FFmtRead();
                        if( p.lgEOL() )
                        {
                                FeII.fe2ener[0] = 0.;
                        }

                        FeII.fe2ener[1] = (realnum)p.FFmtRead();
                        if( p.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 FeIISaveLines routine
                         * that lives in the FeIILevelPops file */
                        sprintf( save.chHeader[save.nsave], 
                                "#FeII ipHi\tipLo\tWL(A)\tlog I\tI/Inorm\t\tTau\n" );
                }

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

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

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

                        /* optional keyword all means do all levels, if not then just do subset */
                        if( p.nMatch(" ALL") )
                        {
                                /* save all levels, calls routine FeIIPunPop */
                                strcpy( save.chSave[save.nsave], "FE2P" );
                                save.punarg[save.nsave][0] = 0.;
                                save.punarg[save.nsave][1] = (realnum)NFE2LEVN;
                        }

                        /* optional keyword range means read lower and upper bound, do these */
                        else if( p.nMatch("RANG") )
                        {
                                /* save range of levels, calls routine FeIIPunPop */
                                strcpy( save.chSave[save.nsave], "FE2P" );
                                save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
                                save.punarg[save.nsave][1] = (realnum)p.FFmtRead();
                                if( p.lgEOL() || save.punarg[save.nsave][0] <0 ||
                                        save.punarg[save.nsave][0]>= save.punarg[save.nsave][1] )
                                {
                                        fprintf( ioQQQ, "There must be two numbers on this save "
                                                "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, %g, must be less than the second, %g.\n",
                                                save.punarg[save.nsave][0],
                                                save.punarg[save.nsave][1]);
                                        fprintf( ioQQQ, "Sorry.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }

                        else
                        {
                                /* save a very few selected levels, calls routine FeIIPunPop */
                                strcpy( save.chSave[save.nsave], "FE2p" );
                        }
                }
                else
                {
                        fprintf( ioQQQ, "There must be a second keyword on this SAVE 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 save continuum command, with many options,
         * the first 3 char of the chSave flag will always be "CON" 
         * with the last indicating which one */
        else if( p.nMatch("CONT") && !p.nMatch("XSPE") )
        {
                /* this flag is checked in PrtComment to generate a caution
                 * if continuum is saved but iterations not performed */
                save.lgPunContinuum = true;

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

                if( p.nMatch("BINS") )
                {
                        /* continuum binning */
                        strcpy( save.chSave[save.nsave], "CONB" );

                        sprintf( save.chHeader[save.nsave], 
                                "#Continuum binning enrOrg/%s\tEnergy\twidth of cells\n",
                                save.chConPunEnr[save.nsave] );
                }

                else if( p.nMatch("DIFF") )
                {
                        /* diffuse continuum, the locally emitted lines and continuum */
                        strcpy( save.chSave[save.nsave], "COND" );

                        /* by default gives lines and continuum separately only for
                         * last zone.  The keyword ZONE says to give the total for every
                         * zone in one very low row */
                        if( p.nMatch("ZONE") )
                        {
                                sprintf( save.chHeader[save.nsave], 
                                        "#energy/%s then emission per zone\n",
                                        save.chConPunEnr[save.nsave] );
                                save.punarg[save.nsave][0] = 2.;

                        }
                        else
                        {
                                sprintf( save.chHeader[save.nsave], 
                                        "#energy/%s\tConEmitLocal\tDiffuseLineEmission\tTotal\n",
                                        save.chConPunEnr[save.nsave] );
                                save.punarg[save.nsave][0] = 1.;
                        }
                }

                else if( p.nMatch("EMIS") )
                {
                        /* continuum volume emissivity and opacity as a function of radius */
                        strcpy( save.chSave[save.nsave], "CONS" );

                        double num = p.FFmtRead();
                        if( p.lgEOL() )
                                p.NoNumb( "continuum emissivity frequency" );
                        save.emisfreq.set( num, save.chConPunEnr[save.nsave] );
                        if( save.emisfreq.Ryd() < rfield.emm || save.emisfreq.Ryd() > rfield.egamry )
                        {
                                fprintf( ioQQQ, " The frequency is outside the Cloudy range\n Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        sprintf( save.chHeader[save.nsave], 
                                 "#Radius\tdepth\tnujnu\tkappa_abs\tkappa_sct\n" );
                }

                else if( p.nMatch("EMIT") )
                {
                        /* continuum emitted by cloud */
                        strcpy( save.chSave[save.nsave], "CONE" );

                        sprintf( save.chHeader[save.nsave], 
                                "#Energy/%s\treflec\toutward\ttotal\tline\tcont\n",
                                save.chConPunEnr[save.nsave] );
                }

                else if( p.nMatch("FINE" ) )
                {
                        rfield.lgSaveOpacityFine = true;
                        /* fine transmitted continuum cloud */
                        strcpy( save.chSave[save.nsave], "CONf" );

                        sprintf( save.chHeader[save.nsave], 
                                "#Energy/%s\tTransmitted\n",
                                save.chConPunEnr[save.nsave] );

                        /* range option - important since so much data */
                        if( p.nMatch("RANGE") ) 
                        {
                                /* get lower and upper range, eventually must be in Ryd */
                                double Energy1 = p.FFmtRead();
                                double Energy2 = p.FFmtRead();
                                if( p.lgEOL() )
                                {
                                        fprintf(ioQQQ,"There must be two numbers, the lower and upper energies in Ryd.\nSorry.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                                if( p.nMatch("UNIT" ) )
                                {
                                        // apply units to range option
                                        const char *energyUnits = p.StandardEnergyUnit();
                                        Energy unitChange;
                                        unitChange.set(Energy1, energyUnits );
                                        Energy1 = unitChange.Ryd();
                                        unitChange.set(Energy2, energyUnits );
                                        Energy2 = unitChange.Ryd();
                                }
                                /* get lower and upper rang in Ryd */
                                save.punarg[save.nsave][0] = (realnum)MIN2( Energy1 , Energy2 );
                                save.punarg[save.nsave][1] = (realnum)MAX2( Energy1 , Energy2 );
                                //fprintf(ioQQQ , "DEBUG units change fine %.3e %.3e\n" , save.punarg[save.nsave][0] ,
                                //              save.punarg[save.nsave][1] );
                                //cdEXIT(EXIT_FAILURE);
                        }
                        else
                        {
                                /* these mean full energy range */
                                save.punarg[save.nsave][0] = 0.;
                                save.punarg[save.nsave][1] = 0.;
                        }
                        /* optional last parameter - how many points to bring together */
                        save.punarg[save.nsave][2] = (realnum)p.FFmtRead();

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

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

                else if( p.nMatch("GRAI") )
                {
                        /* save grain continuum in optically thin limit */
                        strcpy( save.chSave[save.nsave], "CONG" );

                        sprintf( save.chHeader[save.nsave], 
                                "#energy\tgraphite\trest\ttotal\n" );
                }

                else if( p.nMatch("INCI") )
                {
                        /* incident continuum */
                        strcpy( save.chSave[save.nsave], "CONC" );

                        sprintf( save.chHeader[save.nsave], 
                                "#Incident Continuum, Enr\tnFn \n" );
                }

                else if( p.nMatch("INTE") )
                {
                        /* continuum interactions */
                        strcpy( save.chSave[save.nsave], "CONi" );

                        sprintf( save.chHeader[save.nsave], 
                                "#Continuum interactions, inc, otslin. otscon, ConInterOut, outlin \n" );
                        /* this is option for lowest energy, if nothing then zero */
                        save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
                }

                else if( p.nMatch("IONI") )
                {
                        /* save ionizing continuum*/
                        strcpy( save.chSave[save.nsave], "CONI" );

                        /* this is option for lowest energy, if nothing then zero */
                        save.punarg[save.nsave][0] = (realnum)p.FFmtRead();

                        /* this is option for smallest interaction to save, def is 1 percent */
                        save.punarg[save.nsave][1] = (realnum)p.FFmtRead();
                        if( p.lgEOL() )
                                save.punarg[save.nsave][1] = 0.01f;

                        /* "every" option to save this on every zone -
                         * not present then only last zone is saved */
                        if( p.nMatch("EVER" ) )
                        {
                                /* save every zone */
                                save.lgSaveEveryZone[save.nsave] = true;
                                save.nSaveEveryZone[save.nsave] = 1;
                        }
                        else
                        {
                                /* only save last zone */
                                save.lgSaveEveryZone[save.nsave] = false;
                                save.nSaveEveryZone[save.nsave] = 1;
                        }

                        /* put the header at the top of the file */
                        sprintf( save.chHeader[save.nsave], 
                                "#cell(on C scale)\tnu\tflux\tflx*cs\tFinc\totsl\totsc\toutlin\toutcon\trate/tot\tintegral\tline\tcont\n" );
                }
#ifdef USE_NLTE7
                else if( p.nMatch("NLTE") )
                {
                        /* continuum emitted by cloud */
                        strcpy( save.chSave[save.nsave], "CONl" );

                        sprintf( save.chHeader[save.nsave],
                                "   spectrum1   NeXY6   XNUMX\n");
                }
#endif

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

                else if( p.nMatch("TRAN") )
                {
                        /* transmitted continuum */
                        strcpy( save.chSave[save.nsave], "CONT" );

                        sprintf( save.chHeader[save.nsave], 
                                "#ener\tTran Contin\ttrn coef\n" );
                }

                else if( p.nMatch(" TWO") )
                {
                        /* total two photon continua rfield.TotDiff2Pht */
                        strcpy( save.chSave[save.nsave], "CON2" );

                        sprintf( save.chHeader[save.nsave], 
                                "#energy\t n_nu\tnuF_nu \n" );
                }

                else if( p.nMatch(" RAW") )
                {
                        /* "raw" continua */
                        strcpy( save.chSave[save.nsave], "CORA" );

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

                else if( p.nMatch("REFL") )
                {
                        /* reflected continuum */
                        strcpy( save.chSave[save.nsave], "CONR" );

                        sprintf( save.chHeader[save.nsave], 
                                "#Reflected\tcont\tline\ttotal\talbedo\tConID\n" );
                }

                else
                {
                        /* this is the usual save continuum command,
                         * ipType is index for continuum array to set either
                         * iteration or cumulative output */
                        int ipType = 0;
                        if( p.nMatch( "CUMU" ) )
                                ipType = 1;
                        save.punarg[save.nsave][0] = (realnum)ipType;
                        strcpy( save.chSave[save.nsave], "CON " );
                        char chHold[100];
                        strcpy( chHold, "#Cont " );
                        if( ipType > 0 )
                                strcpy( chHold , "#Cumul " );
                        sprintf( save.chHeader[save.nsave], 
                                "%s nu\tincident\ttrans\tDiffOut\tnet trans\treflc\ttotal\treflin\toutlin\tlineID\tcont\tnLine\n" ,
                                chHold );

                        /* >>chng 06 apr 03, add "every" option to save this on every zone -
                         * if every is not present then only last zone is saved */
                        if( p.nMatch("EVER" ) )
                        {
                                /* save every zone */
                                save.lgSaveEveryZone[save.nsave] = true;
                                /* option to say how many to skip */
                                save.nSaveEveryZone[save.nsave] = (long)p.FFmtRead();
                                if( p.lgEOL() )
                                        save.nSaveEveryZone[save.nsave] = 1;
                        }
                        else
                        {
                                /* only save last zone */
                                save.lgSaveEveryZone[save.nsave] = false;
                                save.nSaveEveryZone[save.nsave] = 1;
                        }
                }
        }

        /* save 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( p.nMatch("CONV") )
        {
                if( p.nMatch("REAS") )
                {
                        /* this does not count as a save option (really) */
                        save.lgPunConv = true;
                        /* this is done below */
                        strcpy( save.chSave[save.nsave], "" );
                        save.lgRealSave[save.nsave] = false;
                }
                else if( p.nMatch("ERRO") )
                {
                        /* save zone by zone errors in pressure, electron density, and heating-cooling */
                        /* convergence error */
                        strcpy( save.chSave[save.nsave], "CNVE" );
                        sprintf( save.chHeader[save.nsave], 
                                "#depth\tnPres2Ioniz\tP(cur)\tP%%error\tNE(cor)\tNE(cur)\tNE%%error\tHeat\tCool\tHC%%error\n" );
                }
                else if( p.nMatch("BASE") )
                {
                        /* save converged quantities in Converge base for each pass through
                         * solvers - not one pass per zone */
                        strcpy( save.chSave[save.nsave], "CNVB" );
                        strcpy( save.chSave[save.nsave], "" );
                        save.lgRealSave[save.nsave] = false;
                }
                else
                {
                        fprintf( ioQQQ, "There must be a second keyword on this command.\n" );
                        fprintf( ioQQQ, "The ones I know about are REASON, ERROR, and BASE.\n" );
                        fprintf( ioQQQ, "Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( p.nMatch(" DR ") )
        {
                /* first occurrence of save dr to follow choice in change of zoning */
                save.lgDROn = true;
                strcpy( save.chSave[save.nsave], "" );
                save.lgRealSave[save.nsave] = false;
        }

        else if( p.nMatch("ELEM") && !p.nMatch("GAMMA") && !p.nMatch("COOL") ) // do not trip on SAVE COOLING EACH ELEMENT
        {
                /* option to save ionization structure of some element
                 * will give each stage of ionization, vs depth */
                strcpy( save.chSave[save.nsave], "ELEM" );

                /* this returns element number on c scale */
                /* >>chng 04 nov 23, had converted to f scale, leave on c */
                nelem = p.GetElem();
                if( nelem < 0 || nelem >= LIMELM )
                {
                        fprintf( ioQQQ, " I could not recognize a valid element name on this line.\n" );
                        fprintf( ioQQQ, " Please check your input script. Bailing out...\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* this is the atomic number on the c scale */
                save.punarg[save.nsave][0] = (realnum)nelem;

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

                /* start printing header line - first will be the depth in cm */
                sprintf( save.chHeader[save.nsave], "#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( save.chHeader[save.nsave], 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( save.chHeader[save.nsave], chTemp );
                }
                if( nelem==ipCARBON )
                {
                        sprintf( chTemp, "\tC1\tC1*\tC1**\tC2\tC2*\tCO");
                        strcat( save.chHeader[save.nsave], chTemp );
                }
                else if( nelem==ipOXYGEN )
                {
                        sprintf( chTemp, "\tO1\tO1*\tO1**");
                        strcat( save.chHeader[save.nsave], chTemp );
                }

                /* finally the new line */
                strcat( save.chHeader[save.nsave], "\n");
        }

        else if( p.nMatch("FITS") )
        {

#ifdef FLT_IS_DBL
                fprintf( ioQQQ, "Saving 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 */
                save.lgFITS[save.nsave] = true;
                /* concatenating files in a grid run would be illegal FITS */
                save.lgSaveToSeparateFiles[save.nsave] = true;
                save.lgPunLstIter[save.nsave] = true;
                save.FITStype[save.nsave] = NUM_OUTPUT_TYPES;

                strcpy( save.chSave[save.nsave], "FITS" );
#endif

        }

        else if( p.nMatch("FRED") )
        {
                /* save out some stuff for Fred's dynamics project */
                sprintf( save.chHeader[save.nsave], 
                        "#Radius\tDepth\tVelocity(km/s)\tdvdr(cm/s)\thden\teden\tTemperature\tRadAccel line\tRadAccel con\t"
                        "Force multiplier\ta(e thin)\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\tOVI(1034) TauIn\tTauCon\n");

                strcpy( save.chSave[save.nsave], "FRED" );
        }

        else if( p.nMatch("GAMM") )
        {
                /* save all photoionization rates for all subshells */
                sprintf( save.chHeader[save.nsave], 
                        "#Photoionization rates \n" );
                if( p.nMatch("ELEMENT") )
                {
                        /* element keyword, find element name and stage of ionization, 
                         * will print photoionization rates for valence of that element */
                        strcpy( save.chSave[save.nsave], "GAMe" );

                        /* this returns element number on c scale */
                        nelem = p.GetElem();
                        /* this is the atomic number on the C scale */
                        save.punarg[save.nsave][0] = (realnum)nelem;

                        /* this will become the ionization stage on C scale */
                        save.punarg[save.nsave][1] = (realnum)p.FFmtRead() - 1;
                        if( p.lgEOL() )
                                p.NoNumb("element ionization stage" );
                        if( save.punarg[save.nsave][1]<0 || save.punarg[save.nsave][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( save.chSave[save.nsave], "GAMt" );
                }

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

                        strcpy( save.chSave[save.nsave], "DUSO" );
                        /* save grain opacity command in twice, here and above in opacity */
                        sprintf( save.chHeader[save.nsave], 
                                "#grain\tnu/%s\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n",
                                save.chConPunEnr[save.nsave] );
                }
                else if( p.nMatch("ABUN") )
                {
                        /* save grain abundance */
                        strcpy( save.chSave[save.nsave], "DUSA" );
                        sprintf( save.chHeader[save.nsave], 
                                 "#grain\tdepth\tabundance (g/cm^3)\n" );
                }
                else if( p.nMatch("D/G ") )
                {
                        /* save grain dust/gas mass ratio */
                        strcpy( save.chSave[save.nsave], "DUSD" );
                        sprintf( save.chHeader[save.nsave], 
                                 "#grain\tdepth\tdust/gas mass ratio\n" );
                }
                else if( p.nMatch("PHYS") )
                {
                        /* save grain physical conditions */
                        strcpy( save.chSave[save.nsave], "DUSP" );
                        sprintf( save.chHeader[save.nsave], 
                                "#grain\tdepth\tpotential\n" );
                }
                else if( p.nMatch(" QS ") )
                {
                        strcpy( save.chSave[save.nsave], "DUSQ" );
                        sprintf( save.chHeader[save.nsave], 
                                "#grain\tnu\tQ_abs\tQ_scat\n" );
                }
                else if( p.nMatch("TEMP") )
                {
                        /* save temperatures of each grain species */
                        strcpy( save.chSave[save.nsave], "DUST" );
                        /* cannot save grain labels since they are not known yet */
                        sprintf( save.chHeader[save.nsave], 
                                "#grain temperature\n" );
                }
                else if( p.nMatch("DRIF") )
                {
                        /* save drift velocity of each grain species */
                        strcpy( save.chSave[save.nsave], "DUSV" );
                        /* cannot save grain labels since they are not known yet */
                        sprintf( save.chHeader[save.nsave], 
                                "#grain drift velocity\n" );
                }
                else if( p.nMatch("EXTI") )
                {
                        /* save grain extinction */
                        strcpy( save.chSave[save.nsave], "DUSE" );
                        /* cannot save grain labels since they are not known yet */
                        sprintf( save.chHeader[save.nsave], 
                                "#depth\tA_V(extended)\tA_V(point)\n" );
                }
                else if( p.nMatch("CHAR") )
                {
                        /* save charge per grain (# elec/grain) for each grain species */
                        strcpy( save.chSave[save.nsave], "DUSC" );
                        /* cannot save grain labels since they are not known yet */
                        sprintf( save.chHeader[save.nsave], 
                                "#grain charge\n" );
                }
                else if( p.nMatch("HEAT") )
                {
                        /* save heating due to each grain species */
                        strcpy( save.chSave[save.nsave], "DUSH" );
                        /* cannot save grain labels since they are not known yet */
                        sprintf( save.chHeader[save.nsave], 
                                "#grain heating\n" );
                }
                else if( p.nMatch("POTE") )
                {
                        /* save floating potential of each grain species */
                        strcpy( save.chSave[save.nsave], "DUSP" );
                        /* cannot save grain labels since they are not known yet */
                        sprintf( save.chHeader[save.nsave], 
                                "#grain\tdepth\tpotential\n" );
                }
                else if( p.nMatch("H2RA") )
                {
                        /* save grain H2rate - H2 formation rate for each type of grains */
                        strcpy( save.chSave[save.nsave], "DUSR" );
                        /* cannot save grain labels since they are not known yet */
                        sprintf( save.chHeader[save.nsave], 
                                "#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( p.nMatch("GAUN") )
        {
                strcpy( save.chSave[save.nsave], "GAUN" );
                sprintf( save.chHeader[save.nsave], 
                        "#Gaunt factors.\n" );
        }
        else if( p.nMatch("GRID") )
        {
                strcpy( save.chSave[save.nsave], "GRID" );
                /* automatically generate no hash option */
                save.lgHashEndIter[save.nsave] = false;
        }
        else if( p.nMatch( "HIST" ) )
        {
                /* save pressure history of current zone */
                if( p.nMatch( "PRES") )
                {
                        /* save pressure history - density - pressure for this zone */
                        strcpy( save.chSave[save.nsave], "HISp" );
                        sprintf( save.chHeader[save.nsave], 
                                "#iter zon\tdensity\tpres cur\tpres error\n" );
                }
                /* save temperature history of current zone */
                else if( p.nMatch( "TEMP" ) )
                {
                        /* save pressure history - density - pressure for this zone */
                        strcpy( save.chSave[save.nsave], "HISt" );
                        sprintf( save.chHeader[save.nsave], 
                                "#iter zon\ttemperature\theating\tcooling\n" );
                }
        }

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

        /* QHEAT has to come before HEAT... */
        else if( p.nMatch("QHEA") )
        {
                /* 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( p.nMatch("HEAT") )
        {
                /* save heating */
                strcpy( save.chSave[save.nsave], "HEAT" );
                /*>>chng 06 jun 06, revise to be same as save cooling */
                sprintf( save.chHeader[save.nsave], 
                        "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\theat fracs\n" );
        }

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

        else if( p.nMatch("HUMM") )
        {
                strcpy( save.chSave[save.nsave], "HUMM" );
                sprintf( save.chHeader[save.nsave], 
                        "#input to DHs routine.\n" );
        }

        else if( p.nMatch("HYDR") )
        {
                /* save hydrogen physical conditions */
                if( p.nMatch("COND") )
                {
                        strcpy( save.chSave[save.nsave], "HYDc" );
                        sprintf( save.chHeader[save.nsave], 
                                "#depth\tTe\tHDEN\tEDEN\tHI/H\tHII/H\tH2/H\tH2+/H\tH3+/H\tH-/H\n" );
                        /* save hydrogen ionization */
                }

                /* save information on 21 cm excitation processes - accept either keyword 21cm or 21 cm */
                else if( p.nMatch("21 CM") ||p.nMatch("21CM"))
                {
                        /* save information about 21 cm line */
                        strcpy( save.chSave[save.nsave], "21CM" );
                        sprintf( save.chHeader[save.nsave], 
                                "#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( p.nMatch("IONI") )
                {
                        /* save hydrogen ionization */
                        strcpy( save.chSave[save.nsave], "HYDi" );
                        sprintf( save.chHeader[save.nsave], 
                                "#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( p.nMatch("POPU") )
                {
                        /* save hydrogen populations */
                        strcpy( save.chSave[save.nsave], "HYDp" );
                        sprintf( save.chHeader[save.nsave], 
                                "#depth\tn(H0)\tn(H+)\tn(1s)\tn(2s)\tn(2p)\tetc\n" );
                }
                else if( p.nMatch("LINE") )
                {
                        /* save hydrogen lines
                         * hydrogen line intensities and optical depths  */
                        strcpy( save.chSave[save.nsave], "HYDl" );
                        sprintf( save.chHeader[save.nsave], 
                                "#nHi\tlHi\tnLo\tlLo\tE(ryd)\ttau\n" );
                }
                else if( p.nMatch(" LYA") )
                {
                        /* save hydrogen Lya some details about Lyman alpha  */
                        strcpy( save.chSave[save.nsave], "HYDL" );
                        sprintf( save.chHeader[save.nsave], 
                                "#depth\tTauIn\tTauTot\tn(2p)/n(1s)\tTexc\tTe\tTex/T\tPesc\tPdes\tpump\topacity\talbedo\n" );
                }
                else
                {
                        fprintf( ioQQQ, "Save hydrogen has options: CONDitions, 21 CM, LINE, POPUlations, and IONIzation.\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( p.nMatch("IONI") )
        {
                if( p.nMatch("RATE") )
                {
                        /* save ionization rates, search for the name of an element */
                        if( (nelem = p.GetElem() ) < 0 )
                        {
                                fprintf( ioQQQ, "There must be an element name on the ionization rates command.  Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        save.punarg[save.nsave][0] = (realnum)nelem;
                        strcpy( save.chSave[save.nsave], "IONR" );
                        sprintf( save.chHeader[save.nsave], 
                                "#%s depth\teden\tdynamics.Rate\tabund\tTotIonize\tTotRecom\tSource\t ... \n",
                                elementnames.chElementSym[nelem]);
                }
                else
                {
                        /* save table giving ionization means */
                        strcpy( save.chSave[save.nsave], "IONI" );
                        sprintf( save.chHeader[save.nsave], 
                                "#Mean ionization distribution\n" );
                }
        }

        else if( p.nMatch(" IP ") )
        {
                strcpy( save.chSave[save.nsave], " IP " );
                sprintf( save.chHeader[save.nsave], 
                        "#ionization potentials, valence shell\n" );
        }

        else if( p.nMatch("LEID") )
        {
                if( p.nMatch( "LINE" ) )
                {
                        /* save Leiden lines
                         * final intensities of the Leiden PDR models */
                        strcpy( save.chSave[save.nsave], "LEIL" );
                        sprintf( save.chHeader[save.nsave], "#ion\twl\tInt\trel int\n");
                }
                else
                {
                        /* save Leiden structure
                        * structure of the Leiden PDR models */
                        strcpy( save.chSave[save.nsave], "LEIS" );
                        sprintf( save.chHeader[save.nsave], 
                                /* 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)\tN(OH)\tN(e)\tN(He+)\tN(H+)\tN(H3+)\t"
                                /* 31 - 32 */
                                "H2(Sol)\tH2(FrmGrn)\tH2(photodiss)\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" );
                }
        }


        // save results for NLTE series of plasma comparison meetings, 
        // specifically NLTE7 2011 Dec
        else if( p.nMatch("NLTE") )
        {
# ifdef USE_NLTE7
                strcpy( save.chSave[save.nsave], "NLTE" );
# else
                fprintf(ioQQQ," PROBLEM You must enable the USE_NLTE7 macro at compile-time to use this command.\n");
                fprintf(ioQQQ," To do so, add EXTRA=\"-DUSE_NLTE7\" to the end of the make command.\n");
                fprintf(ioQQQ," An example for a quad core machine:\n make -j 4 EXTRA=\"-DUSE_NLTE7\" \n");
                fprintf(ioQQQ," in the sys_XXX folder that you want to use.\n\n\n");
                cdEXIT(EXIT_FAILURE);
# endif
        }

        /* FE2NRG, FE2TP, and FE2COLL write the internal Fe II data into Stout format */
        else if (p.nMatch("FE2NRG"))
        {
                strcpy( save.chSave[save.nsave], "LY1" );
        }

        else if (p.nMatch("FE2TP"))
        {
                strcpy( save.chSave[save.nsave], "LY2" );
        }

        else if (p.nMatch("FE2COLL"))
        {
                strcpy( save.chSave[save.nsave], "LY3" );
        }

        else if( (p.nMatch("LINE") && p.nMatch("LIST")) || p.nMatch("LINELIST") )
        {
                /* save line list "output file" "Line List file" */
                strcpy( save.chSave[save.nsave], "LLST" );

                /* 
                 * 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 SAVE 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( save.ipPnunit[save.nsave] == NULL )
                {
                        /* make sure we free any allocated space from a previous call */
                        save.SaveLineListFree(save.nsave);

                        save.nLineList[save.nsave] = cdGetLineList(chSecondFilename,
                                                                   save.chLineListLabel[save.nsave],
                                                                   save.wlLineList[save.nsave]);

                        if( save.nLineList[save.nsave] < 0 )
                        {
                                fprintf(ioQQQ,"DISASTER could not open SAVE LINE LIST file %s \n",
                                        chSecondFilename );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                // check whether intrinsic or emergent line emissivity
                save.lgEmergent[save.nsave] = false;
                if( p.nMatch("EMER") )
                        save.lgEmergent[save.nsave] = true;

                // check whether cumulative or specific line emission
                save.lgCumulative[save.nsave] = false;
                if( p.nMatch("CUMU") )
                        save.lgCumulative[save.nsave] = true;

                /* ratio option, in which pairs of lines form ratios, first over
                 * second */
                if( p.nMatch("RATI") )
                {
                        save.lgLineListRatio[save.nsave] = true;
                        if( save.nLineList[save.nsave]%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", save.nLineList[save.nsave]);
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        /* no ratio */
                        save.lgLineListRatio[save.nsave] = false;
                }

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

                // check whether column or row (default)
                if( p.nMatch("COLUMN") )
                {
                        save.punarg[save.nsave][1] = 1;
                }
                else
                {
                        save.punarg[save.nsave][1] = 0;
                }

                /* give header line */
                sprintf( save.chHeader[save.nsave], "#lineslist" );
                // do header now if reporting rows of lines
                if( !save.punarg[save.nsave][1] )
                {
                        for( long int j=0; j<save.nLineList[save.nsave]; ++j )
                        {
                                /* if taking ratio then put div sign between pairs */
                                if( save.lgLineListRatio[save.nsave] && is_odd(j) )
                                        strcat( save.chHeader[save.nsave] , "/" );
                                else
                                        strcat( save.chHeader[save.nsave] , "\t" );
                                sprintf( chTemp, "%s ", save.chLineListLabel[save.nsave][j] );
                                strcat( save.chHeader[save.nsave], chTemp );
                                sprt_wl( chTemp, save.wlLineList[save.nsave][j] );
                                strcat( save.chHeader[save.nsave], chTemp );
                        }
                }
                strcat( save.chHeader[save.nsave], "\n" );
        }

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

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

                else if( p.nMatch("PRES") )
                {
                        /* save contributors to line pressure */
                        strcpy( save.chSave[save.nsave], "PREL" );
                        sprintf( save.chHeader[save.nsave], 
                                "#P depth\tPtot\tPline/Ptot\tcontributors to line pressure\n" );
                }

                else if( p.nMatch("EMIS") )
                {
                        /* this used to be the save lines structure command, is now
                         * the save lines emissivity command 
                         * give line emissivity vs depth */
                        // check whether intrinsic or emergent line emissivity
                        save.lgEmergent[save.nsave] = false;
                        if( p.nMatch("EMER") )
                                save.lgEmergent[save.nsave] = true;
                        strcpy( save.chSave[save.nsave], "LINS" );
                        sprintf( save.chHeader[save.nsave],
                                "#");
                        /* read in the list of lines to examine */
                        parse_save_line(p,false, chTemp );
                        strcat( save.chHeader[save.nsave], chTemp );
                }

                else if( p.nMatch(" RT " ) )
                {
                        /* save line RT */
                        strcpy( save.chSave[save.nsave], "LINR" );
                        /* save some details needed for line radiative transfer 
                         * routine in save_line.cpp */
                        Parse_Save_Line_RT(p);
                }

                else if( p.nMatch("CUMU") )
                {
                        bool lgEOL;
                        /* save lines cumulative 
                         * this will be integrated line intensity, function of depth */
                        strcpy( save.chSave[save.nsave], "LINC" );
                        // option for intrinsic (default) or emergent
                        save.lgEmergent[save.nsave] = false;
                        if( p.nMatch("EMER") )
                                save.lgEmergent[save.nsave] = true;
                        /* option for either relative intensity or abs luminosity */
                        if( p.nMatch("RELA") )
                        {
                                lgEOL = true;
                                sprintf( save.chHeader[save.nsave], "#" );
                        }
                        else
                        {
                                sprintf( save.chHeader[save.nsave], "#" );
                                lgEOL = false;
                        }
                        /* read in the list of lines to examine */
                        parse_save_line(p, lgEOL, chTemp );
                        strcat( save.chHeader[save.nsave], chTemp );
                }

                else if( p.nMatch("DATA") )
                {
                        /* save line data, done in SaveLineData */

                        /* the default will be to make wavelengths like in the printout, called labels,
                         * if units appears then other units will be used instead */
                        save.chConPunEnr[save.nsave] = "labl";

                        /* check for keyword UNITS on line, then scan wavelength or energy units if present,
                         * units are copied into save.chConPunEnr */
                        if( p.nMatch("UNIT") )
                                ChkUnits(p);
                        strcpy( save.chSave[save.nsave], "LIND" );
                        sprintf( save.chHeader[save.nsave], 
                                "#Emission line data.\n" );
                }

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

                else if( p.nMatch("LABE") )
                {
                        /* save line labels */
                        strcpy( save.chSave[save.nsave], "LINL" );
                        sprintf( save.chHeader[save.nsave], 
                                "#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( p.nMatch("LONG") )
                                save.punarg[save.nsave][0] = 1;
                        else
                                save.punarg[save.nsave][0] = 0;
                }

                else if( p.nMatch("OPTI") )
                {
                        /* save line optical depths, done in SaveLineStuff */
                        strcpy( save.chSave[save.nsave], "LINO" );

                        /* the default will be to make wavelengths line in the printout, called labels,
                         * if units appears then other units will be used instead */
                        save.chConPunEnr[save.nsave] = "labl";

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

                        sprintf( save.chHeader[save.nsave], 
                                "#species\tenergy/%s\topt depth\tdamp\n",
                                save.chConPunEnr[save.nsave] );

                        /* this is optional limit to smallest optical depths */
                        save.punarg[save.nsave][0] = (realnum)pow(10.,p.FFmtRead());
                        /* this is default of 0.1 napier */
                        if( p.lgEOL() )
                        {
                                save.punarg[save.nsave][0] = 0.1f;
                        }
                }

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

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

                        if( p.nMatch(" OFF") )
                        {
                                /* no lower limit - print all lines */
                                save.punarg[save.nsave][0] = -1.f;
                        }
                }

                else if( p.nMatch("INTE") )
                {
                        /* this will be full set of line intensities */
                        strcpy( save.chSave[save.nsave], "LINI" );
                        sprintf( save.chHeader[save.nsave], 
                                "#Emission line intrinsic intensities per unit inner area\n" );
                        if( p.nMatch("COLU") )
                                /* column is key to save single column */
                                strcpy( save.chPunRltType, "column" );
                        else
                                /* array is key to save large array */
                                strcpy( save.chPunRltType, "array " );

                        save.punarg[save.nsave][0] = 0.;
                        // ALL option - all lines, even zero intensities
                        if( p.nMatch( " ALL" ) )
                                save.punarg[save.nsave][0] = -1.;

                        // check whether intrinsic or emergent line emissivity
                        save.lgEmergent[save.nsave] = false;
                        if( p.nMatch("EMER") )
                                save.lgEmergent[save.nsave] = true;

                        if( p.nMatch("EVER") )
                        {
                                save.LinEvery = (long int)p.FFmtRead();
                                save.lgLinEvery = true;
                                if( p.lgEOL() )
                                {
                                        fprintf( ioQQQ, 
                                                "There must be a second number, the number of zones to print.\nSorry.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else
                        {
                                save.LinEvery = geometry.nend[0];
                                save.lgLinEvery = false;
                        }
                }
                else
                {
                        fprintf( ioQQQ, 
                                "This option for SAVE LINE is something that I do not understand.  Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( p.nMatch(" MAP") )
        {
                strcpy( save.chSave[save.nsave], "MAP " );
                sprintf( save.chHeader[save.nsave], 
                        "#te, heating, cooling.\n" );
                /* do cooling space map for specified zones
                 * if no number, or <0, do map and save out without doing first zone
                 * does map by calling punt(" map") 
                 */
                hcmap.MapZone = (long)p.FFmtRead();
                if( p.lgEOL() )
                {
                        hcmap.MapZone = 1;
                }

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

                        hcmap.RangeMap[1] = (realnum)p.FFmtRead();
                        if( lgLogOn )
                                hcmap.RangeMap[1] = (realnum)pow((realnum)10.f,hcmap.RangeMap[1]);

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

        else if( p.nMatch("MOLE") )
        {
                /* molecules, especially for PDR calculations */
                strcpy( save.chSave[save.nsave], "MOLE" );
        }

        else if( p.nMatch("MONI") )
        {
                /* save monitors */
                strcpy( save.chSave[save.nsave], "MONI" );
        }

        else if( p.nMatch("OPTICAL") && p.nMatch("DEPTH") )
        {
                /* check for keyword UNITS on line, then scan wavelength or energy units if present,
                 * units are copied into save.chConPunEnr */
                ChkUnits(p);

                /* "every" option to save this on every zone -
                 * not present then only last zone is saved */
                if( p.nMatch("EVER" ) )
                {
                        /* save every zone */
                        save.lgSaveEveryZone[save.nsave] = true;
                        save.nSaveEveryZone[save.nsave] = 1;
                }
                else
                {
                        /* only save last zone */
                        save.lgSaveEveryZone[save.nsave] = false;
                        save.nSaveEveryZone[save.nsave] = 1;
                }
                
                if( p.nMatch("FINE") )
                {
                        /* save fine continuum optical depths */
                        rfield.lgSaveOpacityFine = true;
                        strcpy( save.chSave[save.nsave], "OPTf" );
                        sprintf( save.chHeader[save.nsave], "#energy/%s\tTau tot\topacity\n",
                                save.chConPunEnr[save.nsave] );
                        /* range option - important since so much data */
                        if( p.nMatch("RANGE") ) 
                        {
                                /* get lower and upper range, eventually must be in Ryd */
                                double Energy1 = p.FFmtRead();
                                double Energy2 = p.FFmtRead();
                                if( p.lgEOL() )
                                {
                                        fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                                if( p.nMatch("UNIT" ) )
                                {
                                        // apply units to range option
                                        const char *energyUnits = p.StandardEnergyUnit();
                                        Energy unitChange;
                                        unitChange.set(Energy1, energyUnits );
                                        Energy1 = unitChange.Ryd();
                                        unitChange.set(Energy2, energyUnits );
                                        Energy2 = unitChange.Ryd();
                                }
                                /* get lower and upper rang in Ryd */
                                save.punarg[save.nsave][0] = (realnum)MIN2( Energy1 , Energy2 );
                                save.punarg[save.nsave][1] = (realnum)MAX2( Energy1 , Energy2 );
                                //fprintf(ioQQQ , "DEBUG units change fine %.3e %.3e\n" , save.punarg[save.nsave][0] ,
                                //              save.punarg[save.nsave][1] );
                                //cdEXIT(EXIT_FAILURE);
                        }
                        else
                        {
                                /* these mean full energy range */
                                save.punarg[save.nsave][0] = 0.;
                                save.punarg[save.nsave][1] = 0.;
                        }
                        /* optional last parameter - how many points to bring together */
                        save.punarg[save.nsave][2] = (realnum)p.FFmtRead();
                        /* default is to bring together ten */
                        if( p.lgEOL() )
                                save.punarg[save.nsave][2] = 10;
                        if( save.punarg[save.nsave][2] < 1 )
                        {
                                fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        /* save coarse continuum optical depths */
                        strcpy( save.chSave[save.nsave], "OPTc" );
                        sprintf( save.chHeader[save.nsave], 
                                "#energy/%s\ttotal\tabsorp\tscat\n",
                                save.chConPunEnr[save.nsave] );
                }

        }
        else if( p.nMatch(" OTS") )
        {
                strcpy( save.chSave[save.nsave], " OTS" );
                sprintf( save.chHeader[save.nsave], 
                        "#otscon, lin, conOpac LinOpc\n" );
        }

        else if( p.nMatch("OVER") && p.nMatch(" OVE") )
        {
                /* save overview of model results */
                strcpy( save.chSave[save.nsave], "OVER" );
                sprintf( save.chHeader[save.nsave], 
                        "#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\tH2O/O\tAV(point)\tAV(extend)\n" );
        }

        else if( p.nMatch(" PDR") )
        {
                strcpy( save.chSave[save.nsave], " PDR" );
                sprintf( save.chHeader[save.nsave], 
                        "#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( p.nMatch("PERF") )
        {
                /* output performance characteristics per zone */
                strcpy( save.chSave[save.nsave], "PERF" );
                sprintf( save.chHeader[save.nsave],
                        "#zone\tdTime\tElapsed t\tnPres2Ioniz\n" );
        }

        else if( p.nMatch("PHYS") )
        {
                /* save physical conditions */
                strcpy( save.chSave[save.nsave], "PHYS" );
                sprintf( save.chHeader[save.nsave], 
                        "#PhyC depth\tTe\tn(H)\tn(e)\tHtot\taccel\tfillfac\n" );
        }

        else if( p.nMatch("POIN") )
        {
                /* save out the pointers */
                save.lgPunPoint = true;
                /* this does not count as a save option (really) */
                strcpy( save.chSave[save.nsave], "" );
                save.lgRealSave[save.nsave] = false;
        }

        else if( p.nMatch("PRES") )
        {
                /* the save pressure command */
                strcpy( save.chSave[save.nsave], "PRES" );
                sprintf( save.chHeader[save.nsave], 
                        "#P depth\tPerror%%\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)\tPgr_Int\tint thin elec\tconv?\n" );
        }

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

        else if( p.nMatch("RECO") )
        {
                if( p.nMatch("COEF") )
                {
                        /* recombination coefficients for everything */

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

                else if( p.nMatch("EFFI") )
                {
                        /* save recombination efficiency */
                        strcpy( save.chSave[save.nsave], "RECE" );
                        sprintf( save.chHeader[save.nsave], 
                                "#Recom effic H, Heo, He+\n" );
                }

                else
                {
                        fprintf( ioQQQ, "No option recognized on this save recombination command\n" );
                        fprintf( ioQQQ, "Valid options are COEFFICIENTS, AGN, and EFFICIENCY\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* save results command, either as single column or wide array */
        else if( p.nMatch("RESU") )
        {
                strcpy( save.chSave[save.nsave], "RESU" );
                if( p.nMatch("COLU") )
                {
                        /* column is key to save single column */
                        strcpy( save.chPunRltType, "column" );
                }
                else
                {
                        /* array is key to save large array */
                        strcpy( save.chPunRltType, "array " );
                }

                /* do not change following, is used as flag in getlines */
                sprintf( save.chHeader[save.nsave], 
                        "#results of calculation\n" );
        }

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

        else if( p.nMatch("SOUR") )
        {

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

                if( p.nMatch("DEPT") )
                {
                        /* print continuum source function as function of depth */
                        strcpy( save.chSave[save.nsave], "SOUD" );
                        sprintf( save.chHeader[save.nsave], 
                                "#continuum source function vs depth\n" );
                }
                else if( p.nMatch("SPEC") )
                {
                        /* print spectrum continuum source function at 1 depth */
                        strcpy( save.chSave[save.nsave], "SOUS" );
                        sprintf( save.chHeader[save.nsave], 
                                "#continuum source function nu/%s\tConEmitLocal/widflx"
                                "\tabs opac\tConSourceFcnLocal\tConSourceFcnLocal/plankf\tConSourceFcnLocal/flux\n",
                                save.chConPunEnr[save.nsave] );
                }
                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);
                }
        }


        /* save spectrum the new form of the save continuum, will eventually replace the standard
         * save continuum command */
        else if( p.nMatch("SPECTRUM") && !p.nMatch("XSPE") )
        {
                /* this flag is checked in PrtComment to generate a caution
                 * if continuum is saved but iterations not performed */
                save.lgPunContinuum = true;

                /* set flag for spectrum */
                strcpy( save.chSave[save.nsave], "CONN" );

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

                sprintf( save.chHeader[save.nsave], 
                        "#Cont Enr/%s\tincid nFn\ttrans\tdiff\tlines \n",
                                save.chConPunEnr[save.nsave] );
        }

        else if( p.nMatch("SPECIAL") )
        {
                /* save special, will call routine SaveSpecial */
                strcpy( save.chSave[save.nsave], "SPEC" );
                sprintf( save.chHeader[save.nsave], "#Special.\n" );
        }

        else if( p.nMatch("SPECIES") )
        {
                strcpy( save.chSave[save.nsave], "SPCS" );

                // option to save information about a particular species,
                // the "second filename" may really be the species label.  Rename here for clarity
                strcpy( chLabel, chSecondFilename );
                if( !lgSecondFilename )
                {
                        strcpy( save.chSaveSpecies[save.nsave], "" );
                }
                else if( strlen(chLabel) >= CHARS_SPECIES )
                {
                        fprintf( ioQQQ,"Species string is limited to %li characters.\nSorry.\n", (long)CHARS_SPECIES );
                        cdEXIT(EXIT_FAILURE);
                }
                else
                        strncpy( save.chSaveSpecies[save.nsave], chLabel, CHARS_SPECIES );

                if (p.nMatch( "COLUMN" ) )
                {
                        /* column densities*/
                        strcpy( save.chSaveArgs[save.nsave], "COLU" );
                }
                else if (p.nMatch( "ENERG" ) )
                {
                        /* energy levels, default Rydbergs but option to change units */
                        ChkUnits(p);
                        strcpy( save.chSaveArgs[save.nsave], "ENER" );
                }
                else if( p.nMatch("LABELS") )
                {
                        strcpy( save.chSaveArgs[save.nsave], "LABE" );
                }
                else if (p.nMatch( "LEVELS" ) )
                {
                        /* the number of levels in this zone */
                        strcpy( save.chSaveArgs[save.nsave], "LEVL" );
                }
                else if (p.nMatch( "POPUL" ) )
                {
                        /* save species population fraction for 1 level*/
                        strcpy( save.chSaveArgs[save.nsave], "POPU" );
                }
                else
                {
                        fprintf( ioQQQ, "ParseSave cannot find a recognized keyword on this SAVE SPECIES command line.\n" );
                        fprintf( ioQQQ, "I know about the keywords COLUMN DENSITIES, LABELS, LEVELS, and POPULATIONS.\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

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

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

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

        else if( p.nMatch("WIND") )
        {
                strcpy( save.chSave[save.nsave], "WIND" );
                sprintf( save.chHeader[save.nsave], 
                        "#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( p.nMatch( "TERM" ) )
                {
                        /* only save for last zone, the terminal velocity, for grids */
                        save.punarg[save.nsave][0] = 0.;
                }
                else
                {
                        /* one means save every zone */
                        save.punarg[save.nsave][0] = 1.;
                }
        }

        else if( p.nMatch("XSPE") )
        {
                /* say that this is a FITS file output */
                save.lgFITS[save.nsave] = true;

                /* the save xspec commands */
                save.lgPunLstIter[save.nsave] = true;

                /* remember that a save xspec command was entered */
                grid.lgSaveXspec = true;

                /* range option - important since so much data */
                if( p.nMatch("RANGE") ) 
                {
                        /* get lower and upper range, must be in keV */
                        save.punarg[save.nsave][0] = (realnum)p.FFmtRead();
                        save.punarg[save.nsave][1] = (realnum)p.FFmtRead();
                        if( p.lgEOL() )
                        {
                                fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in keV.\nSorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                        if( save.punarg[save.nsave][0] >=save.punarg[save.nsave][1] )
                        {
                                fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }

                        grid.LoEnergy_keV = save.punarg[save.nsave][0];
                        grid.HiEnergy_keV = save.punarg[save.nsave][1];
                }
                else
                {
                        /* these mean full energy range */
                        save.punarg[save.nsave][0] = 0;
                        save.punarg[save.nsave][1] = 0;
                }

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

        /* save 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( p.nMatch("COLU") && p.nMatch("DENS") )
        {
                if( p.nMatch("SOME" ))
                {
                        /* flag saying save some column densities */
                        strcpy( save.chSave[save.nsave], "COLS" );
                        parse_save_colden( p, save.chHeader[save.nsave] );
                }
                else
                {
                        /* save column densities table */
                        strcpy( save.chSave[save.nsave], "COLU" );
                }
        }
        else
        {
                fprintf( ioQQQ, 
                        "ParseSave cannot find a recognized keyword on this SAVE command line.\nSorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* only open if file has not already been opened during a previous call */
        if( save.ipPnunit[save.nsave] == NULL )
        {
                string file_name;
                file_name += chFilename;
                string mode = "w";
                if( save.lgFITS[save.nsave] )
                        mode += "b";

                /* open the file with the name and mode generated above */
                save.ipPnunit[save.nsave] = open_data( file_name.c_str(), mode.c_str(), 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( p.nMatch("NO BUFFER") )
                        setbuf( save.ipPnunit[save.nsave] , NULL );
        }

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

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

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

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

        else if( p.nMatch("QHEA") )
        {
                gv.QHSaveFile = save.ipPnunit[save.nsave];
                gv.lgQHPunLast = save.lgPunLstIter[save.nsave];
                save.lgQHSaveFile_noclobber = save.lgNoClobber[save.nsave];
                fprintf( gv.QHSaveFile, 
                        "#Probability distributions from quantum heating routine.\n" );
                save.lgRealSave[save.nsave] = false;
        }

        else if( p.nMatch("POIN") )
        {
                /* save out the pointers */
                save.ipPoint = save.ipPnunit[save.nsave];
                save.lgPunPoint_noclobber = save.lgNoClobber[save.nsave];
                save.lgPunPoint = true;
                fprintf( save.ipPoint, 
                        "#pointers. \n" );
                strcpy( save.chSave[save.nsave], "" );
                save.lgRealSave[save.nsave] = false;
        }

        else if( p.nMatch("RECO") && p.nMatch("COEF") )
        {
                /* recombination coefficients for everything
                 * save.lgioRecom set to false in routine zero, non-zero value
                 * is flag to save 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*/
                save.ioRecom = save.ipPnunit[save.nsave];
                save.lgioRecom_noclobber = save.lgNoClobber[save.nsave];
                /* this is logical flag used in routine ion_recom to create the save output */
                save.lgioRecom = true;
                fprintf( save.ioRecom, 
                        "#recombination coefficients cm3 s-1 for current density and temperature\n" );
                strcpy( save.chSave[save.nsave], "" );
                save.lgRealSave[save.nsave] = false;
        }

        else if( p.nMatch("GRID") )
        {
                /* this enables saving GRID output outside the main SaveDo() loop */
                grid.pnunit = save.ipPnunit[save.nsave];
                save.lgSaveGrid_noclobber = save.lgNoClobber[save.nsave];
        }

        else if( p.nMatch(" MAP") )
        {
                /* say output goes to special save */
                ioMAP = save.ipPnunit[save.nsave];
        }

        /* check that string written into save.chHeader[save.nsave] 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_s(save.chHeader[save.nsave] , '\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 - save.chHeader[save.nsave])>=MAX_HEADER_SIZE-1 )
        {
                fprintf( ioQQQ, "DISASTER save.chHeader[%li] has been overwritten "
                        "with a line too long to be read.\n", save.nsave );
                cdEXIT(EXIT_FAILURE);
        }

        /* if lgPunHeader true and cdHeader has been set to a string then print header
         * logic to prevent more than one header in grid calculation */
        if( save.lgPunHeader[save.nsave] && !nMatch(save.chHeader[save.nsave],save.chNONSENSE) )
        {
                fprintf( save.ipPnunit[save.nsave], "%s", save.chHeader[save.nsave] );
                save.lgPunHeader[save.nsave] = false;
        }

        /* increment total number of save commands, */
        ++save.nsave;
        return;
}

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

        DEBUG_ENTRY( "SaveFilesInit()" );

        ASSERT( lgFIRST );
        lgFIRST = false;

        /* set lgNoClobber to not overwrite files, reset with clobber on save 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++ )
        {
                save.lgNoClobber[i] = lgNoClobberDefault;
        }
        save.lgPunConv_noclobber = lgNoClobberDefault;
        save.lgDROn_noclobber = lgNoClobberDefault;
        save.lgTraceConvergeBase_noclobber = lgNoClobberDefault;
        save.lgPunPoint_noclobber = lgNoClobberDefault;
        save.lgioRecom_noclobber = lgNoClobberDefault;
        save.lgQHSaveFile_noclobber = lgNoClobberDefault;
        save.lgSaveGrid_noclobber = lgNoClobberDefault;

        /* initialize chHeader strings with nonsense, compare later to see if we have any actual headers. */
        save.chNONSENSE = "ArNdY38dZ9us4N4e12SEcuQ";

        for( i=0; i < LIMPUN; i++ )
        {
                save.ipPnunit[i] = NULL;

                // is this a real save command?  set false with the dummy 
                // save commands like save dr
                save.lgRealSave[i] = true;

                // do we need to save header?
                save.lgPunHeader[i] = true;
                strcpy( save.chHeader[i], save.chNONSENSE );
        }

        save.lgTraceConvergeBase = false;

        save.ipDRout = NULL;
        save.lgDROn = false;

        save.ipTraceConvergeBase = NULL;
        save.lgTraceConvergeBase = false;

        save.ipPunConv = NULL;
        save.lgPunConv = false;

        save.ipPoint = NULL;
        save.lgPunPoint = false;

        gv.QHSaveFile = NULL;

        save.ioRecom = NULL;
        save.lgioRecom = false;

        grid.pnunit = NULL;

        ioMAP = NULL;

        return;
}

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

        DEBUG_ENTRY( "CloseSaveFiles()" );

        /* close all save units cloudy opened with save command,
         * lgNoClobber is set false with CLOBBER option on save, says to
         * overwrite the files */
        for( i=0; i < save.nsave; i++ )
        {
                /* if lgFinal is true, we close everything, no matter what. 
                 * this means ignoring "no clobber" options */
                if( save.ipPnunit[i] != NULL && ( !save.lgNoClobber[i] || lgFinal ) )
                {
                        /* Test that any FITS files are the right size! */ 
                        if( save.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(save.ipPnunit[i], 0, SEEK_END);
                                long file_size = ftell(save.ipPnunit[i]);
                                if( file_size%2880 )
                                {
                                        fprintf( ioQQQ, " PROBLEM  FITS file is wrong size!\n" );
                                }
                        }

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

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

        if( save.ipTraceConvergeBase != NULL && ( !save.lgTraceConvergeBase_noclobber || lgFinal ) )
        {
                save.ipTraceConvergeBase = NULL;
                save.lgTraceConvergeBase = false;
        }

        if( save.ipPunConv != NULL && ( !save.lgPunConv_noclobber || lgFinal ) )
        {
                save.ipPunConv = NULL;
                save.lgPunConv = false;
        }
        if( save.ipPoint != NULL && ( !save.lgPunPoint_noclobber || lgFinal ) )
        {
                save.ipPoint = NULL;
                save.lgPunPoint = false;
        }
        if( gv.QHSaveFile != NULL  && ( !save.lgQHSaveFile_noclobber || lgFinal ) )
        {
                gv.QHSaveFile = NULL;
        }
        if( save.ioRecom != NULL  && ( !save.lgioRecom_noclobber || lgFinal ) )
        {
                save.ioRecom = NULL;
                save.lgioRecom = false;
        }
        if( grid.pnunit != NULL && ( !save.lgSaveGrid_noclobber || lgFinal ) )
        {
                grid.pnunit = NULL;
        }
        ioMAP = NULL;

        return;
}

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

        DEBUG_ENTRY( "ChkUnits()" );

        /* option to set units for continuum energy in save output */
        if( p.nMatch("UNITS") )
        {
                // p.StandardEnergyUnit() will terminate if no unit was recognized
                save.chConPunEnr[save.nsave] = p.StandardEnergyUnit();
        }
        else
        {
                save.chConPunEnr[save.nsave] = StandardEnergyUnit(" RYD ");
        }
        return;
}

---