/home66/gary/public_html/cloudy/c08_branch/source/parse_commands.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 */
/*ParseCommands main command line parser, decode command, then call other routines to read */
#include "cddefines.h"
#include "physconst.h"
#include "optimize.h"
#include "doppvel.h"
#include "stopcalc.h"
#include "abund.h"
#include "geometry.h"
#include "dense.h"
#include "plot.h"
#include "grid.h"
#include "rfield.h"
#include "grainvar.h"
#include "dynamics.h"
#include "magnetic.h"
#include "trace.h"
#include "h2.h"
#include "mole.h"
#include "hmi.h"
#include "rt.h"
#include "thermal.h"
#include "opacity.h"
#include "atomfeii.h"
#include "called.h"
#include "wind.h"
#include "hextra.h"
#include "iterations.h"
#include "radius.h"
#include "input.h" 
#include "assertresults.h"
#include "phycon.h"
#include "fudgec.h"
#include "version.h"
#include "conv.h"
#include "parse.h"

void ParseCommands(void)
{
        char chCard[INPUT_LINE_LENGTH], 
          chKey2[3], 
          chKey3[4], 
          chKey4[5], 
          chKey5[6];

        bool lgDSet, 
          lgEOF, 
          lgEOL, 
          lgNu2;
        bool lgStop ,
                lgStop_not_enough_info;

        long int i, 
          j, 
          nqh;
        long nInitFile=0;/* used to count number of init files read in */

        realnum a, 
                /* \todo        1       ar1 saves initial radius but does not appear
                 * to be used for anything.  all uses are on left side of =
                 * this could be removed across the code.  But first check
                 * that it is indeed not ever used */
          ar1, 
          teset, 
          z;
        double thickness;

        DEBUG_ENTRY( "ParseCommands()" );

        /* following says abundances are solar  */
        abund.lgAbnSolar = true;

        /* there are no plots desired yet */
        plotCom.lgPlotON = false;
        plotCom.nplot = 0;

        /* this flag remembers whether grains have ever been turned on.  It is passed
         * to routine ParseAbundances - there grains will be turned on with commands
         * such as abundances ism, unless grains were previously set 
         * with a grains command.  this way abundances will not clobber explicitly set
         * grains. */
        lgDSet = false;

        radius.Radius = 0.;
        radius.rinner = 0.;
        nqh = 0;
        rfield.nspec = 0;

        /* will be set true if no buffering command entered */
        input.lgSetNoBuffering = false;

        /* initialize some code variables in case assert command used in input stream */
        InitAssertResults();

        for( i=0; i < LIMSPC; i++ )
        {
                strcpy( rfield.chRSpec[i], "UNKN" );
        }
        optimize.nparm = 0;

        /* this is an option to turn on trace printout on the nth
         * call from the optimizer */
        if( optimize.lgTrOpt )
        {
                /* nTrOpt was set with "optimize trace" command,
                 * is iteration to turn on trace */
                if( optimize.nTrOpt == optimize.nOptimiz )
                {
                        trace.lgTrace = true;
                        called.lgTalk = true;
                        trace.lgTrOvrd = true;
                        fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
                }
        }

        /* say this is a beta version if we are talking and it is the truth */
        if( t_version::Inst().nBetaVer > 0 && called.lgTalk )
        {
                fprintf( ioQQQ, 
                         "\n                               This is a beta release of Cloudy, and is intended for testing only.\n\n" );
        }

        if( called.lgTalk )
        {
                /* this code prints pretty lines at top of output box */
                int indent = (int)((122 - strlen(t_version::Inst().chVersion))/2);
                fprintf( ioQQQ, "%*cCloudy %s\n", indent, ' ', t_version::Inst().chVersion );

                fprintf( ioQQQ, "%57cwww.nublado.org\n\n", ' ' );

                /* now print box and date of version, before printing commands */
                fprintf( ioQQQ, "%23c", ' ' );
                fprintf( ioQQQ, "**************************************");
                fprintf( ioQQQ, "%7.7s", t_version::Inst().chDate);
                fprintf( ioQQQ, "**************************************\n");

                fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
        }

        /* read in commands and print them   */

        /* following signals which read is in progress,
         * 1 is forward read, of input commands
         * -1 is reverse read from bottom of line array, of cloudy.ini file */
        input.iReadWay = 1;

        /* initialize array reader, this sub does nothing but set
         * initial value of a variable, depending on value of iReadWay */
        input_init();

        /* read until eof or blank line, then return control back to main program */

        /*
         * input_readarray puts the next stored line image into chCard
         * it will be <=INPUT_LINE_LENGTH char long, with eol at end
         * it will be in mixed upper and lower case
         * lgEOF is set true if we hit eof and no more lines present
         * code moving over to info contained in structures,
         *
         * input.chOrgCard == original image of command line
         * input.chCARDCAPS == original image converted to caps
         */
        input_readarray(chCard,&lgEOF);

        /* line starting with blank is eof, this checks we are not at eof */
        while( !lgEOF && chCard[0] != ' ' )
        {
                /* echo the line before we cap it, but only if it does
                 * not contain the keyword HIDE */
                /* >>chng 04 jan 21, add HIDE option, mostly for print quiet command */
                if( called.lgTalk && !nMatch("HIDE",input.chCARDCAPS) )
                        fprintf( ioQQQ, "%23c* %-80s*\n", ' ', chCard );

                /* change chCard to all caps */
                caps( chCard );

                /* now set up several partial keys */
                /* first two characters into chKey2 */
                strncpy( chKey2 , chCard , 2 );
                chKey2[2] = '\0';

                /* first three characters into chKey3 */
                strncpy( chKey3 , chCard , 3 );
                chKey3[3] = '\0';

                /* first four characters into chKey4 */
                strncpy( chKey4 , chCard , 4 );
                chKey4[4] = '\0';

                /* first four characters into chKey4 */
                strncpy( chKey5 , chCard , 5 );
                chKey5[5] = '\0';

                /* check whether VARY is on line */
                if( nMatch("VARY",chCard) )
                {
                        optimize.lgVarOn = true;
                        if( optimize.nparm + 1 > LIMPAR )
                        {
                                fprintf( ioQQQ, " Too many VARY lines entered; the limit is%4ld\n", 
                                  LIMPAR );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else
                {
                        optimize.lgVarOn = false;
                }

                /* look for "C " comment - this form of test should not be necessary since
                 * input line is padded with extra space - but works and would continue to
                 * work if padding were removed */
                if( chCard[0]=='C' && (chCard[1]==' ' || chCard[1]== 0) )
                {
                        /* this is null test since lines starting with "C " are comments,
                         * the char after the c can be a space or a newline */
                        i = 1;
                }

                /* start to look for keywords */
                else if( strcmp(chKey4,"ABSO") == 0 )
                {
                        /* enter luminosity in absolute magnitudes, in reads2 */
                        ParseAbsMag(chCard,&nqh);
                }

                else if( strcmp(chKey3,"AGE") == 0 )
                {
                        /* enter age of cloud so we can check for time-steady reads2 */
                        ParseAge(chCard);
                }

                else if( strcmp(chKey4,"AGN ") == 0 )
                {
                        /* enter generic style AGN continuum, in reads2 */
                        ParseAgn(chCard);
                }

                else if( strcmp(chKey4,"ABUN") == 0 )
                {
                        /* chemical abundances, readsun, lgDSet is set true if grains command
                         * ever entered, so that abundances command does not override grains command */
                        ParseAbundances(chCard,lgDSet);
                        /* abundances no longer solar */
                        abund.lgAbnSolar = false;
                }

                else if( strcmp(chKey4,"APER") == 0 )
                {
                        /* aperture command to simulate pencil beam or long slit */

                        /* if the "BEAM" or "SLIT" option is specified then only part 
                         * of the geometry is observed, and intensities
                         * should not be weighted by r^2.  There are two limiting cases, SLIT in which
                         * the slit is longer than the diameter of the nebula and the contribution to the
                         * detected luminosity goes as r^1, and BEAM when the contribution is r^0, 
                         * the same as plane parallel 
                         * 
                         * default value of geometry.iEmissPower is 2 (set in zero.c) for full geometry  
                         */
                        if( nMatch("SLIT",chCard) )
                        {
                                /* long slit is case where slit is longer than diameter, so emissivity contributes
                                 * r^1 to the observed luminosity */
                                geometry.iEmissPower = 1;
                        }
                        else if( nMatch("BEAM",chCard) )
                        {
                                /* long slit is case where slit is longer than diameter, so emissivity contributes
                                 * r^1 to the observed luminosity */
                                /* this is the default or SHORT case, where r^0 is dependence */
                                geometry.iEmissPower = 0;
                        }
                        else
                        {
                                fprintf( ioQQQ, " One of the keywords SLIT or BEAM must appear.\n" );
                                fprintf( ioQQQ, " Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( strcmp(chKey4,"ASSE") == 0 )
                {
                        /* assert that code predicts certain results, in assertresults.h */
                        ParseAssertResults();
                }

                else if( strcmp(chKey4,"ATOM") == 0 )
                {
                        /* accept both forms of feii */
                        if( nMatch("FEII",chCard) || nMatch("FE II",chCard) )
                        {
                                /* parse the atom FeII command - routine in atom_feii.c */
                                ParseAtomFeII(chCard);
                        }

                        else if( nMatch("H-LI",chCard) )
                        {
                                /* parse the atom h-like command */
                                ParseAtomISO(ipH_LIKE, chCard);
                        }

                        else if( nMatch("HE-L",chCard) )
                        {
                                /* parse the atom he-like command */
                                ParseAtomISO(ipHE_LIKE, chCard);
                        }

                        else if( nMatch("ROTO",chCard) )
                        {
                                /* ATOM ROTOR same as ATOM CO */
                                fprintf(ioQQQ,"PROBLEM - the atom rotor command is now the ATOM CO command.  "
                                        "Please use that instead.  I will accept this for now but may not in future versions.\n");
                                ParseAtomCO(chCard);
                        }

                        else if( nMatch(" CO ",chCard) )
                        {
                                /* parameters for the rigid rotor CO molecules (yes, its not an atom)
                                 * command is ATOM CO */
                                ParseAtomCO(chCard);
                        }

                        else if( nMatch(" H2 ",chCard) )
                        {
                                /* parameters for the rigid rotor CO molecules (yes, its not an atom)
                                 * command is ATOM ROTOR, routine is stored in parse_atomh2 */
                                ParseAtomH2(chCard);
                        }

                        else
                        {
                                fprintf( ioQQQ, " I could not recognize a keyword on this atom command.\n");
                                fprintf( ioQQQ, " The available keys are FeII, H-Like, He-like, rotor and H2.\n");
                                fprintf( ioQQQ, " Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( strcmp(chKey4,"BACK") == 0 )
                {
                        /* cosmic background, in parse_backgrd */
                        ParseBackgrd(&nqh,chCard,&ar1);
                }

                else if( strcmp(chKey4,"BLAC") == 0 )
                {
                        /* black body, in reads2 */
                        ParseBlackbody(chCard,&nqh,&ar1);

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* no luminosity options on vary */
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "BLACKbody=%f" );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /* log of temp stored here  */
                                /* >>chng 02 feb 11, log was not here, don't know what happened to it */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.slope[rfield.nspec-1]);
                                /* the increment in the first steps away from the original value */
                                optimize.vincr[optimize.nparm] = 0.5;
                                optimize.nparm += 1;
                        }
                }

                else if( strcmp(chKey4,"BREM") == 0 )
                {
                        /* bremsstrahlung continuum from central object */
                        strcpy( rfield.chSpType[rfield.nspec], "BREMS" );
                        i = 5;
                        rfield.slope[rfield.nspec] = 
                                (realnum)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        /* temperature is interpreted as log if <=10, linear otherwise*/
                        if( rfield.slope[rfield.nspec] <= 10. )
                        {
                                rfield.slope[rfield.nspec] =  
                                        pow(10.,rfield.slope[rfield.nspec]);
                        }
                        rfield.cutoff[rfield.nspec][0] = 0.;

                        /* option for vary keyword */
                        if( optimize.lgVarOn )
                        {
                                /* only one parameter */
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "BREMS, T=%f" );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /* log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] =  (realnum)rfield.slope[rfield.nspec];
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                        ++rfield.nspec;
                        if( rfield.nspec >= LIMSPC )
                        {
                                /* too many continua were entered */
                                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( strcmp(chKey4,"CASE") == 0 )
                {
                        /* do hydrogen case b */
                        ParseCaseB( chCard );
                }

                else if( strcmp(chKey4,"CEXT") == 0 )
                {
                        /* add "extra" cooling, power law temp dependence */
                        thermal.lgCExtraOn = true;
                        i = 5;
                        thermal.CoolExtra = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        thermal.cextpw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                }

                /* replace with CMB command, both will continue to be parsed */
                else if( (strcmp(chKey3,"CMB") == 0) || (strcmp(chKey4,"FIRE") == 0) )
                {
                        /* cosmic thermal background radiation, argument is redshift */
                        i = 5;
                        /* if no number on line then (realnum)FFmtRead returns z=0; i.e., now */
                        z = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        /* in readsun */
                        ParseCMB(z,&nqh,&ar1);
                }

                else if( strcmp(chKey4,"COMP") == 0 )
                {
                        /* compile ascii version of stellar atmosphere continua in volk */
                        ParseCompile(chCard);
                }

                else if( strcmp(chKey4,"CONS") == 0 )
                {
                        /* constant temperature, pressure, density, or gas pressure
                         * in readsun */
                        ParseConstant(chCard);
                }

                else if( strcmp(chKey4,"CORO") == 0 )
                {
                        /* coronal equilibrium; set constant temperature to number on line
                         *  in readsun */
                        ParseCoronal( chCard,&nqh,&ar1);
                }

                else if( strcmp(chKey4,"COSM") == 0 )
                {
                        /* cosmic ray density, log of rate relative to background, or density of rel. electrons,
                         * quantity is log unless keyword linear appears */
                        ParseCosmicRays( chCard );
                }

                else if( strcmp(chKey4,"COVE") == 0 )
                {
                        /* covering factor for gas */
                        i = 5;
                        /* The geometric covering factor accounts for how much of 4\pi is
                         * covered by gas, and so linearly multiplies the predicted intensities */
                        geometry.covgeo = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        /* if negative then log, convert to linear */
                        if( geometry.covgeo <= 0. )
                        {
                                geometry.covgeo = (realnum)pow((realnum)10.f,geometry.covgeo);
                        }

                        if( geometry.covgeo > 1. )
                        {
                                fprintf( ioQQQ, " A covering factor greater than 1 makes no physical sense.  Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* radiative transfer covering factor will be equal to the geometric one */
                        geometry.covrt = geometry.covgeo;
                }

                else if( strcmp(chKey4,"CRAS") == 0 )
                {
                        /* any of several tests to cause the code to crash */
                        ParseCrashDo(chCard);
                }

                else if( strcmp(chKey4,"CYLI") == 0 )
                {
                        /* height of cylinder in cm */
                        i = 5;
                        radius.lgCylnOn = true;
                        radius.CylindHigh = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                }

                else if( strcmp(chKey4,"DIEL") == 0 )
                {
                        /* >>chng 05 dec 21, change dielectronic command to set dielectronic recombination command */
                        fprintf( ioQQQ, " The DIELectronic command has been replaced with the SET DIELectronic recombination command.\n" );
                        fprintf( ioQQQ, " Please have a look at Hazy.\n Sorry.\n\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                else if( strcmp(chKey4,"DIFF") == 0 )
                {
                        /* set method of transferring diffuse fields,
                         * default is outward only, cdDffTrns label "OU2", set in zero.c */
                        if( nMatch(" OTS",chCard) )
                        {
                                if( nMatch("SIMP",chCard) )
                                {
                                        /* this is simple ots, which never allows photons to escape */
                                        strcpy( rfield.chDffTrns, "OSS" );
                                }
                                else
                                {
                                        /* this is regular ots, which allows photons to escape */
                                        strcpy( rfield.chDffTrns, "OTS" );
                                }
                                rfield.lgOutOnly = false;
                        }
                        else if( nMatch(" OUT",chCard) )
                        {
                                rfield.lgOutOnly = true;
                                i = 4;
                                j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                if( lgEOL )
                                {
                                        /* this is the default set in zero */
                                        strcpy( rfield.chDffTrns, "OU2" );
                                }
                                else
                                {
                                        if( j > 0 && j < 10 )
                                        {
                                                sprintf( rfield.chDffTrns, "OU%1ld", j );
                                        }
                                        else
                                        {
                                                fprintf( ioQQQ, " must be between 1 and 9 \n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }

                        else
                        {
                                fprintf( ioQQQ, " There should have been OUTward or OTS on this line.  Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( strcmp(chKey4,"DIST") == 0 )
                {
                        /* distance to the object */
                        i = 5;
                        radius.distance = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* default is for quantity to be log of distance, linear keyword
                         * overrides this - is LINEAR is not on line then exp */
                        if( !nMatch("LINE",chCard ) )
                        {
                                radius.distance = pow(10., radius.distance );
                        }
                        /* default is radius in cm - if parsecs appears then must 
                         * convert to cm */
                        if( nMatch("PARS",chCard ) )
                        {
                                radius.distance *= PARSEC;
                        }
                }

                else if( strcmp(chKey4,"DLAW") == 0 )
                {
                        /* either use dense_fabden routine, or read in table of points */
                        ParseDLaw(chCard);
                }

                else if( strcmp(chKey4,"DOUB") == 0 )
                {
                        /* double optical depth scale after each iteration */
                        rt.DoubleTau = 2.;
                }

                else if( strcmp(chKey4,"DRIV") == 0 )
                {
                        /* call one of several drivers, readsun */
                        ParseDrive(chCard);
                }

                else if( strcmp(chKey4,"EDEN") == 0 )
                {
                        /* option to add "extra" electrons, to test Compton limit
                         *  for very low T(star) - option is log of eden */
                        i = 5;
                        dense.EdenExtra = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* warn that this model is meaningless */
                        phycon.lgPhysOK = false;
                }

                else if( strcmp(chKey4,"ELEM") == 0 )
                {
                        /* element command;
                         * scale or abundance options, to change abundance of specific element
                         * read option to change order of elements
                         * in reads2.f */
                        ParseElement(chCard);
                }

                else if( strcmp(chKey4,"ENER") == 0 )
                {
                        /* energy density (degrees K) of this continuum source */
                        i = 5;
                        if( nqh >= LIMSPC )
                        {
                                /* too many continua were entered */
                                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* silly, but calms down the lint */
                        ASSERT( nqh < LIMSPC );
                        strcpy( rfield.chRSpec[nqh], "SQCM" );
                        teset = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* set radius to very large value if not already set */
                        /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
                        if( !radius.lgRadiusKnown )
                        {
                                /* set radius to large value */
                                ar1 = (realnum)radius.rdfalt;
                                radius.Radius = pow(10.,radius.rdfalt);
                        }

                        /* convert temp to log, recognize linear option */
                        if( nMatch("LINE",chCard) || teset > 10. )
                        {
                                /* option for linear temperature, must store log */
                                teset = (realnum)log10(teset);
                        }

                        if( teset > 5. )
                        {
                                fprintf( ioQQQ, " This intensity may be too large.  The code may crash due to overflow.  Was log intended?\n" );
                        }

                        /* teset is not log of temp, now get log of total luminosity */
                        strcpy( rfield.chSpNorm[nqh], "LUMI" );

                        /* full range of continuum will be used */
                        rfield.range[nqh][0] = rfield.emm;
                        rfield.range[nqh][1] = rfield.egamry;
                        rfield.totpow[nqh] = (4.*teset - 4.2464476 + 0.60206);

                        /* >>chng 06 mar 22, add time option to vary only some continua with time */
                        if( nMatch( "TIME" , chCard ) )
                                rfield.lgTimeVary[nqh] = true;

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                strcpy( optimize.chVarFmt[optimize.nparm], "ENERGY DENSITY %f log " );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.totpow[nqh]);
                                optimize.vincr[optimize.nparm] = 0.1f;
                                optimize.nvarxt[optimize.nparm] = 1;
                                optimize.nparm += 1;
                        }

                        /* finally increment number of continua */
                        ++nqh;
                }

                else if( strcmp(chKey4,"EXTI") == 0 )
                {
                        /* extinguish ionizing continuum by absorbing column AFTER
                         * setting luminosity or Q(H).  First number is the column
                         * density (log), second number is leakage (def=0%)
                         * last number is lowest energy (ryd), last two may be omitted
                         * from right to left 
                         * 
                         * extinction is actually done in extin, which is called by ContSetIntensity */
                        ParseExtinguish( chCard );
                }

                else if( strcmp(chKey4,"FAIL") == 0 )
                {
                        /* reset number of temp failures allowed, default=20 */
                        i = 5;

                        /* save previous value */
                        j = conv.LimFail;
                        conv.LimFail = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        /* >>chng 01 mar 14, switch logic on maps */
                        /* ' map' flag, make map when failure, default is no map,
                         * second check is so that no map does not trigger a map */
                        if( nMatch(" MAP",chCard) && !nMatch(" NO ",chCard) )
                        {
                                conv.lgMap = true;
                        }

                        /* complain if failures was increased above default */
                        if( conv.LimFail > j )
                        {
                                fprintf( ioQQQ, " This command should not be necessary.\n" );
                                fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this simulation.\n" );
                        }
                }

                else if( strcmp(chKey4,"FILL") == 0 )
                {
                        /* filling factor, power law exponent (default=1., 0.) */
                        i = 5;
                        a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        if( a <= 0. )
                        {
                                /* number less than or equal to 0, must have been entered as log */
                                geometry.FillFac = (realnum)pow((realnum)10.f,a);
                        }
                        else
                        {
                                /* number greater than zero, must have been the real thing */
                                geometry.FillFac = a;
                                if( geometry.FillFac > 1.0 )
                                {
                                        if( called.lgTalk )
                                        {
                                                fprintf( ioQQQ, " Filling factor > 1, reset to 1\n" );
                                        }
                                        geometry.FillFac = 1.;
                                }
                        }
                        geometry.fiscal = geometry.FillFac;

                        /* option to have power law dependence, filpow set to 0 in zerologic */
                        geometry.filpow = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                strcpy( optimize.chVarFmt[optimize.nparm], "FILLING FACTOR= %f power=%f" );

                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(geometry.FillFac);
                                optimize.vincr[optimize.nparm] = 0.5;

                                /* power law dependence here, but cannot be varied */
                                optimize.vparm[1][optimize.nparm] = geometry.filpow;
                                optimize.nvarxt[optimize.nparm] = 2;

                                /* do not allow filling factor to go positive */
                                optimize.varang[optimize.nparm][0] = -1e38f;
                                optimize.varang[optimize.nparm][1] = 0.f;
                                optimize.nparm += 1;
                        }
                }

                else if( strcmp(chKey4,"FLUC") == 0 )
                {
                        /* rapid density fluctuations, in readsun */
                        ParseFluc(chCard);
                }

                else if( strcmp(chKey4,"F(NU") == 0 )
                {
                        /* this is the specific flux at nu
                         *  following says F(nu) not nuF(nu) */
                        lgNu2 = false;
                        /* in reads2 */
                        ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2);
                }

                else if( strcmp(chKey4,"FORC") == 0 )
                {
                        /* force temperature of first zone, but don't keep constant
                         * allow to then go to nearest equilibrium
                         *  log if < 10 */
                        i = 5;
                        thermal.ConstTemp = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        if( nMatch(" LOG",chCard) || (thermal.ConstTemp <= 10. && 
                          !nMatch("LINE",chCard)) )
                        {
                                thermal.ConstTemp = (realnum)pow((realnum)10.f,thermal.ConstTemp);
                        }

                        /* low energy bounds of continuum array do not permit too-low a temp */
                        if( thermal.ConstTemp < 3. )
                        {
                                fprintf( ioQQQ, " TE reset to 3K: entered number too small.\n" );
                                thermal.ConstTemp = 3.;
                        }
                }

                else if( strcmp(chKey4,"FUDG") == 0 )
                {
                        /* enter a fudge factor */
                        i = 5;
                        for( j=0; j < NFUDGC; j++ )
                        {
                                fudgec.fudgea[j] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                /* fudgec.nfudge is number of factors on the line */
                                if( !lgEOL )
                                        fudgec.nfudge = j+1; 
                        }

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* no luminosity options on vary */
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "FUDGE=%f" );
                                /* enter remaining parameters as constants */
                                char chHold[1000];
                                for( j=1; j<fudgec.nfudge; ++j )
                                {
                                        sprintf( chHold , " %f" , fudgec.fudgea[j] );
                                        strcat( optimize.chVarFmt[optimize.nparm] , chHold );
                                }

                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /* fudge factor stored here  */
                                optimize.vparm[0][optimize.nparm] = fudgec.fudgea[0];
                                /* the increment in the first steps away from the original value */
                                optimize.vincr[optimize.nparm] = 0.5;
                                optimize.nparm += 1;
                        }
                }

                else if( strcmp(chKey4,"GLOB") == 0 )
                {
                        /* globule with density increasing inward
                         * parameters are outer density, radius of globule, and density power */
                        ParseGlobule(chCard);
                }

                else if( (strcmp(chKey4,"GRAI") == 0 )||( strcmp(chKey4,"PGRA") == 0 ) )
                {
                        /* read parameters dealing with grains, in reads2 */
                        ParseGrain(chCard,&lgDSet);
                }

                else if( strcmp(chKey4,"GRID") == 0 )
                {
                        /* option to run grid of models by varying certain parameters
                         * in reads2 */
                        ParseGrid(chCard);
                }

                else if( strcmp(chKey4,"HDEN") == 0 )
                {
                        /* parse the hden command to set the hydrogen density, in reads2 */
                        ParseHDEN(chCard);
                }

                else if( strcmp(chKey4,"HELI") == 0 )
                {
                        fprintf(ioQQQ,"Sorry, this command is replaced with ATOM HE-LIKE\n");
                        cdEXIT(EXIT_FAILURE);
                }

                else if( strcmp(chKey4,"HEXT") == 0 )
                {
                        /* "extra" heating rate, so that first= log(erg(cm-3, s-1))
                         * second optional number is scale radius, so that HXTOT = TurbHeat*SEXP(DEPTH/SCALE)
                         * if missing then constant heating.
                         * third option is depth from shielded face, to mimic irradiation from both sides*/
                        i = 5;
                        hextra.TurbHeat = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                        if( lgEOL )
                                NoNumb( chCard );
                        /* save initial value in case reset with time dependent command */
                        hextra.TurbHeatSave = hextra.TurbHeat;

                        /* remember which type of scale dependence we will use */
                        const char *chHextraScale;
                        /* these are optional numbers on depth or density option */
                        realnum par1=0. , par2=0.;
                        long int nPar;
                        if( nMatch("DEPT" , chCard ) )
                        {
                                /* depend on depth */
                                hextra.lgHextraDepth = true;
                                chHextraScale = "DEPTH";
                                /* optional scale radius */
                                a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                if( lgEOL )
                                        NoNumb(chCard);
                                else
                                        hextra.turrad = (realnum)pow((realnum)10.f,a);

                                /* depth from shielded face, to mimic illumination from both sides 
                                 * may not be specified */
                                realnum b = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                if( lgEOL )
                                {
                                        hextra.turback = 0.;
                                        nPar = 2;
                                }
                                else
                                {
                                        hextra.turback = (realnum)pow((realnum)10.f,b);
                                        nPar = 3;
                                }
                                par1 = a;
                                par2 = b;
                        }
                        else if( nMatch("DENS" , chCard ) )
                        {
                                /* depends on density */
                                chHextraScale = "DENSITY";
                                hextra.lgHextraDensity = true;

                                /* the log of the density */
                                a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                /* if number not entered then unity is used, heating 
                                 * proportional to density */
                                hextra.HextraScaleDensity = (realnum)pow((realnum)10.f,a);
                                par1 = a;
                                par2 = 0;
                                nPar = 2;
                        }
                        else
                        {
                                /* constant heating */
                                chHextraScale = "";
                                nPar = 1;
                        }

                        /* option to have heating vary with time */
                        if( nMatch( "TIME" , chCard ) )
                                hextra.lgTurbHeatVaryTime = true;

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* 1 to 3 options on hextra vary */
                                optimize.nvarxt[optimize.nparm] = nPar;
                                optimize.vparm[0][optimize.nparm] = log10(hextra.TurbHeat);
                                optimize.vparm[1][optimize.nparm] = par1;
                                optimize.vparm[2][optimize.nparm] = par2;

                                strcpy( optimize.chVarFmt[optimize.nparm], "HEXTra %f " );
                                strcat( optimize.chVarFmt[optimize.nparm], chHextraScale );
                                while( nPar > 1 )
                                {
                                        /* add extra spots to write parameters */
                                        --nPar;
                                        strcat( optimize.chVarFmt[optimize.nparm], " %f" );
                                }
                                if( hextra.lgTurbHeatVaryTime )
                                        strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /* the increment in the first steps away from the original value */
                                optimize.vincr[optimize.nparm] = 0.1f;
                                optimize.nparm += 1;
                        }
                }

                else if( strcmp(chKey4,"HIGH") == 0 )
                {
                        /* approach equilibrium from high te */
                        thermal.lgTeHigh = true;
                }

                else if( strncmp( chCard ,"HYDROGEN",8) == 0 )
                {
                        fprintf(ioQQQ," Sorry, this command has been replaced with the ATOM H-LIKE command.\n");
                        cdEXIT(EXIT_FAILURE);
                }

                else if( strcmp(chKey4,"ILLU") == 0 )
                {
                        /* illumination angle */
                        i = 5;
                        geometry.AngleIllumRadian = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* default is degrees, but radian key also there */
                        if( nMatch("RADI",chCard) )
                        {
                                /* entered as radians, convert to degrees first */
                                geometry.AngleIllumRadian /= (realnum)(PI/180.);
                        }
                        /* we now have degrees, make sure between 0 and 90 */
                        if( geometry.AngleIllumRadian < 0. || geometry.AngleIllumRadian >= 90. )
                        {
                                fprintf( ioQQQ, " Angle of illumination must be between 0 and 90 degrees "
                                        "or 0 and pi/2 radians.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* we really want 1. / COS( theta ) with theta in radians 
                         * first convert angle to radians */
                        geometry.AngleIllumRadian = (realnum)(geometry.AngleIllumRadian*PI/180.);

                        /* this is effective path - dTau_eff = dTau_normal * geometry.DirectionalCosin */
                        geometry.DirectionalCosin = (realnum)(1./cos(geometry.AngleIllumRadian));

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* no luminosity options on vary */
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "ILLUminate %f radians " );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                optimize.vparm[0][optimize.nparm] = geometry.AngleIllumRadian;
                                /* the increment in the first steps away from the original value */
                                optimize.vincr[optimize.nparm] = 0.1f;
                                optimize.nparm += 1;
                        }
                }

                else if( strcmp(chKey4,"INIT") == 0 )
                {
                        /* read cloudy.ini initialization file
                         * need to read in file every time, since some vars
                     * get reset in zero - would be unsafe not to read in again */
                        ParseInit(chCard);

                        /* check that only one init file was in the input stream -
                         * we cannot now read more than one */
                        ++nInitFile;
                        if( nInitFile > 1 )
                        {
                                fprintf( ioQQQ, 
                                        " This is the second init file, I can only handle one.\nSorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* we will put the data for the ini file at the end of the array of
                         * line images, this is the increment for stuffing the lines in - negative */
                        input.iReadWay = -1;

                        /* initialize array reader, this sub does nothing but set
                         * initial value of a variable, depending on value of iReadWay */
                        input_init();
                }

                else if( strcmp(chKey5,"INTEN") == 0 )
                {
                        /* intensity of this continuum source */
                        i = 5;
                        if( nqh >= LIMSPC )
                        {
                                /* too many continua were entered */
                                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* silly, but calms down the lint */
                        ASSERT( nqh < LIMSPC );
                        strcpy( rfield.chRSpec[nqh], "SQCM" );
                        rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        /* set radius to very large value if not already set */
                        /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
                        if( !radius.lgRadiusKnown )
                        {
                                /* set radius to large value */
                                ar1 = (realnum)radius.rdfalt;
                                radius.Radius = pow(10.,radius.rdfalt);
                        }
                        if( nMatch("LINE",chCard) )
                        {
                                /* silly, but calms down the lint */
                                ASSERT( nqh < LIMSPC );
                                /* option for linear input parameter */
                                rfield.totpow[nqh] = log10(rfield.totpow[nqh]);
                        }
                        strcpy( rfield.chSpNorm[nqh], "LUMI" );
                        /* ParseRangeOption in readsun */
                        ParseRangeOption(nqh,chCard);

                        /* >>chng 06 mar 22, add time option to vary only some continua with time */
                        if( nMatch( "TIME" , chCard ) )
                                rfield.lgTimeVary[nqh] = true;

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* create the format we will use to write the command */
                                /* >>chng 03 apr 30, had used log for range option, but if one or other
                                 * was 1 ryd, then became 0 in the log, which was misinterpreted as
                                 * one of the continuum bounds */
                                /*strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY%f log range %f %f" );*/
                                strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY %f range %f %f " );
                                /* array index for this command */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh];
                                optimize.vincr[optimize.nparm] = 0.5;
                                /* range option, but cannot be varied */
                                /*optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[nqh][0]);
                                optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[nqh][1]);*/
                                /* >>change to linear option */
                                optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
                                optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
                                optimize.nvarxt[optimize.nparm] = 3;
                                ++optimize.nparm;
                        }
                        ++nqh;
                }

                else if( strcmp(chKey5,"INTER") == 0 )
                {
                        /* interpolate on input tables for continuum, set of power  laws used
                         * input ordered pairs nu( ryd or log(Hz) ), log(fnu)
                         * additional lines begin CONTINUE
                         * first check that this is the one and only INTERP command
                         * in readsun */
                        ParseInterp(chCard,&lgEOF);
                }

                else if( strcmp(chKey4,"IONI") == 0 )
                {
                        /* inter ionization parameter U=Q/12 R*R N C;
                         * defined per hydrogen, not per electron (as before)
                         * radius must also be entered if spherical, but not needed if plane */
                        ParseIonPar(&nqh,chCard,&ar1);
                }

                else if( strcmp(chKey4,"ITER") == 0 )
                {
                        /* iterate to get optical depths and diffuse field properly */
                        i = 5;
                        iterations.itermx = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) - 1;
                        iterations.itermx = MAX2(iterations.itermx,1);
                        /* >>chng 06 mar 19, malloc itrdim arrays 
                         * iterations.iter_malloc is -1 before space allocated and set to
                         * itermx if not previously set */
                        if( iterations.itermx > iterations.iter_malloc - 1 )
                        {
                                long int iter_malloc_save = iterations.iter_malloc;
                                /* add one for mismatch between array dim and iterations defn,
                                 * another few for safety */
                                iterations.iter_malloc = iterations.itermx+3;
                                iterations.IterPrnt = (long int*)REALLOC(iterations.IterPrnt ,
                                        (size_t)iterations.iter_malloc*sizeof(long int) );
                                geometry.nend = (long int*)REALLOC(geometry.nend ,
                                        (size_t)iterations.iter_malloc*sizeof(long int) );
                                radius.router = (double*)REALLOC(radius.router ,
                                        (size_t)iterations.iter_malloc*sizeof(double) );
                                /* >>chng 06 jun 07, need to set IterPrnt, router, and nend */
                                for(j=iter_malloc_save; j<iterations.iter_malloc; ++j )
                                {
                                        iterations.IterPrnt[j] = iterations.IterPrnt[iter_malloc_save-1];
                                        geometry.nend[j] = geometry.nend[iter_malloc_save-1];
                                        radius.router[j] = radius.router[iter_malloc_save-1];
                                }
                        }

                        if( nMatch("CONV",chCard) )
                        {
                                /* option to keep iterating until it converges, checks on convergence
                                 * are done in update, and checked again in prtcomment */
                                conv.lgAutoIt = true;
                                /* above would have been legitimate setting of ITERMX, set to default 10 */
                                if( lgEOL )
                                {
                                        iterations.itermx = 10 - 1;
                                }
                                a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                /* change convergence criteria, preset in zero */
                                if( !lgEOL )
                                {
                                        conv.autocv = a;
                                }
                        }
                }

                else if( strcmp(chKey4,"L(NU") == 0 )
                {
                        /* this is the specific luminosity at nu
                         * following says n(nu) not nuF(nu) */
                        lgNu2 = false;
                        /* in reads2 */
                        ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2);
                }

                else if( strcmp(chKey4,"LASE") == 0 )
                {
                        /* mostly one frequency (a laser) to check gamma's,*/
                        /* say the continuum type */
                        strcpy( rfield.chSpType[rfield.nspec], "LASER" );

                        /* scan off the laser's energy */
                        i = 5;
                        rfield.slope[rfield.nspec] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);

                        /* negative energies are logs */
                        if( rfield.slope[rfield.nspec] <= 0. )
                        {
                                rfield.slope[rfield.nspec] = 
                                pow(10.,rfield.slope[rfield.nspec]);
                        }
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }

                        /* next number is optional relative width of laser */
                        rfield.cutoff[rfield.nspec][0] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                /* default width is 0.05 */
                                rfield.cutoff[rfield.nspec][0] = 0.05;
                        }

                        /* increment counter making sure we still have space enough */
                        ++rfield.nspec;
                        if( rfield.nspec >= LIMSPC )
                        {
                                /* too many continua were entered */
                                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( strcmp(chKey4,"LUMI") == 0 )
                {
                        /* luminosity of this continuum source */
                        i = 5;
                        if( nqh >= LIMSPC )
                        {
                                /* too many continua were entered */
                                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        strcpy( rfield.chRSpec[nqh], "4 PI" );
                        rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        if( nMatch("LINE",chCard) )
                        {
                                /* option for linear input parameter */
                                rfield.totpow[nqh] = log10(rfield.totpow[nqh]);
                        }

                        strcpy( rfield.chSpNorm[nqh], "LUMI" );

                        /* the solar option - number is log total solar luminosity */
                        if( nMatch("SOLA",chCard) )
                        {
                                /* option to use log of total luminosity in solar units */
                                rfield.range[nqh][0] = rfield.emm;
                                rfield.range[nqh][1] = rfield.egamry;
                                /* log of solar luminosity in erg s-1 */
                                rfield.totpow[nqh] += 33.5827f;
                        }
                        else
                        {
                                /* check if range is present and parse it if it is - ParseRangeOption in readsun 
                                 * this includes TOTAL keyword for total luminosity */
                                ParseRangeOption(nqh,chCard);
                        }

                        /* >>chng 06 mar 22, add time option to vary only some continua with time */
                        if( nMatch( "TIME" , chCard ) )
                                rfield.lgTimeVary[nqh] = true;

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* >>chng 03 apr 30, had used log for range option, but if one or other
                                 * was 1 ryd, then became 0 in the log, which was misinterpreted as
                                 * one of the continuum bounds */
                                strcpy( optimize.chVarFmt[optimize.nparm], "LUMINOSITY %f range %f %f " );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh];
                                optimize.vincr[optimize.nparm] = 0.5;
                                /* range option in, but cannot be varied */
                                optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
                                optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
                                optimize.nvarxt[optimize.nparm] = 3;
                                optimize.nparm += 1;
                        }
                        ++nqh;
                }

                else if( strcmp(chKey4,"MAGN") == 0 )
                {
                        /* parse the magnetic field command, routine in magnetic.c */
                        ParseMagnet( chCard );
                }

                else if( strcmp(chKey4,"MAP ") == 0 )
                {
                        /* do cooling space map for specified zones, in reads2 */
                        ParseMap(chCard);
                }

                else if( strcmp(chKey4,"META") == 0 )
                {
                        /* read depletion for metals, all elements heavier than He
                         * in reads2 */
                        ParseMetal(chCard);
                        /* abundances no longer solar */
                        abund.lgAbnSolar = false;
                }

                else if( strcmp(chKey4,"NEUT") == 0 )
                {
                        /* heating and ionization due to fast neutrons */
                        hextra.lgNeutrnHeatOn = true;

                        /* first number is neutron luminosity
                         * relative to bolometric luminosity */
                        i = 5;
                        hextra.frcneu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                                NoNumb(chCard);

                        /* save as log of fraction */
                        if( hextra.frcneu > 0. )
                                hextra.frcneu = (realnum)log10(hextra.frcneu);

                        /* second number is efficiency */
                        hextra.effneu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                hextra.effneu = 1.0;
                        }
                        else
                        {
                                if( hextra.effneu <= 0. )
                                        hextra.effneu = (realnum)pow((realnum)10.f,hextra.effneu);
                        }
                }

                else if( strcmp(chKey3,"NO ") == 0 )
                {
                        /* don't do something, in readsun */
                        ParseDont(chCard);
                }

                else if( strcmp(chKey4,"NORM") == 0 )
                {
                        /* normalize lines to this rather than h-b, sec number is scale factor */
                        ParseNorm(chCard);
                }

                else if( strcmp(chKey4,"NUF(") == 0 )
                {
                        /* flux density of this continuum source, at optional frequency
                         *  expressed as product nu*f_nu */
                        lgNu2 = true;
                        /* in reads2 */
                        ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2);
                }

                else if( strcmp(chKey4,"NUL(") == 0 )
                {
                        /* specific luminosity density of this continuum source, at opt freq
                         * expressed as product nu*f_nu */
                        lgNu2 = true;
                        /* in reads2 */
                        ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2);
                }

                else if( strcmp(chKey4,"OPTI") == 0 )
                {
                        /* option to optimize the model by varying certain parameters
                         * in reads2 */
                        ParseOptimize(chCard);
                }

                else if( strcmp(chKey4,"PHI(") == 0 )
                {
                        /* enter phi(h), the number of h-ionizing photons/cm2 */
                        i = 5;
                        if( nqh >= LIMSPC )
                        {
                                /* too many continua were entered */
                                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* silly, but calms down the lint */
                        ASSERT( nqh < LIMSPC );
                        strcpy( rfield.chRSpec[nqh], "SQCM" );
                        strcpy( rfield.chSpNorm[nqh], "PHI " );
                        rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* set radius to very large value if not already set */
                        /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
                        if( !radius.lgRadiusKnown )
                        {
                                /* set radius to large value */
                                ar1 = (realnum)radius.rdfalt;
                                radius.Radius = pow(10.,radius.rdfalt);
                        }
                        /* check if continuum so intense that crash is likely */
                        if( rfield.totpow[nqh] > 29. )
                        {
                                fprintf( ioQQQ, " Is the flux for this continuum correct?\n" );
                                fprintf( ioQQQ, " It appears too bright to me.\n" );
                        }
                        /* ParseRangeOption in readsun xx parse_rangeoption*/
                        ParseRangeOption(nqh,chCard);

                        /* >>chng 06 mar 22, add time option to vary only some continua with time */
                        if( nMatch( "TIME" , chCard ) )
                                rfield.lgTimeVary[nqh] = true;

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* >>chng 03 apr 30, had used log for range option, but if one or other
                                 * was 1 ryd, then became 0 in the log, which was misinterpreted as
                                 * one of the continuum bounds */
                                strcpy( optimize.chVarFmt[optimize.nparm], "phi(h) %f range %f %f" );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh];
                                optimize.vincr[optimize.nparm] = 0.5;
                                /* range option in, but cannot be varied */
                                optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
                                optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
                                optimize.nvarxt[optimize.nparm] = 3;
                                optimize.nparm += 1;
                        }
                        ++nqh;
                }

                else if( strcmp(chKey4,"PLOT") == 0 )
                {
                        /* make plot of log(nu.f(nu)) vs log(nu) or opacity
                         * in file plot */
                        ParsePlot(chCard);
                }

                else if( strcmp(chKey4,"POWE") == 0 )
                {
                        /* power law with cutoff, in reads2 */
                        ParsePowerlawContinuum(chCard);
                }

                else if( strcmp(chKey4,"PRIN") == 0 )
                {
                        /* adjust print schedule, in readsun */
                        ParsePrint(chCard);
                }

                else if( strcmp(chKey4,"PUNC") == 0 )
                {
                        /* punch something, in punch */
                        ParsePunch(chCard);
                }

                else if( strcmp(chKey4,"Q(H)") == 0 )
                {
                        /* log of number of ionizing photons */
                        i = 5;
                        if( nqh >= LIMSPC )
                        {
                                /* too many continua were entered */
                                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* silly, but calms down the lint */
                        ASSERT( nqh < LIMSPC );
                        strcpy( rfield.chRSpec[nqh], "4 PI" );
                        strcpy( rfield.chSpNorm[nqh], "Q(H)" );
                        rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( rfield.totpow[nqh] > 100. && called.lgTalk )
                        {
                                fprintf( ioQQQ, " Is this reasonable?\n" );
                        }
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                        /* ParseRangeOption in readsun */
                        ParseRangeOption(nqh,chCard);

                        /* >>chng 06 mar 22, add time option to vary only some continua with time */
                        if( nMatch( "TIME" , chCard ) )
                                rfield.lgTimeVary[nqh] = true;

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* >>chng 03 apr 30, had used log for range option, but if one or other
                                 * was 1 ryd, then became 0 in the log, which was misinterpreted as
                                 * one of the continuum bounds */
                                strcpy( optimize.chVarFmt[optimize.nparm], "Q(H) %f range %f %f" );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh];
                                optimize.vincr[optimize.nparm] = 0.5;
                                /* range option in, but cannot be varied */
                                optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0];
                                optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1];
                                optimize.nvarxt[optimize.nparm] = 3;
                                ++optimize.nparm;
                        }
                        /* increment number of luminosity sources specified */
                        ++nqh;
                }

                else if( strcmp(chKey4,"RATI") == 0 )
                {
                        /* enter a continuum luminosity as a ratio of
                         * nuFnu for this continuum relative to a previous continuum
                         * format; first number is ratio of second to first continuum
                         * second number is energy for this ratio
                         * if third number on line, then 2nd number is energy of
                         * first continuum, while 3rd number is energy of second continuum
                         * in reads2 */
                        ParseRatio(chCard,&nqh);
                }

                else if( strcmp(chKey4,"RADI") == 0 )
                {
                        /* log of inner and outer radii, default second=infinity,
                         * if R2<R1 then R2=R1+R2
                         * there is an optional keyword, "PARSEC" on the line
                         * to use PC as units, reads2 */
                        ParseRadius(chCard,&ar1);
                }

                else if( strcmp(chKey4,"ROBE") == 0 )
                {
                        /* this is the Roberto Terlivich command, no telling if it still works */
                        radius.dRadSign = -1.;
                }

                else if( strcmp(chKey4,"SET ") == 0 )
                {
                        /* set something, in reads2 */
                        ParseSet(chCard);
                }

                else if( strcmp(chKey4,"SPEC") == 0 )
                {
                        /* special key, can do anything */
                        cdEXIT(EXIT_FAILURE);
                }

                else if( strcmp(chKey4,"SPHE") == 0 )
                {
                        /* compute a spherical model, diffuse field from other side in
                         * in reads2 */
                        ParseSphere(chCard);
                }

                else if( strcmp(chKey4,"STAT") == 0 )
                {
                        /* either get or put the code's state as a file on disk */
                        ParseState(chCard);
                }

                else if( strcmp(chKey4,"STOP") == 0 )
                {
                        /* stop model at desired zone, temperature, column density or tau-912
                         * in readsun */
                        ParseStop(chCard);
                }

                else if( strcmp(chKey4,"TABL") == 0 )
                {
                        /* interpolate on input tables for continuum, set of power  laws used
                         * input stored in big BLOCK data
                         * first check that this is the one and only INTERP command
                         * in readsun */
                        ParseTable(&nqh,chCard,&ar1);
                }

                else if( strcmp(chKey4,"TAUM") == 0 )
                {
                        /* minimum optical depths for lines (normally 1e-20)
                         *  (used in STARTER to set lines up) */
                        i = 5;
                        opac.taumin = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
                        if( lgEOL )
                        {
                                NoNumb(chCard);
                        }
                }

                else if( strcmp(chKey4 , "TEST" ) == 0 )
                {
                        /* parse the test command and its options */
                        ParseTest(chCard,&nqh , &ar1, lgDSet);
                }

                else if( strcmp(chKey4,"TIME") == 0 )
                {
                        /* parse the time dependent command, in dynamics.c */
                        ParseDynaTime(chCard);
                }

                else if( strcmp(chKey4,"TITL") == 0 )
                {
                        /* read in title of model starting in col 5 */
                        strcpy( input.chTitle , &input.chOrgCard[5] );
                }

                else if( strcmp(chKey4,"TLAW") == 0 )
                {
                        /* some temperature vs depth law */
                        ParseTLaw(chCard);
                }

                else if( strcmp(chKey4,"TOLE") == 0 )
                {
                        fprintf(ioQQQ,
                                "Sorry, this command has been replaced with the SET TEMPERATURE TOLERANCE command.\n");
                        cdEXIT(EXIT_FAILURE);
                }

                else if( strcmp(chKey4,"TRAC") == 0 )
                {
                        /* turn on trace output, in reads2 */
                        ParseTrace(chCard);
                }

                else if( strcmp(chKey4,"VLAW") == 0 )
                {

                        /* velocity power law as a function of radius. */
                        i = 5;
                        DoppVel.TurbVelLaw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                        DoppVel.lgTurbLawOn = true;
                        /* for now, insist upon non-positive power,
                         * so that velocity decreases with increasing radius. */
                        ASSERT( DoppVel.TurbVelLaw <= 0.f );
                }

                else if( strcmp(chKey4,"TURB") == 0 )
                {
                        /* line has turbulent velocity in km/s */
                        i = 5;
                        DoppVel.TurbVel = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                        /* include turbulent pressure in equation of state?
                         * >>chng 06 mar 24, include turbulent pressure in gas equation of state by default,
                         * but not if NO PRESSURE occurs */
                        if( nMatch(" NO ",chCard) && nMatch("PRES",chCard) )
                        {
                                DoppVel.lgTurb_pressure = false;
                        }
                        else
                        {
                                /* default is to include turbulent pressure in equation of state */
                                DoppVel.lgTurb_pressure = true;
                        }

                        /* true if no number on line - check for equipartition */
                        DoppVel.lgTurbEquiMag = false;

                        if( nMatch("EQUI",chCard) && nMatch("PART",chCard) )
                        {
                                /* turbulence equipartition option */
                                DoppVel.lgTurbEquiMag = true;
                        }

                        if( nMatch(" LOG",chCard) )
                        {
                                DoppVel.TurbVel = (realnum)pow((realnum)10.f,DoppVel.TurbVel);
                        }
                        /* now convert from km/s to cm/s */
                        DoppVel.TurbVel *= 1e5;

                        /* this is optional F parameter from equation 34 of 
                         *>>refer       pressure        turb    Heiles, C. & Troland, T.H. 2005, ApJ, 624, 773 
                         * turbulent energy density will be rho F v^2 / 2 */
                        DoppVel.Heiles_Troland_F = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                /* this is the default value of 3 for isotropic turbulence */
                                DoppVel.Heiles_Troland_F = 3.f;
                        }

                        /* option to have turbulence be dissipative, keyword is dissipate,
                         * argument is log of scale length in cm */
                        if( nMatch("DISS",chCard) )
                        {
                                DoppVel.DispScale = (realnum)pow(10., FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) );
                                if( lgEOL )
                                {
                                        NoNumb(chCard);
                                }
                        }

                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                /* only one parameter to vary */
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "TURBULENCE= %f " );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /* turbulent velocity */
                                optimize.vparm[0][optimize.nparm] = DoppVel.TurbVel/1e5f;
                                optimize.vincr[optimize.nparm] = 0.2f*optimize.vparm[0][optimize.nparm];
                                ++optimize.nparm;
                        }

                        /* set initial turbulence */
                        DoppVel.TurbVelZero = DoppVel.TurbVel;
                }

                else if( strcmp(chKey4,"WIND") == 0 )
                {
                        /* NB - advection and wind commands are now a single command */
                        /* parse the wind command, in dynamics.c */
                        ParseDynaWind(chCard);
                }

                else if( strcmp(chKey2,"XI") == 0 )
                {
                        /* inter x-ray ionization parameter xi=L/r^2 n;
                         * defined per hydrogen
                         * radius must also be entered if spherical, but not needed if plane */
                        ParseIonPar(&nqh,chCard,&ar1);
                }

                /* a line starting with these characters is a comment and so ok 
                 * >>chng 06 sep 04 use routine to check for comments */
                else if( !lgInputComment(chCard) /*(chK1 != '#' && chK1 != '*') && chK1 != '%'*/ )
                {
                        /* end of keys - if we get here key is unrecognized---------- */
                        fprintf( ioQQQ, "  Unrecognized command. Key=\"%4.4s\".  This is routine ParseCommands.\n", 
                          chKey4 );
                        fprintf( ioQQQ, " The line image was==%s==\n", 
                          chCard );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* get next line image, this is tail of while loop that will 
                 * look for lgEOF or blank card */
                input_readarray(chCard,&lgEOF);
        }

        /*------------------------------------------------------------------- */
        /* fall through - hit lgEOF or blank line */

        /* >>chng 00 mar 31, removed dummy call to ParseQuantHeat, PvH */

        for( i=0; i < INPUT_LINE_LENGTH; i++ )
        {
                chCard[i] = ' ';
        }
        chCard[INPUT_LINE_LENGTH-1] = '\0';

        if( called.lgTalk )
        {
                fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
                fprintf( ioQQQ, "%23c***********************************************************************************\n\n\n\n", ' ' );
        }

        /* this is an option to turn on trace printout on the nth
         * call from the optimizer - only allow trace if
         * this is the case and nOptimiz 1 below nTrOpt */
        if( optimize.lgTrOpt )
        {
                /* nTrOpt was set with "optimize trace" command,
                 * is iteration to turn on trace */
                if( optimize.nTrOpt != optimize.nOptimiz )
                {
                        trace.lgTrace = false;
                        /* following overrides turning on trace elsewhere */
                        trace.lgTrOvrd = false;
                }
                else
                {
                        trace.lgTrace = true;
                        called.lgTalk = true;
                        trace.lgTrOvrd = true;
                        fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
                }
        }

        /* set density from DLAW command, must be done here since it may depend on later input commands */
        if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
        {
                dense.gas_phase[ipHYDROGEN] = (realnum)dense_fabden(radius.Radius,radius.depth);

                if( dense.gas_phase[ipHYDROGEN] <= 0. )
                {
                        fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
        {
                dense.gas_phase[ipHYDROGEN] = (realnum)dense_tabden(radius.Radius,radius.depth);

                if( dense.gas_phase[ipHYDROGEN] <= 0. )
                {
                        fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* start checks on parameters set properly - this begins with same line saying start .. */

        /* lgStop_not_enough_info says that not enough info for model, so stop 
         * set true in following tests if anything missing */
        lgStop_not_enough_info = false;
        lgStop = false;

        /* check whether hydrogen density has been set - this value was set to 0 in zero */
        if( dense.gas_phase[ipHYDROGEN] <= 0. )
        {
                fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density MUST be specified.\n" );
                lgStop_not_enough_info = true;
                lgStop = true;
                /* need to set something since used below - will abort
                 * since lgStop is set */
                dense.gas_phase[ipHYDROGEN] = 1.;
        }

        radius.rinner = radius.Radius;

        /* mass flux for wind model - used for mass conservation */
        wind.emdot = dense.gas_phase[ipHYDROGEN]*wind.windv0;

        /* set converge criteria - limit number of iterations and zones */
        if( iterations.lgConverge_set)
        {
                iterations.itermx = MIN2( iterations.itermx , iterations.lim_iter );
                for( j=0; j < iterations.iter_malloc; j++ )
                {
                        geometry.nend[j] = MIN2( geometry.nend[j] , iterations.lim_zone );
                }
        }

        /* NSAVE is number of lines saved, =0 if no commands entered */
        if( input.nSave < 0 )
        {
                fprintf( ioQQQ, " PROBLEM DISASTER No commands were entered - whats up?\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* iterate to convergence and wind models are mutually exclusive */
        if( wind.windv0 > 0. && conv.lgAutoIt )
        {
                if( called.lgTalk )
                {
                        fprintf( ioQQQ, " NOTE PROBLEM Due to the nature of the Sobolev approximation, it makes no sense to converge a windy model.\n" );
                        fprintf( ioQQQ, " NOTE Iterate to convergence is turned off\n\n\n" );
                }
                conv.lgAutoIt = false;
                iterations.itermx = 0;
        }

        /* iterate to convergence and case b do not make sense together */
        /* WJH 22 May 2004: unless we are doing i-front dynamics (negative wind.windv0) */
        if( opac.lgCaseB && conv.lgAutoIt && wind.windv0 >= 0. )
        {
                if( called.lgTalk )
                {
                        fprintf( ioQQQ, " NOTE Case B is an artificial test, it makes no sense to converge this model.\n" );
                        fprintf( ioQQQ, " NOTE Iterate to convergence is turned off.\n\n\n" );
                }
                conv.lgAutoIt = false;
                iterations.itermx = 0;
        }

        /* specifying a density power and constant pressure makes no sense */
        if( dense.DensityPower!=0. && strcmp( dense.chDenseLaw, "CPRE" )==0 )
        {
                if( called.lgTalk )
                {
                        fprintf( ioQQQ, " NOTE Specifying both a density power law and constant pressure is impossible.\n" );
                }
                lgStop = true;
        }

        if( !rfield.lgIonizReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 )
        {
                if( called.lgTalk )
                {
                        fprintf( ioQQQ, " NOTE NO REEVALUATE IONIZATION can only be used with constant density.\n" );
                        fprintf( ioQQQ, " NOTE Resetting to reevaluate ionization.\n\n" );
                }
                rfield.lgIonizReevaluate = true;
        }

        /* check if the combination of stopping column density and H density are physically plausible */
        thickness = min( MIN3( StopCalc.tauend, StopCalc.colpls, StopCalc.colnut ),
                         MIN3( StopCalc.col_h2, StopCalc.col_h2_nut, StopCalc.HColStop ) );
        if( thickness < COLUMN_INIT )
        {
                /* a stop column density was specified - check on physical thickness this corresponds to */
                thickness /= (dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
                /* >> chng 06 dec 21, don't complain if outer radius set small with `stop thickness' command. */
                if( thickness > 1e25 && radius.router[0] > 1e25 )
                {
                        fprintf( ioQQQ, 
                                "NOTE The specified column density and hydrogen density correspond to a thickness of %.2e cm.\n",
                                thickness);
                        fprintf( ioQQQ, 
                                "NOTE This seems large to me.\n");
                        fprintf(ioQQQ,"NOTE a very large radius may cause overflow.\n\n");
                }
        }

        if( gv.lgDColOn && thermal.ConstGrainTemp>0 && called.lgTalk )
        {
                /* warn if constant grain temperature but gas-grain thermal effects
                 * are still included */
                fprintf( ioQQQ, 
                        "NOTE The grain temperatures are set to a constant value with the "
                        "CONSTANT GRAIN TEMPERATURE command, but "
                        "energy exchange \n");
                fprintf( ioQQQ, 
                        "NOTE is still included.  The grain-gas heating-cooling will be incorrect.  "
                        "Consider turning off gas-grain collisional energy\n");
                fprintf( ioQQQ, 
                        "NOTE exchange with the NO GRAIN GAS COLLISIONAL ENERGY EXCHANGE command.\n\n\n");
        }

        /* >>chng 04 nov 27, test on these two */
        if( !rfield.lgDoLineTrans && rfield.lgOpacityFine )
        {
                if( called.lgTalk )
                {
                        fprintf( ioQQQ, " NOTE NO LINE TRANSER set but fine opacities still computed.\n" );
                        fprintf( ioQQQ, " NOTE Turning off fine opacities.\n\n" );
                }
                rfield.lgOpacityFine = false;
        }

        /* >>chng 04 nov 27, test of these two */
        if( h2.lgH2ON && (!rfield.lgDoLineTrans || !rfield.lgOpacityFine) )
        {
                if( called.lgTalk )
                {
                        fprintf( ioQQQ, " NOTE Large H2 molecule turned on but line transfer and fine opacities are not.\n" );
                        fprintf( ioQQQ, " NOTE Turning on line transfer and fine opacities.\n\n" );
                }
                rfield.lgOpacityFine = true;
                rfield.lgDoLineTrans = true;
        }

        if( rfield.lgMustBlockHIon && !rfield.lgBlockHIon )
        {
                /* one of the input continua needs to have H-ionizing radiation
                 * blocked with extinguish command, but it was not done */
                fprintf( ioQQQ, "\n NOTE\n"
                        " NOTE One of the incident continuum is a form used when no H-ionizing radiation is produced.\n" );
                fprintf( ioQQQ, " NOTE You must also include the EXTINGUISH command to make sure this is done.\n" );
                fprintf( ioQQQ, " NOTE The EXTINGUISH command was not included.\n" );
                fprintf( ioQQQ, " NOTE YOU MAY BE MAKING A BIG MISTAKE!!\n NOTE\n\n\n\n" );
        }

        /* if stop temp set below default then we are going into cold and possibly molecular
         * gas - check some parameters in this case */
        if( called.lgTalk && (StopCalc.tend < phycon.TEMP_STOP_DEFAULT || 
                /* thermal.ConstTemp def is zero, set pos when constant temperature is set */
                (thermal.ConstTemp > 0. && thermal.ConstTemp < phycon.TEMP_STOP_DEFAULT ) ) )
        {
                /* print warning if temperature set below default but cosmic rays not turned on */
                /* >>chng 06 mar 02, do not print if molecules are off */
                if( (hextra.cryden == 0.) && !co.lgNoCOMole  && !hmi.lgNoH2Mole )
                {
                        fprintf( ioQQQ, "\n NOTE\n"
                                                        " NOTE The simulation is going into neutral gas but cosmic rays are not included.\n" );
                        fprintf( ioQQQ, " NOTE Ion-molecule chemistry will not occur without a source of ionization.\n" );
                        fprintf( ioQQQ, " NOTE The chemistry network may collapse deep in molecular regions.\n" );
                        fprintf( ioQQQ, " NOTE Consider adding galactic background cosmic rays with the COSMIC RAYS BACKGROUND command.\n" );
                        fprintf( ioQQQ, " NOTE You may be making a BIG mistake.\n NOTE\n\n\n\n" );
                }
        }

        /* dense.gas_phase[ipHYDROGEN] is linear hydrogen density (cm-3) */
        /* test for hydrogen density properly set has already been done above */
        if( called.lgTalk && dense.gas_phase[ipHYDROGEN] < 1e-4 )
        {
                fprintf( ioQQQ, " NOTE Is the entered value of the hydrogen density (%.2e) reasonable?\n",
                        dense.gas_phase[ipHYDROGEN]);
                fprintf( ioQQQ, " NOTE It seems pretty low to me.\n\n\n" );
        }
        else if( called.lgTalk && dense.gas_phase[ipHYDROGEN] > 1e15 )
        {
                fprintf( ioQQQ, " NOTE Is this value of the hydrogen density reasonable?\n" );
                fprintf( ioQQQ, " NOTE It seems pretty high to me.\n\n\n" );
        }

        /* is the model going to crash because of extreme density? */
        if( called.lgTalk && !lgStop && !lgStop_not_enough_info )
        {
                if( dense.gas_phase[ipHYDROGEN] < 1e-6 || dense.gas_phase[ipHYDROGEN] > 1e19 )
                {
                        fprintf( ioQQQ, " NOTE Simulation may crash because of extreme "
                                "density.  The value was %.2e\n\n" , 
                                dense.gas_phase[ipHYDROGEN] );
                }
        }

        if( rfield.nspec == 0 )
        {
                fprintf( ioQQQ, " PROBLEM DISASTER No incident radiation field was specified - "
                        "at least put in the CMB.\n" );
                lgStop = true;
                lgStop_not_enough_info = true;
        }

        if( nqh == 0 )
        {
                fprintf( ioQQQ, " PROBLEM DISASTER Luminosity of continuum MUST be specified.\n" );
                lgStop = true;
                lgStop_not_enough_info = true;
        }

        /* Testing on 0 is safe since this it is initialized that way
         * only print comment if continuum has been specified,
         * if continuum not given then we are aborting anyway */
        if( radius.Radius == 0. && rfield.nspec> 0)
        {
                fprintf( ioQQQ, " PROBLEM DISASTER Starting radius MUST be specified.\n" );
                lgStop = true;
                lgStop_not_enough_info = true;
        }

        if( rfield.nspec != nqh )
        {
                fprintf( ioQQQ, " PROBLEM DISASTER There were not the same number of continuum shapes and luminosities entered.\n" );
                lgStop = true;
        }

        /* we only want to do this test on the first call to the command
         * parser - it will be called many more times but with no grid command
         * during the actual grid calculation */
        static bool lgFirstPass = true;

        /* the NO VARY option sets this flag, and can be used to turn off
         * the grid command, as well as the optimizer */
        if( optimize.lgNoVary && grid.lgGrid )
        {
                /* ignore grids */
                grid.lgGrid = false;
                optimize.nparm = 0;
                grid.nGridCommands = 0;
        }

        if( lgFirstPass && grid.lgGrid && (optimize.nparm!=grid.nGridCommands) )
        {
                /* number of grid vary options do match */
                fprintf( ioQQQ, " PROBLEM DISASTER The GRID command was entered "
                        "but there were %li GRID commands and %li commands with a VARY option.\n" ,
                        grid.nGridCommands , optimize.nparm);
                fprintf( ioQQQ, " There must be the same number of GRIDs and VARY.\n" );
                lgStop = true;
        }
        lgFirstPass = false;

        if( lgStop_not_enough_info )
        {
                fprintf( ioQQQ, " PROBLEM DISASTER I do not have enough information to do the simulation, I cannot go on.\n" );
                fprintf( ioQQQ, "\n\n Sorry.\n\n\n" );
                cdEXIT(EXIT_FAILURE);
        }

        if( lgStop )
        {
                cdEXIT(EXIT_FAILURE);
        }

        /* end checks on parameters set properly - this begins with same line saying start .. */
        return;
}

---