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

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ParseSet scan parameters off SET command */
#include "cddefines.h"
#include "physconst.h"
#include "geometry.h"
#include "input.h"
#include "prt.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 "iterations.h"
#include "conv.h"
#include "secondaries.h"
#include "rfield.h"
#include "ionbal.h"
#include "numderiv.h"
#include "dynamics.h"
#include "iso.h"
#include "punch.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 "parse.h"
#include "lines.h"
#include "assertresults.h"
#include "thirdparty.h"

void ParseSet(char *chCard )
{
        bool lgEOL;
        long int i, 
                ip;
        char chString_quotes_lowercase[INPUT_LINE_LENGTH];
        bool lgQuotesFound;

        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*/
        lgQuotesFound = true;
        if( GetQuote( chString_quotes_lowercase , chCard , false ) )
                lgQuotesFound = false;

        /* commands to set certain variables, "SET XXX=" */
        if( nMatch("ASSE",chCard) && nMatch("ABOR",chCard) )
        {
                /* set assert abort command, to tell the code to raise SIGABRT when an
                 * assert is thrown - this can be caught by the debugger. The keyword
                 * _FPE is accepted for backward compatibility */
                cpu.setAssertAbort( true );
        }

        else if( nMatch("ASSE",chCard) &&nMatch("SCIE",chCard) )
        {
                /* print asserts using scientific notation.  Useful when an asserted value is 
                * less than 10^-4 of the normalization line. */
                lgPrtSciNot = true;
        }

        else if( nMatch("ATOM",chCard) && nMatch( "DATA", chCard ) )
        {
                /* set atomic data */
                long int nelem , ion;
                bool lgAs=false , lgCS=false , lgDR=false;
                const char *chMole;

                /* As? */
                if( nMatch( " AS " , chCard ) )
                        lgAs = true;
                /* DR - dielectronic recombination */
                else if( nMatch( " DR " , chCard ) )
                        lgDR = true;
                /* lgCS collision strengths */
                else if( nMatch( " CS " , chCard ) )
                        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);
                }

                /* ion or molecule? */
                if( nMatch(" ION",chCard) )
                {
                        /* set atomic data - must have an element name and possibly an
                         * ionization stage - nelem is atomic number on C scale */
                        if( (nelem = GetElem(chCard))<0 )
                        {
                                fprintf( ioQQQ, " An element name must appear on this line\n Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        i = 5;
                        /* make sure that the ionization stage is ok if it is specified 
                         * this is entered on physics scale, 1 is atom, 2 first ion */
                        ion = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( !lgEOL && (ion< 1 || ion>nelem+1) )
                        {
                                fprintf( ioQQQ, " The ionization stage is not valid\n Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* two possible sets of transition probabilities for OII */
                        if( nelem==ipOXYGEN && ion == 2 && lgAs )
                        {
                                if( nMatch( " NEW" , chCard ) )
                                {
                                        dense.lgAsChoose[nelem][ion-1] = true;
                                }
                                else if( nMatch( " OLD" , chCard ) )
                                {
                                        /*>>chng 06 mar 30, this is the default */
                                        dense.lgAsChoose[nelem][ion-1] = false;
                                }
                                else
                                {
                                        fprintf( ioQQQ, " The keyword old or new must appear\n Sorry.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else if( nelem==ipSULPHUR && lgDR )
                        {
                                /* change DR data set for S 
                                 * three cases for S DR - ionbal.nDR_S_guess
                                 * 0, default larger of guess and Badnell
                                 * 1, pure Badnell
                                 * 3, scaled oxygen */
                                if( nMatch( " MIX" , chCard ) )
                                        /* this is the default, larger of Badnell or guess */
                                        ionbal.nDR_S_guess = 0;
                                else if( nMatch( "PURE" , chCard ) )
                                        /* use only Badnell 1991 */
                                        ionbal.nDR_S_guess = 1;
                                else if( nMatch( "SCAL" , chCard ) )
                                {
                                        /* scaled oxygen */
                                        ionbal.nDR_S_guess = 2;
                                        i = 5;
                                        ionbal.DR_S_scale[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                        if( lgEOL )
                                                NoNumb(chCard);
                                        for( ion=1; ion<5; ++ion )
                                        {
                                                /* optionally get rest - if too few specified then use
                                                 * previously set value */
                                                ionbal.DR_S_scale[ion] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                                if( lgEOL )
                                                        ionbal.DR_S_scale[ion] = ionbal.DR_S_scale[ion-1];
                                        }
                                }
                                else
                                {
                                        /* disaster - no match */
                                        fprintf( 
                                                ioQQQ, " One of the keywords MIX, PURE, or SCALe must"
                                                "appear on the set atomic physics sulphur dr command.\n"
                                                "Sorry.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }
                        else
                        {
                                fprintf( ioQQQ, " None of the valid set atomic data atoms/ions were found\n Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        /* a molecule */
                        if( nMatch(" H2 ",chCard) )
                        {
                                /* 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( nMatch(" H " , chCard ) && lgCS )
                                {
                                        /* change the H - H2 collision data set */
                                        i = 5;
                                        ion = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                        if( ion!=2 )
                                                TotalInsanity();
                                        long int nYear = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                        if( lgEOL )
                                                NoNumb( chCard );
                                        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;
                                        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;
                                        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( nMatch(" He" , chCard ) && lgCS )
                                {
                                        /* change the He - H2 collision data set */
                                        if( nMatch("ORNL" , chCard ) )
                                        {
                                                /* use the coefficients from
                                                *>>refer        H2      He collision    Lee et al in preparation */
                                                mole.lgH2_He_ORNL = true;
                                        }
                                        else if( nMatch( "BOUR",chCard ) )
                                        {
                                                /* use the coefficients from
                                                *>>refer        H2      He collision    Le Bourlot, J., Pineau des Forets, 
                                                *>>refercon     G., & Flower, D.R. 1999, MNRAS, 305, 802*/
                                                mole.lgH2_He_ORNL = false;
                                        }
                                        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
                                TotalInsanity();
                }
        }

        else if( nMatch(" CHA",chCard) && !nMatch( "HO ", chCard ) )
        {
                /* set log of minimum charge transfer rate for high ions and H
                 * default of 1.92e-10 set in zero */
                i = 5;
                atmdat.HCTAlex = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
                if( atmdat.HCTAlex < 0. )
                {
                        atmdat.HCTAlex = pow(10.,atmdat.HCTAlex);
                }
        }

        else if( nMatch("CHEM",chCard ) && !nMatch( "HO ", chCard ) )
        {
                /* turn on Steve Federman's chemistry */
                if( nMatch("FEDE",chCard ) )
                {
                        if( nMatch( " ON " , chCard ) )
                        {
                                /* This turns on the diffuse cloud chemical rates of 
                                 * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
                                co.lgFederman = true;
                        }
                        else if( nMatch( " OFF" , chCard ) )
                        {
                                co.lgFederman = false;
                        }
                        else
                        {
                                /* this is the default when command used - true */
                                co.lgFederman = true;
                        }
                }
                /* >>chng 06 may 30 --NPA.  Turn on non-equilibrium chemistry */
                else if( nMatch(" NON",chCard ) && nMatch( "EQUI", chCard ))
                {

                        /* 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 */

                        co.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( nMatch("NEUT",chCard ) )
                        {
                                if( nMatch( " ON " , chCard ) )
                                {
                                        /* This turns on the diffuse cloud chemical rates of 
                                         * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
                                        co.lgNeutrals = true;
                                }
                                else if( nMatch( " OFF" , chCard ) )
                                {
                                        co.lgNeutrals = false;
                                }
                                else
                                {
                                        /* this is the default when command used - true */
                                        co.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( nMatch("PROT",chCard )  && nMatch( "ELIM", chCard ) )
                {
                        if( nMatch( " ON " , chCard ) )
                        {
                                /* This turns on the diffuse cloud chemical rates of 
                                 * >>refer      CO      chemistry       Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/
                                co.lgProtElim = true;
                        }
                        else if( nMatch( " OFF" , chCard ) )
                        {
                                co.lgProtElim = false;
                        }
                        else
                        {
                                /* this is the default when command used - true */
                                co.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( nMatch("COLL",chCard) && nMatch("STRE",chCard) && nMatch("AVER",chCard) )
        {
                if( nMatch(" OFF",chCard) )
                {
                         iso.lgCollStrenThermAver = false;
                }
                else
                {
                        /* this is default behavior of this command */
                         iso.lgCollStrenThermAver = true;
                }
        }

        else if( nMatch("COVE",chCard) )
        {
                iterations.lgConverge_set = true;
                /* set coverage - limit number of iterations and zones */
                if( nMatch("FAST",chCard) )
                {
                        iterations.lim_zone = 1;
                        iterations.lim_iter = 0;
                }
                else
                {
                        iterations.lim_zone = 10;
                        iterations.lim_iter = 1;
                }
        }

        else if( nMatch("CSUP",chCard) )
        {
                /* force secondary ionization rate, log entered */
                i = 5;
                secondaries.SetCsupra = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                secondaries.lgCSetOn = true;
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
                secondaries.SetCsupra = (realnum)pow((realnum)10.f,secondaries.SetCsupra);
        }

        else if( nMatch(" D/H",chCard) )
        {
                /* set deuterium abundance, D to H ratio */
                i = 5;
                hydro.D2H_ratio = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( hydro.D2H_ratio <= 0. || nMatch( " LOG", chCard ) )
                {
                        hydro.D2H_ratio = pow(10.,hydro.D2H_ratio);
                }
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
        }

        else if( nMatch( " HO " , chCard) && nMatch( "CHAR" , chCard) )
        {
                /* which solver will include H + O charge transfer?  done in
                 * chemistry by default */
                if( nMatch( "CHEM" , chCard ) )
                {
                        /* this is the default */
                        ionbal.lgHO_ct_chem = true;
                }
                else if( nMatch( "IONI" , chCard ) )
                {
                        /* do it in the ionization solver */
                        ionbal.lgHO_ct_chem = false;
                }
                else
                {
                        fprintf( ioQQQ, " I did not recognize a subkey on this SET OH CHARge transfer line.\n" );
                        fprintf( ioQQQ, " Valid keys are CHEMistry and IONIzation.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("12C1",chCard) )
        {
                /* set the 12C to 13C abundance ratio - default is 30 */
                i = 5;
                /* first two numbers on the line are 12 and 13 - we don't want them */
                co.C12_C13_isotope_ratio = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                co.C12_C13_isotope_ratio = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                /* now we can get the ratio */
                co.C12_C13_isotope_ratio = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                        NoNumb(chCard);

                if( co.C12_C13_isotope_ratio <= 0. || nMatch( " LOG", chCard ) )
                        co.C12_C13_isotope_ratio = (realnum)pow((realnum)10.f,co.C12_C13_isotope_ratio);
        }

        /* set dynamics ... */
        else if( nMatch("DYNA",chCard) )
        {
                /* set dynamics advection length */
                if( nMatch("ADVE",chCard) && nMatch("LENG",chCard) )
                {
                        /* <0 => relative fraction of length, + val in cm */
                        i = 5;
                        dynamics.AdvecLengthInit = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb( chCard );
                        /* if fraction is present, then number was linear fraction, if not present
                        * then physical length in cm, log10 */
                        if( nMatch("FRAC",chCard) )
                        {
                                /* 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( nMatch("PRES",chCard) && nMatch("MODE",chCard) )
                {
                        dynamics.lgSetPresMode = true;
                        if( nMatch("SUBS",chCard) )
                        {
                                /* subsonic */
                                strcpy( dynamics.chPresMode , "subsonic" );
                        }
                        else if( nMatch("SUPE",chCard) )
                        {
                                /* supersonic */
                                strcpy( dynamics.chPresMode , "supersonic" );
                        }
                        else if( nMatch("STRO",chCard) )
                        {
                                /* strong d */
                                strcpy( dynamics.chPresMode , "strongd" );
                        }
                        else if( nMatch("ORIG",chCard) )
                        {
                                /* original */
                                strcpy( dynamics.chPresMode , "original" );
                        }
                }
                else if( nMatch("ANTI",chCard) && nMatch("DEPT",chCard) )
                {
                        dynamics.lgSetPresMode = true;
                        strcpy( dynamics.chPresMode , "antishock" );
                        /* shock depth */
                        i = 5;
                        /* get log of shock depth in cm */
                        dynamics.ShockDepth = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb( chCard );
                        dynamics.ShockDepth = pow( 10., dynamics.ShockDepth );
                }
                else if( nMatch("ANTI",chCard) && nMatch("MACH",chCard) )
                {
                        dynamics.lgSetPresMode = true;
                        strcpy( dynamics.chPresMode , "antishock-by-mach" );
                        /* Mach number */
                        i = 5;
                        /* get (isothermal) Mach number where we want antishock to occur */
                        dynamics.ShockMach = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb( chCard );
                }
                else if( nMatch("RELA",chCard) )
                {
                        /* set how many iterations we will start with, before allowing
                         * changes.  This allows the solution to relax to an equilibrium */
                        i = 5;
                        dynamics.n_initial_relax = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb( chCard );
                        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( nMatch("SHOC",chCard) && nMatch("DEPT",chCard) )
                {
                        dynamics.lgSetPresMode = true;
                        strcpy( dynamics.chPresMode , "shock" );
                        /* shock depth */
                        i = 5;
                        /* get log of shock depth in cm */
                        dynamics.ShockDepth = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb( chCard );
                        dynamics.ShockDepth = pow( 10., dynamics.ShockDepth );
                }
                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( nMatch("DIDZ",chCard) )
        {
                /* 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 */
                i = 5;
                radius.drChange = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( radius.drChange <= 0. )
                {
                        radius.drChange = (realnum)pow((realnum)10.f,radius.drChange);
                }
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
        }

        /* something to do with electron density */
        else if( nMatch("EDEN",chCard) )
        {
                /* >>chng 02 apr 20, from ERROR to TOLERANCE to be parallel to set temp equivalent,
                 * >>chng 04 jun 03, but also accept error as keyword */
                if( nMatch("CONV",chCard) || nMatch("ERRO",chCard))
                {
                        /* keyword is eden convergence  
                         * set tolerance in eden match */
                        i = 5;
                        conv.EdenErrorAllowed = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        if( conv.EdenErrorAllowed < 0. )
                        {
                                conv.EdenErrorAllowed = pow(10.,conv.EdenErrorAllowed);
                        }
                }

                else if( nMatch("SOLV",chCard) )
                {
                        /* the electron density solver */
                        if( nMatch("NEW",chCard) ) 
                        {
                                /* new method */
                                strcpy( conv.chSolverEden , "new" );
                        }
                        else
                        {
                                /* default method */
                                strcpy( conv.chSolverEden , "simple" );
                        }
                }
                else 
                {
                        /* no keyword, assume log of electron density */
                        i = 5;
                        /* set the electron density */
                        dense.EdenSet = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                        /* warn that this model is meaningless */
                        phycon.lgPhysOK = false;
                }
        }

        else if( nMatch("FINE",chCard ) && nMatch("CONT",chCard ) )
        {
                /* 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 = GetElem(chCard))<0 )
                {
                        fprintf( ioQQQ, " An element name must appear on this line\n Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
                i = 5;
                /* 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)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                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);
                }
        }

        /* set grain command - but not set H2 grain command */
        else if( nMatch("GRAI",chCard ) && !nMatch(" H2 ",chCard) )
        {
                if( nMatch("HEAT",chCard ) )
                {
                        /* scale factor to change grain heating as per Allers et al. */
                        i = 5;
                        gv.GrainHeatScaleFactor = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        /* warn that this model is not the best we can do */
                        phycon.lgPhysOK = false;
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                }
                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( nMatch("LEID",chCard ) && nMatch("HACK",chCard ) )
        {
                if( nMatch( "H2* " , chCard ) && nMatch( " OFF" , chCard ) )
                {
                        /* 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( nMatch( "CR " , chCard ) && nMatch( " OFF" , chCard ) )
                {
                        /* the CR Leiden hack option - turn off CR excitations of H2 */
                        hmi.lgLeidenCRHack = false;

                }
                else if( nMatch("RATE",chCard ) &&  nMatch("UMIS",chCard ))
                {
                        /* This command will use the rates given in the UMIST database,  It
                         * will set to zero many reactions in hmole_step.c that are not
                         * in UMIST */
                        co.lgUMISTrates = false;
                }
        }

        /* set H2 ... */
        else if( nMatch(" H2 ",chCard) )
        {
                i = 5;
                ip = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                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( nMatch("SOLO",chCard) || ( nMatch("SMAL",chCard)&&nMatch("MODE",chCard) ) )
                {
                        if( nMatch("SOLO",chCard) )
                        {
                                /* 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( nMatch("TH85", chCard ) )
                        {
                                /* rate from is eqn A8 of Tielens and Hollenbach 85a, */
                                hmi.chH2_small_model_type = 'T';
                        }
                        else if( nMatch( " BHT" , chCard ) )
                        {
                                /* 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( nMatch( "BD96" , chCard ) )
                        {
                                /* 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( nMatch( "ELWE" , chCard ) )
                        {
                                /* 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( nMatch("GRAI",chCard ) && nMatch("FORM",chCard ) && nMatch("PUMP",chCard )  )
                {
                        if( nMatch( "DB96" , chCard ) )
                        {
                                /* Draine & Bertoldi 1996 */
                                hmi.chGrainFormPump = 'D';
                        }
                        else if( nMatch( "TAKA" , chCard ) )
                        {
                                /* Takahashi 2001 */
                                hmi.chGrainFormPump = 'T';
                        }
                        else if( nMatch( "THER" , chCard ) )
                        {
                                /* 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( nMatch("JURA",chCard) )
                {
                        if( nMatch("TH85", chCard ) )
                        {
                                /* rate from is eqn A8 of Tielens and Hollenbach 85a*/
                                hmi.chJura = 'T';
                        }
                        else if( nMatch( "CT02" , chCard ) )
                        {
                                /* this rate is equation Cazeux & Tielens */
                                hmi.chJura = 'C';
                        }
                        else if( nMatch( "SN99" , chCard ) )
                        {
                                /* this rate is from Sternberg & Neufeld 99 */
                                hmi.chJura = 'S';
                        }
                        else if( nMatch( "RATE" , chCard ) )
                        {
                                /* 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.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) );
                                if( 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( nMatch( "SCAL" , chCard ) )
                        {
                                /* this is a scale factor to multiply the Jura rate */
                                hmi.ScaleJura = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                /* log or negative number means log was entered */
                                if( nMatch( " LOG" , chCard ) || hmi.ScaleJura<0. )
                                {
                                        hmi.ScaleJura = (realnum)pow( 10., (double)hmi.ScaleJura );
                                }
                                if( lgEOL )
                                {
                                        NoNumb( chCard );
                                }

                                /* option to vary scale factor */
                                if( optimize.lgVarOn )
                                {
                                        optimize.nvarxt[optimize.nparm] = 1;
                                        strcpy( optimize.chVarFmt[optimize.nparm], "SET H2 JURA SCALE %f" );

                                        /* 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( nMatch(" TAD",chCard) )
                {
                        hmi.Tad = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb(chCard);
                        /* log if <= 10. unless linear key appears too */
                        if( hmi.Tad <=10. && !nMatch("LINE",chCard) )
                                hmi.Tad = (realnum)pow((realnum)10.f,hmi.Tad);
                }

                else if( nMatch("FRAC",chCard ) )
                {
                        /* 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 = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb( chCard );

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

        /* 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( nMatch("HCOR",chCard)  )
        {
                i = 5;
                dense.HCorrFac = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                        NoNumb( chCard );
        }

        else if( nMatch(" PAH",chCard)  )
        {
                /* set one of several possible PAH abundance distribution functions */
                if( nMatch(" H0 " , chCard ) )
                {
                        /* the default, abundance depends on H0 / Htot */
                        strcpy( gv.chPAH_abundance_fcn , "H0" );
                }
                else if( nMatch("CONS" , chCard ) )
                {
                        /* constant PAH abundance */
                        strcpy( gv.chPAH_abundance_fcn , "CON" );
                }

                else if(nMatch("BAKE",chCard ) )
                {
                        /* 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 */
                        phycon.lgPhysOK = false;
                }
                else
                {
                        fprintf( ioQQQ, " one of the keywords H0, BAKES, or CONStant must appear.\n Sorry." );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("PRES",chCard) && nMatch("IONI",chCard) )
        {
                /* set limit to number of calls from pressure to ionize solver,
                 * this set limit to conv.nPres2Ioniz */
                i = 5;
                conv.limPres2Ioniz = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
                else if( conv.limPres2Ioniz <= 0 )
                {
                        fprintf( ioQQQ, " The limit must be greater than zero.\n Sorry." );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        /* something to do with pressure */
        else if( nMatch("PRES",chCard) )
        {
                /* tolerance on pressure convergence */
                if( nMatch("CONV",chCard) )
                {
                        /* keyword is tolerance 
                         * set tolerance or relative error in pressure match */
                        i = 5;
                        conv.PressureErrorAllowed = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb(chCard);

                        if( conv.PressureErrorAllowed < 0. )
                                conv.PressureErrorAllowed = (realnum)pow((realnum)10.f,conv.PressureErrorAllowed);
                }

                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.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else if( nMatch("RECO",chCard) && nMatch("MBIN",chCard) )
        {
                /* dielectronic recombination */
                if( nMatch("DIEL",chCard) && nMatch("ECTR",chCard) )
                {
                        /* change various factors for dielectronic recombination */
                        if( nMatch("KLUD",chCard) )
                        {
                                if( nMatch("BADN",chCard) )
                                {
                                        /* set true to use Badnell mean in place of existing hacks */
                                        ionbal.lg_use_DR_Badnell_rate_coef_mean_ion = true;
                                }
                                else
                                {
                                        int j;
                                        /* option to turn off guess of dielectronic rec coefficient for 3rd row elements */
                                        i = 3;
                                        /* this is first call, no number, lgEOL true but zero returned, the intended effect*/
                                        ionbal.GuessDiel[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                        /* >>chng 06 jan 30 also turn off guess if GuessDiel set zero */
                                        if( ionbal.GuessDiel[0]<=0. || nMatch(" OFF",chCard) )
                                                ionbal.lg_guess_coef = false;

                                        for( j=1; j<4; ++j )
                                        {
                                                ionbal.GuessDiel[j] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                                /* here lgEOL means to use previous values */
                                                if( lgEOL )
                                                        ionbal.GuessDiel[j] = ionbal.GuessDiel[j-1];
                                        }
                                }
                                /* option to add log normal noise */
                                if( nMatch("NOISE" , chCard ) ) 
                                {
                                        i = 3;
                                        ionbal.guess_noise = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                        if( 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
                                {
                                        ionbal.guess_noise = 0.;
                                }
                        }

                        else if( nMatch("BURG"  , chCard) )
                        {
                                /* modify suppression of burgess dielectronic recombinations */
                                if( nMatch(" ON ",chCard) )
                                {
                                        /* leave suppression on - this is default state */
                                        ionbal.lgSupDie[0] = true;
                                }

                                else if( nMatch(" OFF",chCard) )
                                {
                                        /* turn suppression off */
                                        ionbal.lgSupDie[0] = false;
                                }

                                else
                                {
                                        fprintf( ioQQQ, " flag ON or OFF must appear.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }

                        else if( nMatch("NUSS"  , chCard) )
                        {
                                /* modify suppression of Nussbaumer and Storey dielectronic recombination */
                                if( nMatch(" ON ",chCard) )
                                {
                                        /* turn suppression on */
                                        ionbal.lgSupDie[1] = true;
                                }
                                else if( nMatch(" OFF",chCard) )
                                {
                                        /* leave suppression off - this is default state */
                                        ionbal.lgSupDie[1] = false;
                                }
                                else
                                {
                                        fprintf( ioQQQ, " flag ON or OFF must appear.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                        }

                        else if( nMatch("BADN" , chCard ) )
                        {
                                /* use the new (in early 2006) Badnell DR rates */
                                if( nMatch( " OFF" , chCard ) )
                                {
                                        ionbal.lgDR_recom_Badnell_use = false;
                                }
                                else
                                {
                                        ionbal.lgDR_recom_Badnell_use = true;
                                }

                                /* option to print rates then exit */
                                if( nMatch("PRIN" , chCard ) )
                                        ionbal.lgRecom_Badnell_print = true;
                        }
                        else
                        {
                                fprintf( ioQQQ, " key KLUDge, BURGess, NUSSbaumer, or BADNell must appear.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                /* radiative recombination */
                else if( nMatch("RADI",chCard) && nMatch("ATIV",chCard) )
                {
                        if( nMatch("BADN" , chCard ) )
                        {
                                /* use the new (in early 2006) Badnell DR rates */
                                if( nMatch( " OFF" , chCard ) )
                                {
                                        ionbal.lgRR_recom_Badnell_use = false;
                                }
                                else
                                {
                                        ionbal.lgRR_recom_Badnell_use = true;
                                }

                                /* option to print rates then exit */
                                if( nMatch("PRIN" , chCard ) )
                                        ionbal.lgRecom_Badnell_print = true;
                        }
                        else
                        {
                                fprintf( ioQQQ, " key BADNell must appear.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        fprintf( ioQQQ, " key RADIATIve or DIELECTRonic must appear on set recombination command.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("SPEC",chCard) )
        {
                /* "set spectrum" for optional parameters for punch spectrum command */
                if( nMatch("RANG" , chCard ) )
                {
                        /* energy range option - argument is energy in Rydbergs */
                        i = 5;
                        /* >>chng 05 aug 25, rm pow(10., - was wrong */
                        punch.cp_range_min[punch.cp_npun] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        punch.cp_range_max[punch.cp_npun] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* confirm that wavelengths are in correct order */
                        if( punch.cp_range_min[punch.cp_npun] >= punch.cp_range_max[punch.cp_npun] )
                        {
                                fprintf( ioQQQ, " The limits must be in increasing order.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( nMatch("RESO" , chCard ) )
                {
                        /* resolving power, default is zero, which means leave at native resolution */
                        /* >>chng 05 aug 25, following had pow(number), so crashed for high resolving power
                         * bug caught by Yan Changshou */
                        i = 5;
                        punch.cp_resolving_power[punch.cp_npun] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                }
        }

        else if( nMatch(" DR ",chCard) )
        {
                /* set zone thickness by forcing drmax and drmin */
                i = 5;
                /* at this stage linear, but default is log */
                radius.sdrmax = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( !nMatch("LINE",chCard) ) 
                {
                        /* linear was not on command, so default is log */
                        radius.sdrmax = pow( 10. , radius.sdrmax );
                }
                if( lgEOL )
                {
                        NoNumb(chCard);
                }

                /* NB these being equal are tested in convinittemp to set dr */
                radius.sdrmin = radius.sdrmax;
                if( 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);
                }
        }

        else if( nMatch("DRMA",chCard) )
        {
                /* set maximum zone thickness */
                i = 5;
                radius.sdrmax = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                        NoNumb(chCard);

                /* log if less than 38 or if log keyword occurs */
                if( radius.sdrmax < log10(BIGFLOAT) || nMatch(" LOG",chCard) )
                        radius.sdrmax = pow(10.,radius.sdrmax);
        }

        else if( nMatch("DRMI",chCard) )
        {
                /* option to set minimum zone thickness */
                i = 5;
                radius.sdrmin = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                        NoNumb(chCard);

                /* log if less than 38 or if log keyword occurs */
                if( radius.sdrmin < log10(BIGFLOAT ) || nMatch(" LOG",chCard) )
                        radius.sdrmin = pow(10.,radius.sdrmin);

                radius.lgSMinON = true;
        }

        else if( nMatch("FLXF",chCard) )
        {
                /* faintest continuum flux to consider */
                i = 5;
                rfield.FluxFaint = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
                if( rfield.FluxFaint < 0. )
                {
                        rfield.FluxFaint = (realnum)pow((realnum)10.f,rfield.FluxFaint);
                }
        }

        else if( nMatch("LINE",chCard) && nMatch("PREC",chCard) )
        {
                /* set line precision (number of significant figures)
                 * this affects all aspects of reading and writing lines */
                i = 5;
                LineSave.sig_figs = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
                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( nMatch("NFNU",chCard) )
        {
                /* 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 = nMatch("NT R",chCard) || nMatch("NT_R",chCard);
                /* "incident transmitted" keyword */
                prt.lgSourceTransmitted = nMatch("NT_T",chCard) || nMatch("NT T",chCard);
                /* "diffuse inward" keyword */
                prt.lgDiffuseInward = nMatch("SE_I",chCard) || nMatch("SE I",chCard);
                /* "diffuse outward" keyword */
                prt.lgDiffuseOutward = nMatch("SE_O",chCard) || nMatch("SE O",chCard);

                /* 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);
                }
                /* automatically print the nFnu items in the output - it is not necessary to also include
                 * a print continua command if this is entered */
                prt.lgPrnDiff = true;
        }

        else if( nMatch("IND2",chCard) )
        {
                if( nMatch(" ON ",chCard) )
                {
                        /* set flag saying to off or on induced two photon processes */
                        iso.lgInd2nu_On = true;
                }
                else if( nMatch(" OFF",chCard) )
                {
                        iso.lgInd2nu_On = false;
                }
                else
                {
                        fprintf( ioQQQ, " set ind2 needs either ON or OFF.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("TEMP",chCard) )
        {
                /* set something to do with temperature, currently solver, tolerance, floor */
                if( nMatch("SOLV",chCard) )
                {
                        /* solver */
                        /* the electron density solver */
                        if( nMatch("NEW",chCard) ) 
                        {
                                /* new method */
                                strcpy( conv.chSolverTemp , "new" );
                        }
                        else
                        {
                                /* default method */
                                strcpy( conv.chSolverTemp , "simple" );
                        }
                }

                else if( nMatch("FLOO",chCard) )
                {
                        i = 5;
                        StopCalc.TeFloor = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                        /*  linear option */
                        if( StopCalc.TeFloor <= 10. && !nMatch("LINE",chCard) )
                        {
                                StopCalc.TeFloor = (realnum)pow(10.,StopCalc.TeFloor);
                        }

                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        if( StopCalc.TeFloor < 2.8 )
                        {
                                fprintf( ioQQQ, " TE < 3K, reset to 2.8K.\n" );
                                StopCalc.TeFloor = 2.8f;
                        }
                }

                else if( nMatch("CONV",chCard) || nMatch("TOLE",chCard) )
                {
                        /* error tolerance in heating cooling match, number is error/total */
                        i = 5;
                        conv.HeatCoolRelErrorAllowed = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        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 TEPERATURE command.\n%s\n" , chCard);
                        fprintf( ioQQQ, "The keywords are SOLVer, FLOOr, and CONVergence.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("TEST",chCard) )
        {
                /* set flag saying to turn on some test - this is in cddefines.h in the global namespace */
                lgTestCodeEnabled = true;
        }

        else if( nMatch("TRIM",chCard) )
        {
                /* set trim upper or lower, for ionization stage trimming
                 * ion trimming ionization trimming 
                 * in routine TrimStage */
                i = 5;
                if( nMatch("UPPE",chCard) )
                {
                        /* set trim upper - either set value or turn off */
                        double save = ionbal.trimhi;
                        ionbal.trimhi = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                        if( lgEOL  && nMatch(" OFF",chCard) )
                        {
                                /* turn off upward trimming */
                                lgEOL = false;
                                ionbal.lgTrimhiOn = false;
                                /* reset high limit to proper value */
                                ionbal.trimhi = save;
                        }
                }

                else if( nMatch("LOWE",chCard) )
                {
                        /* set trim lower */
                        ionbal.trimlo = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                }

                /* turn off ionization stage trimming */
                else if( nMatch("SMAL",chCard) || nMatch(" OFF",chCard) )
                {
                        /* set small limits to both upper and lower limit*/
                        ionbal.trimlo = SMALLFLOAT;
                        ionbal.trimhi = SMALLFLOAT;
                        lgEOL = false;
                }

                else
                {
                        /* set trim upper */
                        ionbal.trimhi = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));

                        /* set trim lower to same number */
                        ionbal.trimlo = ionbal.trimhi;
                }

                if( lgEOL )
                {
                        NoNumb(chCard);
                }

                if( ionbal.trimlo >= 1. || ionbal.trimhi >= 1. )
                {
                        fprintf( ioQQQ, " number must be negative since log\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        else if( nMatch("SKIP",chCard) )
        {
                /* skip punch command, for punching every n't point */
                i = 5;
                punch.ncPunchSkip = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
        }

        else if( nMatch(" UTA",chCard) )
        {
                /* set UTA command, to determine which UTA data set to use 
                 * default is to use Gu et al. 2006 data set */
                if( nMatch("GU06" , chCard ) )
                {
                        if( nMatch(" OFF" , chCard ) )
                        {
                                /* 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( nMatch("KISI" , chCard ) )
                {
                        if( nMatch(" OFF" , chCard ) )
                        {
                                /* 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( nMatch("EAKH",chCard) )
        {
                /* set WeakHeatCool, threshold on punch heating and cooling, default 0.05 */
                i = 5;
                punch.WeakHeatCool = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                if( lgEOL )
                {
                        NoNumb(chCard);
                }

                if( punch.WeakHeatCool < 0. )
                {
                        punch.WeakHeatCool = (realnum)pow((realnum)10.f,punch.WeakHeatCool);
                }
        }

        else if( nMatch("KSHE",chCard) )
        {
                /* upper limit to opacities for k-shell ionization */
                i = 5;
                continuum.EnergyKshell = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }

                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( nMatch("NCHR",chCard) )
        {
                /* option to set the number of charge states for grain model */
                long ii = 5;
                double val = FFmtRead(chCard,&ii,INPUT_LINE_LENGTH,&lgEOL);

                if( lgEOL )
                {
                        NoNumb(chCard);
                }
                else
                {
                        long nChrg = nint(val);
                        if( nChrg < 2 || nChrg > NCHS )
                        {
                                fprintf( ioQQQ, " illegal value for number of charge states: %ld\n", nChrg );
                                fprintf( ioQQQ, " choose a value between 2 and %d\n", NCHS );
                                fprintf( ioQQQ, " or increase NCHS in grainvar.h and recompile\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        else
                        {
                                SetNChrgStates(nChrg);
                        }
                }
        }

        else if( nMatch("NEGO",chCard) )
        {
                /* punch negative opacities if they occur, set negopac */
                opac.lgNegOpacIO = true;
        }

        else if( nMatch("NEND",chCard) )
        {
                /* 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 )
                {
                        i = 5;
                        /* this is default limit to number of zones, change it to this value */
                        geometry.nEndDflt = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        geometry.lgEndDflt = false;
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        /* now change all limits, for all iterations, to this value */
                        for( i=0; i < iterations.iter_malloc; i++ )
                        {
                                geometry.nend[i] = geometry.nEndDflt;
                        }
                }
        }

        else if( nMatch("TSQD",chCard) )
        {
                /* upper limit for highest density considered in the 
                 * Peimbert-style t^2 section of the printout */
                i = 5;
                peimbt.tsqden = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                if( lgEOL )
                {
                        NoNumb(chCard);
                }
                peimbt.tsqden = (realnum)pow((realnum)10.f,peimbt.tsqden);
        }

        else if( nMatch("NMAP",chCard) )
        {
                /* how many steps in plot or punch of heating-cooling map */
                i = 5;
                hcmap.nMapStep = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }
        }

        else if( nMatch("NUME",chCard) && nMatch("DERI",chCard) )
        {
                /* this is an option to use numerical derivatives for heating and cooling */
                NumDeriv.lgNumDeriv = true;
        }

        else if( nMatch("PATH",chCard) )
        {
                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( nMatch("PHFI",chCard) )
        {
                /* which version of PHFIT to use, 1995 or 1996 */
                i = 5;
                ip = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        NoNumb(chCard);
                }

                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 punch xxx command */
        else if( nMatch("PUNC",chCard) )
        {
                if( nMatch("HASH",chCard) )
                {
                        /* 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 PUNCH HASH line.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* specify the terminator between punch 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(punch.chHashString , "\n" );
                        }
                        else if( strcmp( chString_quotes_lowercase , "time" ) == 0 )
                        {
                                /* special case where output time between iterations */
                                sprintf(punch.chHashString , "TIME_DEP" );
                        }
                        else 
                        {
                                /* usual case, simply copy what is in quotes */
                                sprintf(punch.chHashString , chString_quotes_lowercase );
                        }
                }

                else if( nMatch("WIDT",chCard) )
                {
                        /* set punch line width for contrast in continuum plots */
                        i = 5;
                        /* units are km/sec */
                        punch.PunchLWidth = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                        /* lgEOL if no number hit, could have been c */
                        if( lgEOL )
                        {
                                if( nMatch(" C " , chCard ) )
                                {
                                        /* no number but _c_ entered, so enter speed of light in km/s */
                                        punch.PunchLWidth = (realnum)(SPEEDLIGHT/1e5);
                                }
                                else
                                {
                                        NoNumb(chCard);
                                }
                        }

                        if( punch.PunchLWidth <= 0. )
                        {
                                fprintf( ioQQQ, " line width must be greater than zero.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* >>chng 06 dec 06, value of PunchLWidth is km/s but SPEEDLIGHT
                         * is cm/s - add 1e5 before compare with SPEEDLIGHT */
                        else if( punch.PunchLWidth*0.9999e5 > SPEEDLIGHT )
                        {
                                fprintf( ioQQQ, " line width must be entered in km/s and be less than c.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* this is the factor that is used to add the lines to the continuum,
                        * 15 converts to cm/s */
                        punch.PunchLWidth = (realnum)(SPEEDLIGHT/(punch.PunchLWidth*1e5));

                }

                else if( nMatch("PREF",chCard) )
                {
                        /* specify a prefix before all punch filenames */
                        /* 
                        * 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 */

                        strcpy( punch.chFilenamePrefix , chString_quotes_lowercase );
                }

                else if( nMatch("FLUS",chCard) )
                {
                        /* flus the output after every iteration */
                        punch.lgFLUSH = true;
                }

                else
                {
                        fprintf( ioQQQ, " There should have been an option on this command.\n" );
                        fprintf( ioQQQ, " Valid options are HASH and PREFIX.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* set continuum .... options */
        else if( nMatch("CONT" , chCard ) )
        {
                if( nMatch("RESO" , chCard ) )
                {
                        /* set resolution, get factor that will multiply continuum resolution that
                        * is contained in the file continuum_mesh.ini */
                        i = 5;
                        continuum.ResolutionScaleFactor = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* negative numbers were logs */
                        if( continuum.ResolutionScaleFactor <= 0. )
                                continuum.ResolutionScaleFactor = pow(10.,continuum.ResolutionScaleFactor);
                }

                else if( nMatch("SHIE" , chCard ) )
                {
                        /* 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( nMatch("PESC" , chCard ) )
                        {
                                /* this uses an inward looking escape probability */
                                rt.nLineContShield = LINE_CONT_SHIELD_PESC;
                        }
                        else if( nMatch("FEDE" , chCard ) )
                        {
                                /* 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( nMatch("FERL" , chCard ) )
                        {
                                rt.nLineContShield = LINE_CONT_SHIELD_FERLAND;
                        }
                        else if( nMatch("NONE" , chCard ) )
                        {
                                /* 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;
}

---