parse_set.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 */
/*ParseSet scan parameters off SET command */
#include "cddefines.h"
#include "physconst.h"
#include "geometry.h"
#include "input.h"
#include "prt.h"
#include "atomfeii.h"
#include "mole.h"
#include "rt.h"
#include "phycon.h"
#include "optimize.h"
#include "hcmap.h"
#include "hmi.h"
#include "dense.h"
#include "h2.h"
#include "co.h"
#include "iterations.h"
#include "conv.h"
#include "secondaries.h"
#include "rfield.h"
#include "ionbal.h"
#include "numderiv.h"
#include "dynamics.h"
#include "iso.h"
#include "mole.h"
#include "predcont.h"
#include "save.h"
#include "stopcalc.h"
#include "opacity.h"
#include "hydrogenic.h"
#include "peimbt.h"
#include "radius.h"
#include "atmdat.h"
#include "continuum.h"
#include "grains.h"
#include "grainvar.h"
#include "parser.h"
#include "lines.h"
#include "monitor_results.h"
#include "thirdparty.h"
 
void ParseSet(Parser &p)
{
        long int ip;
        char chString_quotes_lowercase[INPUT_LINE_LENGTH];
        DEBUG_ENTRY( "ParseSet()" );
 
        /* first check for any strings within quotes - this will get the string
         * and blank it out, so that these are not confused with keywords.  if
         * double quotes not present routine returns unity, zero if found*/
        bool lgQuotesFound = true;
        if (p.GetQuote(chString_quotes_lowercase, false))
                lgQuotesFound = false;
 
        /* commands to set certain variables, "SET XXX=" */
        if (p.nMatch("ASSE") && p.nMatch("ABOR"))
        {
                /* set assert abort command, to tell the code to raise SIGABRT when an
                 * assert is thrown - this can be caught by the debugger. */
                cpu.i().setAssertAbort(true);
        }
 
        else if (p.nMatch("MONI") && p.nMatch("SCIE"))
        {
                /* print monitored values using scientific notation.  Useful when an asserted value is 
                 * less than 10^-4 of the normalization line. */
                lgPrtSciNot = true;
        }
 
        else if (p.nMatch("ATOM") && p.nMatch("DATA"))
        {
                /* set atomic data */
                long int ion;
                bool UNUSED lgAs = false;
                bool lgCS = false;
                bool UNUSED lgDR = false;
                const char *chMole;
 
                /* As? */
                if (p.nMatch(" AS "))
                        lgAs = true;
                /* DR - dielectronic recombination */
                else if (p.nMatch(" DR "))
                        lgDR = true;
                /* lgCS collision strengths */
                else if (p.nMatch(" CS "))
                        lgCS = true;
                else
                {
                        /* no keyword */
                        fprintf(ioQQQ,
                                        " One of the keywords AS DR or CS must appear to change"
                                                " the transition probabilities, dielectronic recombination rate,"
                                                " or collision strength.\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
 
                // only molecules enter at present
                if (p.nMatch(" H2 "))
                {
                        /* H2 */
                        chMole = "H2";
                }
                else
                {
                        /* nothing recognized */
                        fprintf(ioQQQ, " No molecule was on this SET ATOMIC DATA command.\n Sorry.\n");
                        fprintf(ioQQQ, " Use SET ATOMIC DATA ION to change an ion.\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
 
                if (strcmp(chMole, "H2") == 0)
                {
                        if (p.nMatch(" H ") && lgCS)
                        {
                                /* change the H - H2 collision data set */
                                ion = (long int) p.FFmtRead();
                                if (ion != 2)
                                        TotalInsanity();
                                long int nYear = (long) p.FFmtRead();
                                if (p.lgEOL())
                                        p.NoNumb("collison data set (year)");
                                if (nYear == 1999)
                                {
                                        /* use the coefficients from
                                         *>>refer       H2      H collision     Le Bourlot, J., Pineau des Forets,
                                         *>>refercon    G., & Flower, D.R. 1999, MNRAS, 305, 802 */
                                        h2.lgH2_H_coll_07 = false;
                                        h2.coll_source[0].filename = "coll_rates_H_99.dat";
                                }
                                else if (nYear == 2007)
                                {
                                        /* use the coefficients from
                                         *>>refer       H2      H collision     Wrathmall, S. A., Gusdorf A.,
                                         *>>refercon    & Flower, D.R. 2007, MNRAS, 382, 133 */
                                        h2.lgH2_H_coll_07 = true;
                                        h2.coll_source[0].filename = "coll_rates_H_07.dat";
                                }
                                else 
                                {
                                        /* not an option */
                                        fprintf(ioQQQ, " the SET ATOMIC DATA MOLECULE H2"
                                                " H CS command must have year 1999 or 2007.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else if (p.nMatch(" HE ") && lgCS)
                        {
                                /* change the He - H2 collision data set */
                                if (p.nMatch("ORNL"))
                                {
                                        /* use the coefficients from
                                         *>>refer       H2      He collision    Lee et al in preparation */
                                        h2.lgH2_He_ORNL = true;
                                        h2.coll_source[1].filename = "coll_rates_He_ORNL.dat";
                                }
                                else if( p.nMatch( "BOUR") )
                                {
                                        /* use the coefficients from
                                        *>>refer        H2      He collision    Le Bourlot, J., Pineau des Forets, 
                                        *>>refercon     G., & Flower, D.R. 1999, MNRAS, 305, 802*/
                                        h2.lgH2_He_ORNL = false;
                                        h2.coll_source[1].filename = "coll_rates_He_LeBourlot.dat";
                                }
                                else
                                {
                                        /* not an option */
                                        fprintf(ioQQQ, " the SET ATOMIC DATA MOLECULE H2"
                                                "He CS command must have ORNL or Le BOURlot.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else
                        {
                                fprintf( ioQQQ, " the SET ATOMIC DATA H2 CS command requires a collider keyword.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }
 
        /* check energy conservation on every zone - relatively slow */
        else if (p.nMatch("CHECK") && p.nMatch("ENERGY") && p.nMatch("EVERY")
                        && p.nMatch("ZONE"))
        {
                continuum.lgCheckEnergyEveryZone = true;
        }
 
        else if (p.nMatch("UPDA") && p.nMatch("COUP") &&
                                p.nMatch("EVER") && p.nMatch("ION") )
        {
                conv.lgUpdateCouplings = true;
        }
 
 
        /* enable Dere07 collisional ionization rate coeffs */
        else if (p.nMatch(" COLL") && p.nMatch(" IONIZ"))
        {
 
                if (p.nMatch(" HYBRID"))
                {
                        atmdat.CIRCData = t_atmdat::HYBRID;
                }
                else if (p.nMatch(" DIMA"))
                {
                        atmdat.CIRCData = t_atmdat::DIMA;
                }
                else
                {
                        fprintf(ioQQQ,
                                        " \nPROBLEM Unrecognized set collisional ionization data option.\n");
                        fprintf(ioQQQ, " Valid options are Dima or Hybrid.\n");
                        fprintf(ioQQQ, " See Hazy 1 for details.\n");
                        cdEXIT( EXIT_FAILURE );
                }
        }
 
        else if( p.nMatch(" H2 " ) && p.nMatch( "CONT" ) && p.nMatch( "DISS" ))
        {
                if( p.nMatch( "STAN" ) )
                {
                        /* This uses on the H2 direct photodissociation cross sections
                         * calculated by Philip Stancil */
                        mole_global.lgStancil = true;
                }
 
                else if( p.nMatch( "AD69" ) )
                {
                         /* Use constant cross section from
                            >>refer     H2      dissoc  Allison, A.C. & Dalgarno, A. 1969,   
                            Atomic Data, 1, 91 */ 
 
                        mole_global.lgStancil = false;
                }
                else
                {
                        /* should not have happened ... */
                        fprintf( ioQQQ, " There should have been an option on this SET H2 CONTinuum DISSociation command.\n" );
                        fprintf( ioQQQ, " consult Hazy to find valid options.\n Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);                   
                }
 
        }
 
        else if( p.nMatch(" CHA") && !p.nMatch( "HO ") && !p.nMatch( "HHE") )
        {
                /* set log of minimum charge transfer rate for high ions and H
                 * default of 1.92e-10 set in zero */
                atmdat.HCTAlex = p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("minimum charge transfer rate");
                }
                if (atmdat.HCTAlex < 0.)
                {
                        atmdat.HCTAlex = pow(10., atmdat.HCTAlex);
                }
        }
        else if( p.nMatch("CHEM") && !p.nMatch( "HO ") && !p.nMatch(" HHE") )
        {
                /* turn on Steve Federman's chemistry */
                if (p.nMatch("FEDE"))
                {
                        if (p.nMatch(" ON "))
                        {
                                /* This turns on the diffuse cloud chemical rates of 
                                 * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
                                mole_global.lgFederman = true;
                        }
                        else if( p.nMatch( " OFF" ) )
                        {
                                mole_global.lgFederman = false;
                        }
                        else
                        {
                                /* this is the default when command used - true */
                                mole_global.lgFederman = true;
                        }
                }
                /* >>chng 06 may 30 --NPA.  Turn on non-equilibrium chemistry */
                else if (p.nMatch(" NON") && p.nMatch("EQUI"))
                {
 
                        /* option to use effective temperature as defined in
                         * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
                         * By default, this is false - changed with set chemistry command */
 
                        mole_global.lgNonEquilChem = true;
 
                        /* >>chng 06 jul 21 -- NPA.  Option to include non-equilibrium
                         * effects for neutral reactions with a temperature dependent rate.
                         * Reasoning is that non-equilibrium chemistry is caused by MHD, and
                         * if magnetic field is only coupled to ions, then neutrals may not be 
                         * affected.  
                         * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
                         * By default, this is true - changed with set chemistry command */
 
                        if (p.nMatch("NEUT"))
                        {
                                if (p.nMatch(" ON "))
                                {
                                        /* This turns on the diffuse cloud chemical rates of 
                                         * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
                                        mole_global.lgNeutrals = true;
                                }
                                else if( p.nMatch( " OFF" ) )
                                {
                                        mole_global.lgNeutrals = false;
                                }
                                else
                                {
                                        /* this is the default when command used - true */
                                        mole_global.lgNeutrals = true;
                                }
                        }
                }
 
                /* turn off proton elimination rates, which are of the form
                 *
                 *
                 *    A + BH+ -->  AB + H+ or
                 *    AH + B+ -->  AB + H+ 
                 *
                 * the following paper:
                 *
                 * >>refer      CO      chemistry       Huntress, W. T., 1977, ApJS, 33, 495
                 * says reactions of these types are much less likely than
                 * identical reactions which leave the product AB ionized (AB+), 
                 * leaving an H instead of H+ (this is called H atom elimination 
                 * currently we only have one reaction of this type, it is
                 * C+ + OH -> CO + H+ */
                else if (p.nMatch("PROT") && p.nMatch("ELIM"))
                {
                        if (p.nMatch(" ON "))
                        {
                                /* This turns on the diffuse cloud chemical rates of 
                                 * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
                                mole_global.lgProtElim = true;
                        }
                        else if( p.nMatch( " OFF" ) )
                        {
                                mole_global.lgProtElim = false;
                        }
                        else
                        {
                                /* this is the default when command used - true */
                                mole_global.lgProtElim = true;
                        }
                }
 
                else
                {
                        /* should not have happened ... */
                        fprintf(ioQQQ,
                                        " There should have been an option on this SET CHEMISTRY command.\n");
                        fprintf(ioQQQ, " consult Hazy to find valid options.\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        /* set collision strength averaging ON / OFF */
        else if( p.nMatch("COLL") && p.nMatch("STRE") && p.nMatch("AVER") )
        {
                if( p.nMatch(" OFF") )
                {
                         iso_ctrl.lgCollStrenThermAver = false;
                }
                else
                {
                        /* this is default behavior of this command */
                         iso_ctrl.lgCollStrenThermAver = true;
                }
        }
 
        else if (p.nMatch("COVE"))
        {
                iterations.lgConverge_set = true;
                /* set coverage - limit number of iterations and zones */
                if (p.nMatch("FAST"))
                {
                        iterations.lim_zone = 1;
                        iterations.lim_iter = 0;
                }
                else
                {
                        iterations.lim_zone = 10;
                        iterations.lim_iter = 1;
                }
        }
 
        else if (p.nMatch("CSUP"))
        {
                /* force H^0 secondary ionization rate, log entered */
                secondaries.SetCsupra = (realnum) p.FFmtRead();
                secondaries.lgCSetOn = true;
                if (p.lgEOL())
                {
                        p.NoNumb("secondary ionization rate");
                }
                secondaries.SetCsupra = (realnum) pow((realnum) 10.f,
                                secondaries.SetCsupra);
        }
 
        else if (p.nMatch(" D/H"))
        {
                /* set deuterium abundance, D to H ratio */
                hydro.D2H_ratio = p.FFmtRead();
                if (hydro.D2H_ratio <= 0. || p.nMatch(" LOG"))
                {
                        hydro.D2H_ratio = pow(10., hydro.D2H_ratio);
                }
                if (p.lgEOL())
                {
                        p.NoNumb("D to H ratio");
                }
        }
 
        else if( p.nMatch("DENS") && p.nMatch("TOLE") )
        {
                /* set error in total gas-phase density of each element, including molecules */
                conv.GasPhaseAbundErrorAllowed = (realnum)p.FFmtRead();
                if( p.lgEOL() )
                {
                        p.NoNumb("density tolerance");
                }
                if( conv.GasPhaseAbundErrorAllowed <= 0. )
                {
                        conv.GasPhaseAbundErrorAllowed = (realnum)pow((realnum)10.f,conv.GasPhaseAbundErrorAllowed);
                }
        }
 
        else if( p.nMatch( " HO " ) && p.nMatch( "CHAR" ) )
        {
                fprintf(ioQQQ, " The SET HO CHAR command is no longer supported.\n" );
                cdEXIT(EXIT_FAILURE);
        }
 
        else if( p.nMatch( " HHE" ) && p.nMatch( "CHAR" ) )
        {
                fprintf(ioQQQ, " The SET HHE CHAR command is no longer supported.\n" );
                cdEXIT(EXIT_FAILURE);
        }
 
        else if( p.nMatch("12C1") )
        {
                /* set the 12C to 13C abundance ratio - default is 30 */
                /* first two numbers on the line are 12 and 13 - we don't want them */
                co.C12_C13_isotope_ratio_parsed = (realnum) p.FFmtRead();
                co.C12_C13_isotope_ratio_parsed = (realnum) p.FFmtRead();
 
                /* now we can get the ratio */
                co.C12_C13_isotope_ratio_parsed = (realnum) p.FFmtRead();
                if (p.lgEOL())
                        p.NoNumb("12C to 13C abundance ratio");
 
                if (co.C12_C13_isotope_ratio_parsed <= 0. || p.nMatch(" LOG"))
                        co.C12_C13_isotope_ratio_parsed = (realnum) pow((realnum) 10.f,
                                        co.C12_C13_isotope_ratio_parsed);
        }
 
        /* set dynamics ... */
        else if (p.nMatch("DYNA"))
        {
                /* set dynamics advection length */
                if (p.nMatch("ADVE") && p.nMatch("LENG"))
                {
                        /* <0 => relative fraction of length, + val in cm */
                        dynamics.AdvecLengthInit = p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("advection length");
                        /* if fraction is present, then number was linear fraction, if not present
                         * then physical length in cm, log10 */
                        if (p.nMatch("FRAC"))
                        {
                                /* we scanned in the number, if it is a negative then it is the log of the fraction */
                                if (dynamics.AdvecLengthInit <= 0.)
                                        dynamics.AdvecLengthInit = pow(10.,
                                                        dynamics.AdvecLengthInit);
 
                                /* make neg sign as flag for this in dynamics.c */
                                dynamics.AdvecLengthInit *= -1.;
                        }
                        else
                        {
                                /* fraction did not occur, the number is the log of the length in cm -
                                 * convert to linear cm */
                                dynamics.AdvecLengthInit = pow(10., dynamics.AdvecLengthInit);
                        }
                }
                else if (p.nMatch("PRES") && p.nMatch("MODE"))
                {
                        dynamics.lgSetPresMode = true;
                        if (p.nMatch("SUBS"))
                        {
                                /* subsonic */
                                strcpy(dynamics.chPresMode, "subsonic");
                        }
                        else if (p.nMatch("SUPE"))
                        {
                                /* supersonic */
                                strcpy(dynamics.chPresMode, "supersonic");
                        }
                        else if (p.nMatch("STRO"))
                        {
                                /* strong d */
                                strcpy(dynamics.chPresMode, "strongd");
                        }
                        else if (p.nMatch("ORIG"))
                        {
                                /* original */
                                strcpy(dynamics.chPresMode, "original");
                        }
                }
                else if (p.nMatch("ANTI") && p.nMatch("DEPT"))
                {
                        dynamics.lgSetPresMode = true;
                        strcpy(dynamics.chPresMode, "antishock");
                        /* shock depth */
                        /* get log of shock depth in cm */
                        dynamics.ShockDepth = p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("antishock depth");
                        dynamics.ShockDepth = pow(10., dynamics.ShockDepth);
                }
                else if (p.nMatch("ANTI") && p.nMatch("MACH"))
                {
                        dynamics.lgSetPresMode = true;
                        strcpy(dynamics.chPresMode, "antishock-by-mach");
                        /* Mach number */
                        /* get (isothermal) Mach number where we want antishock to occur */
                        dynamics.ShockMach = p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("antishock-by-mach");
                }
                else if (p.nMatch("RELA"))
                {
                        /* set how many iterations we will start with, before allowing
                         * changes.  This allows the solution to relax to an equilibrium */
                        dynamics.n_initial_relax = (long int) p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("relaxation cycles before start of dynamics");
                        else if (dynamics.n_initial_relax < 2)
                        {
                                fprintf(ioQQQ,
                                                " First iteration to relax dynamics must be > 1."
                                                        "It was %li. Sorry.\n", dynamics.n_initial_relax);
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else if (p.nMatch("SHOC") && p.nMatch("DEPT"))
                {
                        dynamics.lgSetPresMode = true;
                        strcpy(dynamics.chPresMode, "shock");
                        /* shock depth */
                        /* get log of shock depth in cm */
                        dynamics.ShockDepth = p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("shock depth");
                        dynamics.ShockDepth = pow(10., dynamics.ShockDepth);
                }
                else if (p.nMatch("POPU") && p.nMatch("EQUI"))
                {
                        // set dynamics populations eqauilibrium
                        // force equilibrium populations 
                        dynamics.lgEquilibrium = true;
                }
                else
                {
                        /* should not have happened ... */
                        fprintf(ioQQQ,
                                        " There should have been an option on this SET DYNAMICS command.\n");
                        fprintf(ioQQQ, " consult Hazy to find valid options.\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("DIDZ"))
        {
                /* set parameter to do with choice of dr;
                 * par is the largest optical depth to allow in the zone. 
                 * >>chng 96 jan 08 had been two numbers - dtau1 removed */
                radius.drChange = (realnum) p.FFmtRead();
                if (radius.drChange <= 0.)
                {
                        radius.drChange = (realnum) pow((realnum) 10.f, radius.drChange);
                }
                if (p.lgEOL())
                {
                        p.NoNumb("largest optical depth allowed in zone");
                }
        }
 
        /* something to do with electron density */
        else if (p.nMatch("EDEN"))
        {
                /* set eden convergence sets convergence criterion, keyword parallel with
                 * set temperature convergence,  but also accept (older) error as keyword */
                if (p.nMatch("CONV") || p.nMatch("ERRO"))
                {
                        /* keyword is eden convergence  
                         * set tolerance in eden match */
                        conv.EdenErrorAllowed = p.FFmtRead();
                        if (p.lgEOL())
                        {
                                p.NoNumb("electron density error allowed");
                        }
 
                        if (conv.EdenErrorAllowed < 0.)
                        {
                                conv.EdenErrorAllowed = pow(10., conv.EdenErrorAllowed);
                        }
                }
                else if(p.nMatch("FRACTION"))
                {
                        /* set eden fraction sets the ratio eden/hden */
                        dense.EdenFraction = p.FFmtRead();
                        if (p.lgEOL())
                        {
                                p.NoNumb("electron density fraction");
                        }
 
                        if (dense.EdenFraction <= 0. )
                        {
                                dense.EdenFraction = pow(10., dense.EdenFraction);
                        }
                        /* warn that this model is meaningless */
                        phycon.lgPhysOK = false;
                }
                else
                {
                        /* no keyword, assume log of electron density */
                        /* set the electron density */
                        dense.EdenSet = (realnum) pow(10., p.FFmtRead());
                        if (p.lgEOL())
                        {
                                p.NoNumb("electron density");
                        }
 
                        /* warn that this model is meaningless */
                        phycon.lgPhysOK = false;
                }
        }
 
        /* Stop using gbar to fill in dBase transitions  */
        else if( p.nMatch("GBAR"))
        {
                if( p.nMatch(" OFF ") )
                {
                        atmdat.lgGbarOn = false;
                }
        }
 
        /* something to do with electron density */
        else if (p.nMatch("IONI") && p.nMatch("TOLE"))
        {
                /* keyword is ionization tolerance */
                conv.IonizErrorAllowed = p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("ionization error allowed");
                }
                
                if (conv.IonizErrorAllowed < 0.)
                {
                        conv.IonizErrorAllowed = pow(10., (double)conv.IonizErrorAllowed);
                }
        }
 
        // enable isotopes and isotopologues
        else if (p.nMatch("ISOTOPE") || p.nMatch("ISOTOPO") )
        {
                if( p.nMatch(" ALL ") )
                {
                        // leave out Deuterium for now, to test scaled isotopologues more cleanly
                        for( long nelem = ipHELIUM ; nelem < LIMELM; ++nelem )
                        {
                                mole_global.lgTreatIsotopes[nelem] = true;
                        }
                }
        }
 
        else if (p.nMatch("FINE") && p.nMatch("CONT"))
        {
                /* set fine continuum resolution - an element name, used to get
                 * thermal width, and how many resolution elements to use to resolve 
                 * a line of this element at 1e4 K 
                 * first get an element name, nelem is atomic number on C scale
                 * default is iron */
                if ((rfield.fine_opac_nelem = p.GetElem()) < 0)
                {
                        fprintf(ioQQQ,
                                        " An element name must appear on this line\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
                /* set the number of resolution elements in HWHM at 1e4 K for turbulent 
                 * velocity field, default is one element */
                rfield.fine_opac_nresolv = (long int) p.FFmtRead();
                if (rfield.fine_opac_nresolv < 1)
                {
                        fprintf(ioQQQ,
                                        " The number of resolution elements within FWHM of line must appear\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
 
                /* option to specify the lower and upper limits of the fine continuum
                 * the upper limit is always optional.  */
                if (!p.lgEOL())
                {
                        realnum lower_limit = (realnum) p.FFmtRead();
                        if (lower_limit < 0)
                                lower_limit = pow((realnum) 10.f, lower_limit);
 
                        if (lower_limit > 0.2f)
                                fprintf(
                                                ioQQQ,
                                                " The fine continuum lower limit is quite high (%f Ryd). Please check.\n",
                                                lower_limit);
 
                        /* reset to coarse limits if arguments are beyond them */
                        rfield.fine_ener_lo = MAX2(rfield.emm, lower_limit);
 
                        if (!p.lgEOL())
                        {
                                realnum upper_limit = (realnum) p.FFmtRead();
                                if (upper_limit < 0)
                                        upper_limit = pow((realnum) 10.f, upper_limit);
 
                                if (upper_limit < 10.f)
                                        fprintf(
                                                        ioQQQ,
                                                        " The fine continuum upper limit is quite low (%f Ryd). Please check.\n",
                                                        upper_limit);
 
                                /* reset to coarse limits if arguments are beyond them */
                                rfield.fine_ener_hi = MIN2(rfield.egamry, upper_limit);
                        }
                }
        }
 
        /* set grain command - but not set H2 grain command */
        else if (p.nMatch("GRAI") && !p.nMatch(" H2 "))
        {
                if (p.nMatch("HEAT"))
                {
                        /* scale factor to change grain heating as per Allers et al. */
                        gv.GrainHeatScaleFactor = (realnum) p.FFmtRead();
                        /* warn that this model is not the best we can do */
                        phycon.lgPhysOK = false;
                        if (p.lgEOL())
                        {
                                p.NoNumb("grain heating");
                        }
                }
                else
                {
                        fprintf(ioQQQ,
                                        " A keyword must appear on the SET GRAIN line - options are HEAT \n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        /* these are the "set Leiden hack" commands, used to turn off physics and
         * sometimes replace with simple approximation */
        else if (p.nMatch("LEID") && p.nMatch("HACK"))
        {
                if (p.nMatch("H2* ") && p.nMatch(" OFF"))
                {
                        /* turn off reactions with H2* in the chemistry network */
                        hmi.lgLeiden_Keep_ipMH2s = false;
                        /* warn that this model is not the best we can do */
                        phycon.lgPhysOK = false;
                }
                else if (p.nMatch("CR ") && p.nMatch(" OFF"))
                {
                        /* the CR Leiden hack option - turn off CR excitations of H2 */
                        hmi.lgLeidenCRHack = false;
                }
                else if( p.nMatch("RATE") &&  p.nMatch("UMIS"))
                {
                        /* This command will use the rates given in the UMIST database,
                         * It will set to zero many reactions that are not in UMIST */
                        mole_global.lgLeidenHack = true;
                }
        }
 
        /* set H2 ... */
        else if (p.nMatch(" H2 "))
        {
                ip = (long int) p.FFmtRead();
                if (ip != 2)
                {
                        fprintf(ioQQQ,
                                        " The first number on this line must be the 2 in H2\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }
 
                /* SET H2 SMALL MODEL or SOLOMON - which approximation to Solomon process */
                if (p.nMatch("SOLO") || (p.nMatch("SMAL") && p.nMatch("MODE")))
                {
                        if (p.nMatch("SOLO"))
                        {
                                /* warn that Solomon will not be used forever */
                                fprintf(ioQQQ,
                                                "PROBLEM - *set H2 Solomon* has been changed to *set H2 small model*."
                                                        "  This is OK for now but it may not work in a future version.\n");
                        }
                        if (p.nMatch("TH85"))
                        {
                                /* rate from is eqn A8 of Tielens and Hollenbach 85a, */
                                hmi.chH2_small_model_type = 'T';
                        }
                        else if (p.nMatch(" BHT"))
                        {
                                /* the improved H2 formalism given by 
                                 *>>refer       H2      dissoc  Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M
                                 >>refcon       1990, ApJ, 365, 620 */
                                hmi.chH2_small_model_type = 'H';
                        }
                        else if (p.nMatch("BD96"))
                        {
                                /* this rate is equation 23 of
                                 *>>refer       H2      dissoc  Bertoldi, F., & Draine, B.T., 1996, 458, 222 */
                                /* this is the default */
                                hmi.chH2_small_model_type = 'B';
                        }
                        else if (p.nMatch("ELWE"))
                        {
                                /* this rate is equation 23 of
                                 *>>refer       H2      dissoc  Elwert et al., in preparation */
                                /* this is the default */
                                hmi.chH2_small_model_type = 'E';
                        }
                        else
                        {
                                fprintf(ioQQQ,
                                                " One of the keywords TH85, _BHT, BD96 or ELWErt must appear.\n Sorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
 
                /* series of commands that deal with grains */
                /* which approximation to grain formation pumping */
                if (p.nMatch("GRAI") && p.nMatch("FORM") && p.nMatch("PUMP"))
                {
                        if (p.nMatch("DB96"))
                        {
                                /* Draine & Bertoldi 1996 */
                                hmi.chGrainFormPump = 'D';
                        }
                        else if (p.nMatch("TAKA"))
                        {
                                /* the default from
                                 * >>refer      H2      pump    Takahashi, Junko, 2001, ApJ, 561, 254-263 */
                                hmi.chGrainFormPump = 'T';
                        }
                        else if (p.nMatch("THER"))
                        {
                                /* thermal distribution, upper right column of page 239 of
                                 *>>refer       H2      formation       Le Bourlot, J, 1991, A&A, 242, 235 */
                                hmi.chGrainFormPump = 't';
                        }
                        else
                        {
                                fprintf(ioQQQ,
                                                " The grain form pump option is wrong.\n Sorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
 
                /* which approximation to Jura rate */
                else if (p.nMatch("JURA"))
                {
                        if (p.nMatch("TH85"))
                        {
                                /* rate from is eqn A8 of Tielens and Hollenbach 85a*/
                                hmi.chJura = 'T';
                        }
                        else if (p.nMatch("CT02"))
                        {
                                /* this rate is equation Cazeux & Tielens */
                                hmi.chJura = 'C';
                        }
                        else if (p.nMatch("SN99"))
                        {
                                /* this rate is from Sternberg & Neufeld 99 */
                                hmi.chJura = 'S';
                        }
                        else if (p.nMatch("RATE"))
                        {
                                /* set H2 rate - enters log of Jura rate - F for fixed
                                 * no dependence on grain properties 
                                 * had been C, a bug since triggered Cazeux & Tielens
                                 * >>chng 07 jan 10, bug caught by Robin Williams & Fixed by Nick Abel */
                                hmi.chJura = 'F';
                                hmi.rate_h2_form_grains_set = pow(10., p.FFmtRead());
                                if (p.lgEOL())
                                {
                                        /* no number on the line so use Jura's value, 3e-17 
                                         * >>refer      H2      Jura    Jura, M. 1975, ApJ, 197, 575 */
                                        hmi.rate_h2_form_grains_set = 3e-17;
                                }
                        }
                        else if (p.nMatch("SCAL"))
                        {
                                /* this is a scale factor to multiply the Jura rate */
                                hmi.ScaleJura = (realnum) p.FFmtRead();
                                /* log or negative number means log was entered */
                                if (p.nMatch(" LOG") || hmi.ScaleJura < 0.)
                                {
                                        hmi.ScaleJura = (realnum) pow(10., (double) hmi.ScaleJura);
                                }
                                if (p.lgEOL())
                                        p.NoNumb("scale for Jura rate");
 
                                /* option to vary scale factor */
                                if (optimize.lgVarOn)
                                {
                                        optimize.nvarxt[optimize.nparm] = 1;
                                        strcpy(optimize.chVarFmt[optimize.nparm],
                                                        "SET H2 JURA SCALE %f LOG");
 
                                        /* pointer to where to write */
                                        optimize.nvfpnt[optimize.nparm] = input.nRead;
 
                                        /* log of Jura rate scale factor will be parameter */
                                        optimize.vparm[0][optimize.nparm] = (realnum) log10(
                                                        hmi.ScaleJura);
                                        optimize.vincr[optimize.nparm] = 0.3f;
 
                                        ++optimize.nparm;
                                }
                        }
                        else
                        {
                                fprintf(ioQQQ, " The Jura rate option is wrong.\n Sorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
 
                /* what temperature to use for binding energy, Tad in Le Bourlot, J., 2000, A&A, 360, 656-662  */
                else if (p.nMatch(" TAD"))
                {
                        hmi.Tad = (realnum) p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("temperature for binding energy");
                        /* log if <= 10. unless linear key appears too */
                        if (hmi.Tad <= 10. && !p.nMatch("LINE"))
                                hmi.Tad = (realnum) pow((realnum) 10.f, hmi.Tad);
                }
 
                else if (p.nMatch("FRAC"))
                {
                        /* this is special option to force H2 abundance to value for testing
                         * this factor will multiply the hydrogen density to become the H2 density
                         * no attempt to conserve particles, or do the rest of the molecular equilibrium
                         * set consistently is made */
                        hmi.H2_frac_abund_set = p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("H2 fractional abundance");
 
                        /* a number <= 0 is the log of the ratio */
                        if (hmi.H2_frac_abund_set <= 0.)
                                hmi.H2_frac_abund_set = pow(10., hmi.H2_frac_abund_set);
                        /* don't let it exceed 0.5 */
                        /* >>chng 03 jul 19, from 0.5 to 0.4999, do not want atomic density exactly zero */
                        hmi.H2_frac_abund_set = MIN2(0.49999, hmi.H2_frac_abund_set);
                }
#if 0
                else if( p.nMatch("FORM") && p.nMatch("SCAL") )
                {
                        /* this is special option to scale H2 formation rate. 
                         * In the fully molecular or fully ionized limits,
                         * this should supercede the above "FRAC" option because
                         * it allows the same thing without breaking the chemistry
                         * or any conservation checks */
                        hmi.H2_formation_scale = p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("H2 formation scale");
 
                        /* a number <= 0 is the log of the ratio */
                        if (hmi.H2_formation_scale <= 0.)
                                hmi.H2_formation_scale = pow(10., hmi.H2_frac_abund_set);
                }
#endif
        }
 
        /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the
         * electron density to account for collisions with atomic H.  it is an order of
         * magnitude guess, so this command provides ability to see whether it affects results */
        else if (p.nMatch("HCOR"))
        {
                dense.HCorrFac = (realnum) p.FFmtRead();
                if (p.lgEOL())
                        p.NoNumb("scale for H0 correction to e- collision rate");
        }
 
        else if (p.nMatch(" PAH"))
        {
                /* set one of several possible PAH abundance distribution functions */
                if (lgQuotesFound)
                {
                        /* abundance depends specified species as fraction of Htot */
                        gv.chPAH_abundance = chString_quotes_lowercase;
                }
                else if (p.nMatch("CONS"))
                {
                        /* constant PAH abundance */
                        gv.chPAH_abundance = "CON";
                }
                else if (p.nMatch("BAKE"))
                {
                        /* turn on simple PAH heating from Bakes & Tielens - this is very approximate */
                        /*>>>refer      PAH     heating Bakes, E.L.O., & Tielens, A.G.G.M. 1994, ApJ, 427, 822 */
                        gv.lgBakesPAH_heat = true;
                        /* warn that this model is not the best we can do: energy conservation is violated */
                        phycon.lgPhysOK = false;
                }
                else
                {
                        fprintf(ioQQQ,
                                        " a string, or one of the keywords BAKES, or CONStant must appear.\n Sorry.");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        /* something to do with pressure */
        else if (p.nMatch("PRES"))
        {
                /* tolerance on pressure convergence */
                if (p.nMatch("CONV"))
                {
                        /* keyword is tolerance 
                         * set tolerance or relative error in pressure match */
                        conv.PressureErrorAllowed = (realnum) p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("pressure convergence tolerance");
 
                        if (conv.PressureErrorAllowed < 0.)
                                conv.PressureErrorAllowed = (realnum) pow((realnum) 10.f,
                                                conv.PressureErrorAllowed);
                }
                else if( p.nMatch("IONI") )
                {
                        /* set limit to number of calls from pressure to ionize solver,
                         * this set limit to conv.nPres2Ioniz */
                        conv.limPres2Ioniz = (long) p.FFmtRead();
                        if (p.lgEOL())
                        {
                                p.NoNumb("number of calls from pressure to ion solver");
                        }
                        else if (conv.limPres2Ioniz <= 0)
                        {
                                fprintf(ioQQQ, " The limit must be greater than zero.\n Sorry.");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
 
                else
                {
                        /* >>chng 04 mar 02, printout had been wrong, saying TOLErange
                         * rather than CONVergence.  Nick Abel */
                        fprintf(ioQQQ,
                                        " I didn\'t recognize a key on this SET PRESSURE line.\n");
                        fprintf(ioQQQ, " The ones I know about are: CONVergence and IONIze.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else if (p.nMatch("RECOMBIN"))
        {
                /* dielectronic recombination */
                if (p.nMatch("DIELECTR"))
                {
                        if (p.nMatch("MEAN"))
                        {
                                /* change various factors for the dielectronic recombination means */
                                if (p.nMatch(" OFF"))
                                {
                                        for (int ion = 0; ion < LIMELM; ++ion)
                                                ionbal.DR_mean_scale[ion] = 0.;
                                }
                                /* multiply guess by scale factor */
                                else if (p.nMatch("SCALE"))
                                {
                                        // there must be at least one scale factor
                                        ionbal.DR_mean_scale[0] = (realnum) p.FFmtRead();
                                        if (p.lgEOL())
                                        {
                                                fprintf(
                                                                ioQQQ,
                                                                " There must be at least one scale factor on the SET RECOMBIANTION MEAN command.\n");
                                                cdEXIT(EXIT_FAILURE);
                                        }
 
                                        for (int ion = 1; ion < LIMELM; ++ion)
                                        {
                                                ionbal.DR_mean_scale[ion] = (realnum) p.FFmtRead();
                                                if (p.lgEOL())
                                                        ionbal.DR_mean_scale[ion]
                                                                        = ionbal.DR_mean_scale[ion - 1];
                                        }
                                        for (int ion = 0; ion < LIMELM; ++ion)
                                        {
                                                if (ionbal.DR_mean_scale[ion] < 0)
                                                {
                                                        fprintf(ioQQQ,
                                                                        " All scale factors on the SET RECOMBIANTION MEAN command must be >=0.\n");
                                                        cdEXIT(EXIT_FAILURE);
                                                }
                                        }
                                }
                                /* option to add log normal noise */
                                else if (p.nMatch("NOISE"))
                                {
                                        ionbal.guess_noise = (realnum) p.FFmtRead();
                                        if (p.lgEOL())
                                                ionbal.guess_noise = 2.;
                                        /* Seed the random-number generator with current time so that
                                         * the numbers will be different every time we run. */
                                        init_genrand((unsigned) time(NULL));
                                }
                                else
                                {
                                        fprintf(ioQQQ, " key OFF or NOISE must appear.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else
                        {
                                fprintf(ioQQQ, " key MEAN must appear.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        fprintf(ioQQQ,
                                        " key DIELECTRonic must appear on set recombination command.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch(" DR "))
        {
                /* set zone thickness by forcing drmax and drmin */
                /* at this stage linear, but default is log */
                radius.sdrmax = p.FFmtRead();
                if (!p.nMatch("LINE"))
                {
                        /* linear was not on command, so default is log */
                        radius.sdrmax = pow(10., radius.sdrmax);
                }
                if (p.lgEOL())
                {
                        p.NoNumb("zone thickness");
                }
 
                radius.lgSdrmaxRel = p.nMatch("RELA");
 
                /* NB these being equal are tested in convinittemp to set dr */
                radius.lgFixed = true;  
                radius.sdrmin = radius.sdrmax;
                radius.lgSdrminRel = radius.lgSdrmaxRel;
                if (!radius.lgSdrmaxRel && radius.sdrmax < DEPTH_OFFSET * 1e4)
                {
                        fprintf(
                                        ioQQQ,
                                        "\n Thicknesses less than about %.0e will NOT give accurate results. If tricking the code\n",
                                        DEPTH_OFFSET * 1e4);
                        fprintf(
                                        ioQQQ,
                                        " into computing emissivities instead of intensities, try to instead use a thickness of unity,\n");
                        fprintf(
                                        ioQQQ,
                                        " and then multiply (divide) the results by the necessary thickness (product of densities).\n\n");
                        cdEXIT(EXIT_FAILURE);
                }
                if (radius.lgSdrmaxRel && (radius.sdrmax <= 0. || radius.sdrmax >= 1.))
                {
                        fprintf(
                                        ioQQQ,
                                        " When using a relative dr, a fraction between 0 and 1 must be entered. Found: %g\n",
                                        radius.sdrmax);
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("DRMA"))
        {
                /* set maximum zone thickness" */
                radius.sdrmax = p.FFmtRead();
                if (p.lgEOL())
                        p.NoNumb("maximum zone thickness");
 
                if ((radius.sdrmax < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
                        radius.sdrmax = pow(10., radius.sdrmax);
 
                // option to set min relative to current radius
                radius.lgSdrmaxRel = p.nMatch("RELA");
                if (radius.lgSdrmaxRel && (radius.sdrmax <= 0. || radius.sdrmax >= 1.))
                {
                        fprintf(
                                        ioQQQ,
                                        " When using a relative drmax, a fraction between 0 and 1 must be entered. Found: %g\n",
                                        radius.sdrmax);
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("DRMI") && p.nMatch("DEPT"))
        {
                /* option to set minimum zone thickness relative to depth */
                radius.sdrmin_rel_depth = p.FFmtRead();
                if (p.lgEOL())
                        p.NoNumb("minimum zone thickness rel to depth");
 
                if ((radius.sdrmin_rel_depth < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
                        radius.sdrmin_rel_depth = pow(10., radius.sdrmin_rel_depth );
 
                if( radius.sdrmin_rel_depth <= 0. || radius.sdrmin_rel_depth >= 1.)
                {
                        fprintf( ioQQQ, " When using a relative drmin, a fraction between 0 and 1 must be entered. Found: %g\n",
                                        radius.sdrmin_rel_depth );
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("DRMI"))
        {
                /* option to set minimum zone thickness */
                radius.sdrmin = p.FFmtRead();
                if (p.lgEOL())
                        p.NoNumb("minimum zone thickness");
 
                if ((radius.sdrmin < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
                        radius.sdrmin = pow(10., radius.sdrmin);
 
                // option to set min relative to current radius
                radius.lgSdrminRel = p.nMatch("RELA");
                radius.lgSMinON = true;
                if (radius.lgSdrminRel && (radius.sdrmin <= 0. || radius.sdrmin >= 1.))
                {
                        fprintf( ioQQQ, " When using a relative drmin, a fraction between 0 and 1 must be entered. Found: %g\n",
                                        radius.sdrmin);
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("FLXF"))
        {
                /* faintest continuum flux to consider */
                rfield.FluxFaint = (realnum) p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("faintest continuum flux to consider");
                }
                if (rfield.FluxFaint < 0.)
                {
                        rfield.FluxFaint = (realnum) pow((realnum) 10.f, rfield.FluxFaint);
                }
        }
 
        else if (p.nMatch("LINE") && p.nMatch("PREC"))
        {
                /* set line precision (number of significant figures)
                 * this affects all aspects of reading and writing lines */
                LineSave.sig_figs = (int) p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("line precision");
                }
                else if (LineSave.sig_figs > 6)
                {
                        fprintf(ioQQQ,
                                        " set line precision currently only works for up to 6 significant figures.\n");
                        cdEXIT(EXIT_FAILURE);
                }
                /* option to print main line array as a single column */
                /* prt.lgPrtLineArray = false; */
        }
 
        /* >>chng 00 dec 08, command added by Peter van Hoof */
        else if (p.nMatch("NFNU"))
        {
                if (p.nMatch(" ADD"))
                {
                        /* option to add a continuum point */
                        double energy = p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("energy");
                        t_PredCont::Inst().add(energy, p.StandardEnergyUnit());
                }
                else
                {
                        /* set nFnu [incident_reflected] [incident_transmitted] [diffuse_inward] [diffuse_outward]
                         *   command for specifying what to include in the nFnu entries in LineSv */
                        /* >>chng 01 nov 12, also accept form with space */
                        /* "incident reflected" keyword */
                        prt.lgSourceReflected = p.nMatch("INCIDENT R") || p.nMatch(
                                        "INCIDENT_R");
                        /* "incident transmitted" keyword */
                        prt.lgSourceTransmitted = p.nMatch("INCIDENT_T") || p.nMatch(
                                        "INCIDENT T");
                        /* "diffuse inward" keyword */
                        prt.lgDiffuseInward = p.nMatch("DIFFUSE_I")
                                        || p.nMatch("DIFFUSE I");
                        /* "diffuse outward" keyword */
                        prt.lgDiffuseOutward = p.nMatch("DIFFUSE_O") || p.nMatch(
                                        "DIFFUSE O");
 
                        /* at least one of these needs to be set ! */
                        if (!(prt.lgSourceReflected || prt.lgSourceTransmitted
                                        || prt.lgDiffuseInward || prt.lgDiffuseOutward))
                        {
                                fprintf(ioQQQ,
                                                " set nFnu expects one or more of the following keywords:\n");
                                fprintf(ioQQQ, " INCIDENT_REFLECTED, INCIDENT_TRANSMITTED, "
                                        "DIFFUSE_INWARD, DIFFUSE_OUTWARD\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }
 
        else if (p.nMatch("IND2"))
        {
                if (p.nMatch(" ON "))
                {
                        /* set flag saying to off or on induced two photon processes */
                        iso_ctrl.lgInd2nu_On = true;
                }
                else if( p.nMatch(" OFF") )
                {
                        iso_ctrl.lgInd2nu_On = false;
                }
                else
                {
                        fprintf( ioQQQ, " set ind2 needs either ON or OFF.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        /* Set the default collision strength for dBase transitions
         * without collision or radiative data */
        else if (p.nMatch("SPECIES"))
        {
                if (p.nMatch(" GBAR"))
                {
                        atmdat.collstrDefault = (long)p.FFmtRead();
                        if (p.lgEOL())
                        {
                                p.NoNumb("the default collision strengths when no collision or radiative data are available");
                        }
                }
                else
                {
                        fprintf(ioQQQ, " SET SPECIES takes option GBAR.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("TEMP"))
        {
                /* set something to do with temperature, currently solver, tolerance, floor */
                if (p.nMatch("FLOO"))
                {
                        StopCalc.TeFloor = (realnum) p.FFmtRead();
                        if (p.lgEOL())
                                p.NoNumb("temperature floor");
 
                        /*  linear option */
                        if (p.nMatch(" LOG") || (StopCalc.TeFloor <= 10. && !p.nMatch(
                                        "LINE")))
                                StopCalc.TeFloor = (realnum) pow(10., StopCalc.TeFloor);
 
                        if (StopCalc.TeFloor < phycon.TEMP_LIMIT_LOW)
                        {
                                fprintf(ioQQQ, " TE < %gK, reset to %gK.\n",
                                                phycon.TEMP_LIMIT_LOW, 1.0001 * phycon.TEMP_LIMIT_LOW);
                                StopCalc.TeFloor = realnum(1.0001 * phycon.TEMP_LIMIT_LOW);
                        }
 
                        if (StopCalc.TeFloor > phycon.TEMP_LIMIT_HIGH)
                        {
                                fprintf(ioQQQ,
                                                " TE > %gK. Cloudy cannot handle this. Bailing out.\n",
                                                phycon.TEMP_LIMIT_HIGH);
                                cdEXIT(EXIT_FAILURE);
                        }
 
                        /* option to vary scale factor */
                        if (optimize.lgVarOn)
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy(optimize.chVarFmt[optimize.nparm],
                                                "SET TEMPERATURE FLOOR %f LOG");
 
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
 
                                /* vary log of temperature  */
                                optimize.vparm[0][optimize.nparm] = (realnum) log10(
                                                StopCalc.TeFloor);
                                optimize.varang[optimize.nparm][0] = (realnum) log10(
                                                1.00001 * phycon.TEMP_LIMIT_LOW);
                                optimize.varang[optimize.nparm][1] = (realnum) log10(
                                                0.99999 * phycon.TEMP_LIMIT_HIGH);
                                optimize.vincr[optimize.nparm] = 0.3f;
 
                                ++optimize.nparm;
                        }
                }
 
                else if (p.nMatch("CONV") || p.nMatch("TOLE"))
                {
                        /* error tolerance in heating cooling match, number is error/total */
                        conv.HeatCoolRelErrorAllowed = (realnum) p.FFmtRead();
                        if (p.lgEOL())
                        {
                                p.NoNumb("heating cooling tolerance");
                        }
                        if (conv.HeatCoolRelErrorAllowed <= 0.)
                        {
                                conv.HeatCoolRelErrorAllowed = (realnum) pow((realnum) 10.f,
                                                conv.HeatCoolRelErrorAllowed);
                        }
                }
 
                else
                {
                        fprintf(ioQQQ,
                                        "\nI did not recognize a keyword on this SET TEMPERATURE command.\n");
                        p.PrintLine(ioQQQ);
                        fprintf(ioQQQ, "The keywords are FLOOr and CONVergence.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("TEST"))
        {
                /* set flag saying to turn on some test - this is in cddefines.h in the global namespace */
                lgTestCodeEnabled = true;
        }
 
        else if (p.nMatch("TRIM"))
        {
                /* set trim upper or lower, for ionization stage trimming
                 * ion trimming ionization trimming 
                 * in routine TrimStage */
                if (p.nMatch("UPPE"))
                {
                        /* set trim upper - either set value or turn off */
                        double save = ionbal.trimhi;
                        ionbal.trimhi = pow(10., p.FFmtRead());
                        if (p.lgEOL() && p.nMatch(" OFF"))
                        {
                                /* turn off upward trimming */
                                p.setEOL(false);
                                ionbal.lgTrimhiOn = false;
                                /* reset high limit to proper value */
                                ionbal.trimhi = save;
                        }
                }
 
                else if (p.nMatch("LOWE"))
                {
                        /* set trim lower */
                        ionbal.trimlo = pow(10., p.FFmtRead());
                }
 
                /* turn off ionization stage trimming */
                else if (p.nMatch("SMAL") || p.nMatch(" OFF"))
                {
                        /* set small limits to both upper and lower limit*/
                        ionbal.trimlo = SMALLFLOAT;
                        ionbal.trimhi = SMALLFLOAT;
                        p.setEOL(false);
                }
 
                else
                {
                        /* set trim upper */
                        ionbal.trimhi = pow(10., p.FFmtRead());
 
                        /* set trim lower to same number */
                        ionbal.trimlo = ionbal.trimhi;
                }
 
                if (p.lgEOL())
                {
                        p.NoNumb("trimming parameter");
                }
 
                if (ionbal.trimlo >= 1. || ionbal.trimhi >= 1.)
                {
                        fprintf(ioQQQ, " number must be negative since log\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("SKIP"))
        {
                /* skip save command, for saving every n't point */
                save.ncSaveSkip = (long) p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("number of points to skip in save");
                }
        }
 
        else if (p.nMatch(" UTA"))
        {
                /* set UTA command, to determine which UTA data set to use 
                 * default is to use Gu et al. 2006 data set */
                if (p.nMatch("GU06"))
                {
                        if (p.nMatch(" OFF"))
                        {
                                /* use Behar et al. 2001 */
                                ionbal.lgInnerShell_Gu06 = false;
                        }
                        else
                        {
                                /* the default, use 
                                 *>>refer       Fe      UTA     Gu, M.F., Holczer, T., Behar, E. & Kahn, S.M. 
                                 *>>refercon    2006, ApJ, 641, 1227 */
                                ionbal.lgInnerShell_Gu06 = true;
                        }
                }
                else if (p.nMatch("KISI"))
                {
                        if (p.nMatch(" OFF"))
                        {
                                /* the default, do not use it */
                                ionbal.lgInnerShell_Kisielius = false;
                        }
                        else
                        {
                                /* use the data from Kisielius et al. 2003, MNRAS, 344, 696 */
                                ionbal.lgInnerShell_Kisielius = true;
                        }
                }
        }
 
        else if (p.nMatch("WEAKH"))
        {
                /* set WeakHeatCool, threshold on save heating and cooling, default 0.05 */
                save.WeakHeatCool = (realnum) p.FFmtRead();
 
                if (p.lgEOL())
                {
                        p.NoNumb("threshold on save heating and cooling");
                }
 
                if (save.WeakHeatCool < 0.)
                {
                        save.WeakHeatCool
                                        = (realnum) pow((realnum) 10.f, save.WeakHeatCool);
                }
        }
 
        else if (p.nMatch("KSHE"))
        {
                /* upper limit to opacities for k-shell ionization */
                continuum.EnergyKshell = (realnum) p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("k-shell ionization opacity limit");
                }
 
                if (continuum.EnergyKshell == 0.)
                {
                        /* arg of 0 causes upper limit to energy array */
                        continuum.EnergyKshell = rfield.egamry;
                }
 
                else if (continuum.EnergyKshell < 194.)
                {
                        fprintf(ioQQQ, " k-shell energy must be greater than 194 Ryd\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else if (p.nMatch("NCHR"))
        {
                /* option to set the number of charge states for grain model */
                double val = p.FFmtRead();
 
                if (p.lgEOL())
                {
                        p.NoNumb("number of charge states");
                }
                else
                {
                        long nChrg = nint(val);
                        if (nChrg < 2 || nChrg > NCHU)
                        {
                                fprintf(ioQQQ,
                                                " illegal value for number of charge states: %ld\n",
                                                nChrg);
                                fprintf(ioQQQ, " choose a value between 2 and %d\n", NCHU);
                                fprintf(ioQQQ,
                                                " or increase NCHS in grainvar.h and recompile\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                        else
                        {
                                SetNChrgStates(nChrg);
                        }
                }
        }
 
        else if (p.nMatch("NEGO"))
        {
                /* save negative opacities if they occur, set negopac */
                opac.lgNegOpacIO = true;
        }
 
        else if (p.nMatch("NEND"))
        {
                /* default limit to number of zones to be computed
                 * only do this if nend is NOT currently left at its default
                 * nend is set to nEndDflt in routine zero
                 * this command only has effect if stop zone not entered */
                if (geometry.lgEndDflt)
                {
                        /* this is default limit to number of zones, change it to this value */
                        geometry.nEndDflt = (long) p.FFmtRead();
                        geometry.lgEndDflt = false;
                        if (p.lgEOL())
                                p.NoNumb("limit to zone number");
 
                        /* now change all limits, for all iterations, to this value */
                        for (long i = 0; i < iterations.iter_malloc; i++)
                                geometry.nend[i] = geometry.nEndDflt;
 
                        if (geometry.nEndDflt > 2000)
                                fprintf(
                                                ioQQQ,
                                                "CAUTION - it will take a lot of memory to save"
                                                        " results for %li zones.  Is this many zones really necessary?\n",
                                                geometry.nEndDflt);
                }
        }
 
        else if (p.nMatch("TSQD"))
        {
                /* upper limit for highest density considered in the 
                 * Peimbert-style t^2 section of the printout */
                peimbt.tsqden = (realnum) p.FFmtRead();
 
                if (p.lgEOL())
                {
                        p.NoNumb("highest density in t^2 section of printout");
                }
                peimbt.tsqden = (realnum) pow((realnum) 10.f, peimbt.tsqden);
        }
 
        else if (p.nMatch("NMAP"))
        {
                /* how many steps in plot or save of heating-cooling map */
                hcmap.nMapStep = (long) p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("steps in heating-cooling map");
                }
        }
 
        else if (p.nMatch("NUME") && p.nMatch("DERI"))
        {
                /* this is an option to use numerical derivatives for heating and cooling */
                NumDeriv.lgNumDeriv = true;
        }
 
        else if (p.nMatch("PATH"))
        {
                fprintf(ioQQQ, " The SET PATH command is no longer supported.\n");
                fprintf(ioQQQ, " Please set the correct path in the file path.h.\n");
                fprintf(ioQQQ,
                                " Or set the environment variable CLOUDY_DATA_PATH.\n Sorry.\n");
                cdEXIT(EXIT_FAILURE);
        }
 
        else if (p.nMatch("PHFI"))
        {
                /* which version of PHFIT to use, 1995 or 1996 */
                ip = (long) p.FFmtRead();
                if (p.lgEOL())
                {
                        p.NoNumb("version of PHFIT");
                }
 
                if (ip == 1995)
                {
                        /* option to go back to old results, pre op */
                        t_ADfA::Inst().set_version(PHFIT95);
                }
                else if (ip == 1996)
                {
                        /* default is to use newer results, including opacity project */
                        t_ADfA::Inst().set_version(PHFIT96);
                }
                else
                {
                        fprintf(ioQQQ, " Two possible values are 1995 and 1996.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        /* set save xxx command */
        else if (p.nMatch("PUNC") || p.nMatch("SAVE"))
        {
                if (p.nMatch("HASH"))
                {
                        /* to use the hash option there must be double quotes on line - were there? */
                        if (!lgQuotesFound)
                        {
                                fprintf(ioQQQ,
                                                " I didn\'t recognize a key on this SET SAVE HASH line.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* specify the terminator between save output sets - normally a set of hash marks */
                        /* 
                         * get any string within double quotes, and return it as
                         * null terminated string
                         * also sets name in OrgCard and chCard to blanks so that
                         * do not trigger off it later */
                        if (strcmp(chString_quotes_lowercase, "return") == 0)
                        {
                                /* special case - return becomes new line */
                                sprintf(save.chHashString, "\n");
                        }
                        else if (strcmp(chString_quotes_lowercase, "time") == 0)
                        {
                                /* special case where output time between iterations */
                                sprintf(save.chHashString, "TIME_DEP");
                        }
                        else
                        {
                                /* usual case, simply copy what is in quotes */
                                strcpy(save.chHashString, chString_quotes_lowercase);
                        }
                }
 
                else if ((p.nMatch("LINE") && p.nMatch("WIDT")) || p.nMatch("RESO")
                                || p.nMatch("LWID"))
                {
                        /* set spectral resolution for contrast in continuum plots */
                        if (p.nMatch(" C "))
                        {
                                fprintf(ioQQQ, " the keyword C is no longer necessary since"
                                        " energy conservation is now supported by default.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                        else if (p.nMatch("SUPP"))
                                /* suppress emission lines in save output */
                                save.Resolution = realnum(0.);
                        else
                        {
                                double number = p.FFmtRead();
                                if (p.lgEOL())
                                        p.NoNumb("line width or resolution");
                                if (number <= 0.)
                                {
                                        fprintf(ioQQQ,
                                                        " line width or resolution must be greater than zero.\n");
                                        cdEXIT(EXIT_FAILURE);
                                }
 
                                if (p.nMatch("RESO"))
                                        /* supply R = lambda/delta(lambda) */
                                        save.Resolution = realnum(number);
                                else
                                        /* FWHM in km/s */
                                        save.Resolution = realnum(SPEEDLIGHT / (number * 1.e5));
                        }
 
                        // option to do exactly the same thing with absorption lines
                        // keyword is absorption 
                        if (p.nMatch("ABSO"))
                                save.ResolutionAbs = save.Resolution;
                }
 
                else if (p.nMatch("PREF"))
                {
                        if( save.nsave > 0 )
                        {
                                fprintf(ioQQQ, " The SET SAVE PREFIX command should precede all save commands.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        // specify a prefix before all save filenames
                        save.chFilenamePrefix = chString_quotes_lowercase;
                }
 
                else if (p.nMatch("FLUS"))
                {
                        /* flus the output after every iteration */
                        save.lgFLUSH = true;
                }
 
                else
                {
                        fprintf(ioQQQ,
                                        " There should have been an option on this command.\n");
                        fprintf(ioQQQ,
                                        " Valid options for SET SAVE are summarized in Hazy 1 "
                                                "Miscellaneous commands.\n");
                        fprintf(ioQQQ,
                                        " The SET PUNCHLWIDTH command is now SET SAVE LINE WIDTH.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        /* set continuum .... options */
        else if (p.nMatch("CONT"))
        {
                if (p.nMatch("FEII") || p.nMatch("FE I"))
                {
                        /* set lower, upper wavelength range, and number of bins, for
                         * save feii continuum command */
                        FeII.feconwlLo = (realnum) p.FFmtRead();
                        FeII.feconwlHi = (realnum) p.FFmtRead();
                        FeII.nfe2con = (long) p.FFmtRead();
                        if (p.lgEOL())
                        {
                                fprintf(ioQQQ,
                                                "PROBLEM The set FeII continuum command must have"
                                                        " three numbers, the lower and upper wavelength range in Angstroms"
                                                        " and the number of bins to divide this into.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
 
                        if (FeII.feconwlLo >= FeII.feconwlHi)
                        {
                                fprintf(ioQQQ, "PROBLEM The first two numbers on the set "
                                        "FeII continuum command must be the lower and upper "
                                        "wavelength range in Angstroms and the first must be less "
                                        "than the second.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                        if (FeII.nfe2con < 2)
                        {
                                fprintf(ioQQQ,
                                                "PROBLEM The third number on the set FeII continuum "
                                                        "command must be the number of bins to divide the range into and"
                                                        " it must be greater than 1.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
 
                }
 
                else if (p.nMatch("RESO"))
                {
                        /* set resolution, get factor that will multiply continuum resolution that
                         * is contained in the file continuum_mesh.ini */
                        continuum.ResolutionScaleFactor = p.FFmtRead();
                        if (p.lgEOL())
                        {
                                p.NoNumb("continuum resolution scale");
                        }
                        /* negative numbers were logs */
                        if (continuum.ResolutionScaleFactor <= 0.)
                                continuum.ResolutionScaleFactor = pow(10.,
                                                continuum.ResolutionScaleFactor);
                }
 
                else if (p.nMatch("SHIE"))
                {
                        /* set continuum shielding function */
                        /* these are all possible values of rt.nLineContShield,
                         * first is default, these are set with set continuum shielding */
                        /*#define LINE_CONT_SHIELD_PESC 1*/
                        /*#define LINE_CONT_SHIELD_FEDERMAN     2*/
                        /*#define LINE_CONT_SHIELD_FERLAND      3*/
                        if (p.nMatch("PESC"))
                        {
                                /* this uses an inward looking escape probability */
                                rt.nLineContShield = LINE_CONT_SHIELD_PESC;
                        }
                        else if (p.nMatch("FEDE"))
                        {
                                /* set continuum shielding Federman,
                                 * this is the default, and uses the appendix of
                                 * >>refer      RT      continuum shielding     Federman, S.R., Glassgold, A.E., & 
                                 * >>refercon   Kwan, J. 1979, ApJ, 227, 466*/
                                rt.nLineContShield = LINE_CONT_SHIELD_FEDERMAN;
                        }
                        else if (p.nMatch("FERL"))
                        {
                                rt.nLineContShield = LINE_CONT_SHIELD_FERLAND;
                        }
                        else if (p.nMatch("NONE"))
                        {
                                /* turn off continuum shielding */
                                rt.nLineContShield = 0;
                        }
                        else
                        {
                                fprintf(ioQQQ,
                                                " I didn\'t recognize a key on this SET CONTINUUM SHIELDing line.\n");
                                fprintf(ioQQQ,
                                                " The ones I know about are: PESC, FEDErman, & FERLand.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
 
                else
                {
                        fprintf(ioQQQ,
                                        " I didn\'t recognize a key on this SET CONTINUUM line.\n");
                        fprintf(ioQQQ,
                                        " The ones I know about are: RESOlution and SHIEld.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
 
        else
        {
                fprintf(ioQQQ, " I didn\'t recognize a key on this SET command.\n");
                cdEXIT(EXIT_FAILURE);
        }
 
        return;
}