root/branches/newmole/source/cddrive.cpp

/* 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 */
/*cdDrive main routine to call cloudy under all circumstances) */
/*cdReasonGeo write why the model stopped and type of geometry on io file */
/*cdWarnings write all warnings entered into comment stack */
/*cdEms obtain the local emissivity for a line, for the last computed zone */
/*cdColm get the column density for a constituent  */
/*cdLine get the predicted line intensity, also index for line in stack */
/*cdLine_ip get the predicted line intensity, using index for line in stack */
/*cdDLine get the predicted emergent line intensity, also index for line in stack */
/*cdCautions print out all cautions after calculation, on arbitrary io unit */
/*cdTemp_last routine to query results and return temperature of last zone */
/*cdDepth_depth get depth structure from previous iteration */
/*cdTimescales returns thermal, recombination, and H2 formation timescales */
/*cdSurprises print out all surprises on arbitrary unit number */
/*cdNotes print stack of notes about current calculation */
/*cdPressure_last routine to query results and return pressure of last zone */
/*cdTalk tells the code whether to print results or be silent */
/*cdOutp redirect output to arbitrary Fortran unit number */
/*cdRead routine to read in command lines when cloudy used as subroutine */
/*cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit */
/*cdIonFrac get ionization fractions for a constituent */
/*cdTemp get mean electron temperature for any element */
/*cdCooling_last routine to query results and return cooling of last zone */
/*cdHeating_last routine to query results and return heating of last zone */
/*cdEDEN_last return electron density of last zone */
/*cdNoExec call this routine to tell code not to actually execute */
/*cdDate - puts date of code into string */
/*cdVersion produces string that gives version number of the code */
/*cdExecTime any routine can call this, find the time [s] since cdInit was called */
/*cdPrintCommands( FILE *) prints all input commands into file */
/*cdDrive main routine to call cloudy under all circumstances) */
/*cdNwcns get the number of cautions and warnings, to tell if calculation is ok */
/*debugLine provides a debugging hook into the main line array  */
/*cdEms_ip obtain the local emissivity for a line with known index */
/*cdnZone gets number of zones */
/*cdClosePunchFiles closes all the punch files that have been used */
/*cdLineListPunch create a file with a list of all emission lines punched, 
 *and their index within the emission line stack */
/*cdB21cm - returns B as measured by 21 cm */
/*cdPrtWL print line wavelengths in Angstroms in the standard format */
#include "cddefines.h"
#include "trace.h"
#include "conv.h"
#include "init.h"
#include "lines.h"
#include "pressure.h"
#include "prt.h"
#include "colden.h"
#include "dense.h"
#include "radius.h"
#include "struc.h"
#include "mole.h"
#include "elementnames.h"
#include "mean.h"
#include "phycon.h"
#include "called.h"
#include "parse.h"
#include "input.h"
#include "taulines.h"
#include "version.h"
#include "thermal.h"
#include "optimize.h"
#include "grid.h"
#include "timesc.h"
#include "cloudy.h"
#include "warnings.h"
#include "lines_service.h"
#include "cddrive.h"

/*************************************************************************
 *
 * cdDrive - main routine to call cloudy - returns 0 if all ok, 1 if problems
 *
 ************************************************************************/

int cdDrive(void )
{
        bool lgBAD;

        DEBUG_ENTRY( "cdDrive()" );
        /*********************************
         * main routine to call cloudy   *
         *********************************/

        /* this is set false when code loaded, set true when cdInit called,
         * this is check that cdInit was called first */
        if( !lgcdInitCalled )
        {
                printf(" cdInit was not called first - this must be the first call.\n");
                cdEXIT(EXIT_FAILURE);
        }

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, 
                        "cdDrive: lgOptimr=%1i lgVaryOn=%1i lgNoVary=%1i input.nSave:%li\n",
                        optimize.lgOptimr , optimize.lgVaryOn , optimize.lgNoVary, input.nSave );
        }

        /* should we call cloudy, or the optimization driver? 
         * possible to have VARY on line without OPTIMIZE being set 
         * lgNoVary set with "no optimize" command */
        if( optimize.lgOptimr && optimize.lgVaryOn && !optimize.lgNoVary )
                optimize.lgVaryOn = true;
        else
                optimize.lgVaryOn = false;

        /* one time initialization of core load - returns if already called 
         * called here rather than in cdInit since at this point we know if
         * single sim or grid */
        InitCoreload();
        /* count now many sims have been done in this coreload, above set to
         * zero if first call, does nothing on subsequent calls
         * this increment brings to 1 on first sim */
        ++nSimThisCoreload;

        if( optimize.lgVaryOn )
        {
                /* this branch called if optimizing or grid calculation */
                if( trace.lgTrace )
                        fprintf( ioQQQ, "cdDrive: calling grid_do\n");
                /* option to drive optimizer set if OPTIMIZE was in input stream */
                lgBAD = grid_do();
        }
        else
        {
                if( trace.lgTrace )
                        fprintf( ioQQQ, "cdDrive: calling cloudy\n");
                try {

                        /* optimize did not occur, only compute one model, call cloudy */
                        lgBAD = cloudy();
                }
                catch( bad_mpi &bad_info )
                {
                        fprintf( ioQQQ, "An MPI run failed.  The fail mode is %li\n", bad_info.failMode() );
                        ClosePunchFiles( false );
                        lgAbort = true;
                        lgBAD = true;
                }
        }

        /* reset flag saying that cdInit has not been called */
        lgcdInitCalled = false;

        if( lgAbort || lgBAD )
        {
                if( trace.lgTrace )
                        fprintf( ioQQQ, "cdDrive: returning failure during call. \n");
                /* lgAbort set true if something wrong, so return lgBAD false. */
                return(1);
        }
        else
        {
                /* everything is ok, return 0 */
                return(0);
        }
}

/*************************************************************************
*
* cdPrtWL write emission line wavelength
*
************************************************************************/

/*cdPrtWL print line wavelengths in Angstroms in the standard format - 
 * a wrapper for prt_wl which must be kept parallel with sprt_wl
 * both of those live in pdt.c */
void cdPrtWL( FILE *io , realnum wl )
{

        DEBUG_ENTRY( "cdPrtWL()" );

        prt_wl( io , wl );
        return;
}


/*************************************************************************
 *
 * cdReasonGeo write why the model stopped and type of geometry on io file 
 *
 ************************************************************************/


/*cdReasonGeo write why the model stopped and type of geometry on io file */
void cdReasonGeo(FILE * ioOUT)
{

        DEBUG_ENTRY( "cdReasonGeo()" );

        /*this is the reason the calculation stopped*/
        fprintf( ioOUT, "%s", warnings.chRgcln[0] );
        fprintf( ioOUT , "\n" );
        /* this is the geometry */
        fprintf( ioOUT, "%s", warnings.chRgcln[1] );
        fprintf( ioOUT , "\n" );
        return;
}


/*************************************************************************
 *
 * cdWarnings write all warnings entered into comment stack 
 *
 ************************************************************************/

/*cdWarnings write all warnings entered into comment stack */

void cdWarnings(FILE *ioPNT )
{
        long int i;

        DEBUG_ENTRY( "cdWarnings()" );

        if( warnings.nwarn > 0 )
                {
                        for( i=0; i < warnings.nwarn; i++ )
                        {
                                /* these are all warnings that were entered in comment */
                                fprintf( ioPNT, "%s", warnings.chWarnln[i] );
                                fprintf( ioPNT, "\n" );
                        }
                }

        return;
}


/*************************************************************************
 *
 * cdCautions print out all cautions after calculation, on arbitrary io unit 
 *
 ************************************************************************/

/*cdCautions print out all cautions after calculation, on arbitrary io unit */

void cdCautions(FILE * ioOUT)
{
        long int i;

        DEBUG_ENTRY( "cdCautions()" );

        if( warnings.ncaun > 0 )
        {
                for( i=0; i < warnings.ncaun; i++ )
                {
                        fprintf( ioOUT, "%s", warnings.chCaunln[i] );
                        fprintf( ioOUT, "\n" );
                }
        }
        return;
}

/*************************************************************************
 *
 * cdTimescales returns thermal, recombination, and H2 formation timescales 
 *
 ************************************************************************/

void cdTimescales(
        /* the thermal cooling timescale */
        double *TTherm , 
        /* the hydrogen recombination timescale */
        double *THRecom , 
        /* the H2 formation timescale */
        double *TH2 )
{

        DEBUG_ENTRY( "cdTimescales()" );

        /* these were all evaluated in AgeCheck, which was called by PrtComment */

        /* thermal or cooling timescale */
        *TTherm = timesc.time_therm_long;

        /* the hydrogen recombination timescale */
        *THRecom = timesc.time_Hrecom_long;

        /* longer of the the H2 formation and destruction timescales */
        *TH2 = MAX2( timesc.time_H2_Dest_longest , timesc.time_H2_Form_longest );
        return;
}


/*************************************************************************
 *
 * cdSurprises print out all surprises on arbitrary unit number 
 *
 ************************************************************************/

/*cdSurprises print out all surprises on arbitrary unit number */

void cdSurprises(FILE * ioOUT)
{
        long int i;

        DEBUG_ENTRY( "cdSurprises()" );

        if( warnings.nbang > 0 )
        {
                for( i=0; i < warnings.nbang; i++ )
                {
                        fprintf( ioOUT, "%s", warnings.chBangln[i] );
                        fprintf( ioOUT, "\n" );
                }
        }

        return;
}


/*************************************************************************
 *
 * cdNotes print stack of notes about current calculation 
 *
 ************************************************************************/

/*cdNotes print stack of notes about current calculation */

void cdNotes(FILE * ioOUT)
{
        long int i;

        DEBUG_ENTRY( "cdNotes()" );

        if( warnings.nnote > 0 )
        {
                for( i=0; i < warnings.nnote; i++ )
                {
                        fprintf( ioOUT, "%s", warnings.chNoteln[i] );
                        fprintf( ioOUT, "\n" );
                }
        }
        return;
}

/*************************************************************************
 *
 * cdCooling_last routine to query results and return cooling of last zone 
 *
 ************************************************************************/

/*cdCooling_last routine to query results and return cooling of last zone */
double cdCooling_last(void) /* return cooling for last zone */
{
        return (thermal.ctot);
}

/*************************************************************************
 *
 * cdVersion - puts version number of code into string 
 * incoming string must have sufficient length and will become null
 * terminated string
 *
 ************************************************************************/

void cdVersion(char chString[])
{
        strcpy( chString , t_version::Inst().chVersion );
        return;
}

/*************************************************************************
 *
 * cdDate - puts date of code into string 
 * incoming string must have at least 8 char and will become null
 * terminated string
 *
 ************************************************************************/

/* cdDate - puts date of code into string  */
void cdDate(char chString[])
{
        strcpy( chString , t_version::Inst().chDate );
        return;
}


/*************************************************************************
 *
 * cdHeating_last routine to query results and return heating of last zone
 *
 ************************************************************************/

/*cdHeating_last routine to query results and return heating of last zone */

double cdHeating_last(void) /* return heating for last zone */
{
        return (thermal.htot);
}


/*************************************************************************
 *
 * cdEDEN_last return electron density of last zone
 *
 ************************************************************************/

/*cdEDEN_last return electron density of last zone */

double cdEDEN_last(void) /* return electron density for last zone */
{
        return ( dense.eden );
}

/*************************************************************************
 *
 * cdNoExec call this routine to tell code not to actually execute
 *
 ************************************************************************/

/*cdNoExec call this routine to tell code not to actually execute */
#include "noexec.h"

void cdNoExec(void)
{

        DEBUG_ENTRY( "cdNoExec()" );

        /* option to read in all input quantiles and NOT execute the actual model
         * only check on input parameters - set by calling cdNoExec */
        noexec.lgNoExec = true;

        return;
}


/*************************************************************************
 *
 * cdSetExecTime routine to initialize variables keeping track of time at start of calculation
 *
 ************************************************************************/

/* set to false initially, then to true when cdSetExecTime() is called to
 * initialize the clock */
static bool lgCalled=false;

/* >>chng 06 dec 19, RP rm "|| defined(__HP_aCC)" to run on hp */
#if defined(_MSC_VER) 
/* _MSC_VER branches assume getrusage isn't implemented by MS 
 * also is not implemented on our HP superdome */
struct timeval {
        long    tv_sec;         /* seconds */
        long    tv_usec;        /* microseconds */
};
#else
#include <sys/time.h>
#include <sys/resource.h>
#endif

/* will be used to save initial time */
static struct timeval before;

/* cdClock stores time since arbitrary datum in clock_dat           */
STATIC void cdClock(struct timeval *clock_dat)
{
        DEBUG_ENTRY( "cdClock()" );

/* >>chng 06 sep 2 rjrw: use long recurrence, fine grain UNIX clock *
 * -- maintain system dependences in a single routine               */
#if defined(_MSC_VER) || defined(__HP_aCC)
        clock_t clock_val;
        /* >>chng 05 dec 21, from above to below due to negative exec times */
        /*return(  (double)(clock() - before) / (double)CLOCKS_PER_SEC );*/
        clock_val = clock();
        clock_dat->tv_sec = clock_val/CLOCKS_PER_SEC;
        clock_dat->tv_usec = 1000000*(clock_val-(clock_dat->tv_sec*CLOCKS_PER_SEC))/CLOCKS_PER_SEC;
        /*>>chng 06 oct 05, this produced very large number, time typically 50% too long 
        clock_dat->tv_usec = 0;*/
#else
        struct rusage rusage;
        if(getrusage(RUSAGE_SELF,&rusage) != 0)
        { 
                fprintf( ioQQQ, "DISASTER cdClock called getrusage with invalid arguments.\n" );
                fprintf( ioQQQ, "Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }
        clock_dat->tv_sec = rusage.ru_utime.tv_sec;
        clock_dat->tv_usec = rusage.ru_utime.tv_usec;
#endif

}
/* cdSetExecTime is called by cdInit when everything is initialized,
 * so that every time cdExecTime is called the elapsed time is returned */
void cdSetExecTime(void)
{
        cdClock(&before);
        lgCalled = true;
}
/* cdExecTime returns the elapsed time cpu time (sec) that has elapsed 
 * since cdInit called cdSetExecTime.*/
double cdExecTime(void)
{
        DEBUG_ENTRY( "cdExecTime()" );

        struct timeval clock_dat;
        /* check that we were properly initialized */
        if( lgCalled)
        {
                cdClock(&clock_dat);
                /*fprintf(ioQQQ,"\n DEBUG sec %.2e usec %.2e\n",
                        (double)(clock_dat.tv_sec-before.tv_sec),
                        1e-6*(double)(clock_dat.tv_usec-before.tv_usec));*/
                return (double)(clock_dat.tv_sec-before.tv_sec)+1e-6*(double)(clock_dat.tv_usec-before.tv_usec);
        }
        else
        {
                /* this is a big problem, we were asked for the elapsed time but
                 * the timer was not initialized by calling SetExecTime */
                fprintf( ioQQQ, "DISASTER cdExecTime was called before SetExecTime, impossible.\n" );
                fprintf( ioQQQ, "Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }
}

/*************************************************************************
 *
 * cdPrintCommands prints all input commands into file
 *
 ************************************************************************/

/* cdPrintCommands( FILE *)
 * prints all input commands into file */
void cdPrintCommands( FILE * ioOUT )
{
        long int i;
        fprintf( ioOUT, " Input commands follow:\n" );
        fprintf( ioOUT, "c ======================\n" );

        for( i=0; i <= input.nSave; i++ )
        {
                fprintf( ioOUT, "%s\n", input.chCardSav[i] );
        }
        fprintf( ioOUT, "c ======================\n" );
}


/*************************************************************************
 *
 * cdEms obtain the local emissivity for a line, for the last computed zone
 *
 ************************************************************************/

/* \todo 2 This routine, cdLine, cdEmis_ip, and cdLine_ip should be consolidated somewhat.
 * in particular so that this routine has the same "closest line" reporting that cdLine has. */
long int cdEmis(
        /* return value will be index of line within stack,
         * negative of number of lines in the stack if the line could not be found*/
        /* 4 char null terminated string label */
        char *chLabel,
        /* line wavelength */
        realnum wavelength, 
        /* the vol emissivity of this line in last computed zone */
        double *emiss )
{
        /* use following to store local image of character strings */
        char chCARD[INPUT_LINE_LENGTH];
        char chCaps[5];
        long int j;
        realnum errorwave;

        DEBUG_ENTRY( "cdEms()" );

        /* routine returns the emissivity in the desired line
         * only used internally by the code, to do punch lines structure */

        strcpy( chCARD, chLabel );

        /* make sure chLabel is all caps */
        caps(chCARD);/* convert to caps */

        /* get the error associated with 4 significant figures */
        errorwave = WavlenErrorGet( wavelength );

        for( j=0; j < LineSave.nsum; j++ )
        {
                /* change chLabel to all caps to be like input chLineLabel */
                cap4(chCaps , (char*)LineSv[j].chALab);

                /* check wavelength and chLabel for a match */
                /*if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength)<errorwave && 
                        strcmp(chCaps,chCARD) == 0 ) */
                if( fabs(LineSv[j].wavelength-wavelength) < errorwave && strcmp(chCaps,chCARD) == 0 )
                {
                        /* match, so set emiss to emissivity in line */
                        *emiss = LineSv[j].emslin[LineSave.lgLineEmergent];
                        /* and announce success by returning line index within stack */
                        return j;
                }
        }

        /* we fell through without finding the line - return false */
        return -LineSave.nsum;
}

/*************************************************************************
 *
 * cdEms_ip obtain the local emissivity for a line with known index
 *
 ************************************************************************/

void cdEmis_ip(
        /* index of the line in the stack */
        long int ipLine, 
        /* the vol emissivity of this line in last computed zone */
        double *emiss )
{
        DEBUG_ENTRY( "cdEmis_ip()" );

        /* returns the emissivity in a line - implements punch lines structure 
         * this version uses previously stored line index to save speed */
        ASSERT( ipLine >= 0 && ipLine < LineSave.nsum );
        *emiss = LineSv[ipLine].emslin[LineSave.lgLineEmergent];
        return;
}

/*************************************************************************
 *
 * cdColm get the column density for a constituent - 0 return if ok, 1 if problems 
 *
 ************************************************************************/

int cdColm(
        /* return value is zero if all ok, 1 if errors happened */
        /* 4-char + eol string that is first 
         * 4 char of element name as spelled by Cloudy, upper or lower case */
        const char *chLabel,    

        /* integer stage of ionization, 1 for atom, 2 for A+, etc, 
         * 0 is special flag for CO, H2, OH, or excited state */
        long int ion,

        /* the column density derived by the code [cm-2] */
        double *theocl )        
{
        long int nelem;
        /* use following to store local image of character strings */
        char chLABEL_CAPS[20];

        DEBUG_ENTRY( "cdColm()" );

        /* check that chLabel[4] is null - supposed to be 4 char + end */
        if( chLabel[4]!=0 )
        {
                fprintf( ioQQQ, " cdColm called with insane ion label, =%s, must be 4 character + end of string.\n", 
                  chLabel );
                return(1);
        }

        strcpy( chLABEL_CAPS, chLabel );
        /* convert element label to all caps */
        caps(chLABEL_CAPS);

        /* zero ionization stage has special meaning.  The quantities recognized are
         * the molecules, "H2  ", "OH  ", "CO  ", etc 
         * "CII*" excited state C+ */
        if( ion < 0 )
        {
                fprintf( ioQQQ, " cdColm called with insane ion, =%li\n", 
                  ion );
                return(1);
        }

        else if( ion == 0 )
        {
                /* this case molecular column density */
                /* want the molecular hydrogen H2 column density */
                if( strcmp( chLABEL_CAPS , "H2  " )==0 )
                {
                        *theocl = colden.colden[ipCOL_H2g] + colden.colden[ipCOL_H2s];
                }

                /* H- column density  */
                else if( strcmp( chLABEL_CAPS , "H-  " )==0 )
                {
                        *theocl = colden.colden[ipCOL_HMIN];
                }

                /* H2+ column density ipCOL_H2p is 4 */
                else if( strcmp( chLABEL_CAPS , "H2+ " )==0 )
                {
                        *theocl = colden.colden[ipCOL_H2p];
                }

                /* H3+ column density */
                else if( strcmp( chLABEL_CAPS , "H3+ " )==0 )
                {
                        *theocl = colden.colden[ipCOL_H3p];
                }

                /* H2g - ground H2 column density */
                else if( strcmp( chLABEL_CAPS , "H2G " )==0 )
                {
                        *theocl = colden.colden[ipCOL_H2g];
                }

                /* H2* - excited H2 - column density */
                else if( strcmp( chLABEL_CAPS , "H2* " )==0 )
                {
                        *theocl = colden.colden[ipCOL_H2s];
                }

                /* HeH+ column density */
                else if( strcmp( chLABEL_CAPS , "HEH+" )==0 )
                {
                        *theocl = colden.colden[ipCOL_HeHp];
                }

                /* carbon monoxide column density */
                else if( strcmp( chLABEL_CAPS , "CO  " )==0 )
                {
                        *theocl = findspecies("CO")->column;
                }

                /* OH column density */
                else if( strcmp( chLABEL_CAPS , "OH  " )==0 )
                {
                        *theocl = findspecies("OH")->column;
                }

                /* H2O column density */
                else if( strcmp( chLABEL_CAPS , "H2O " )==0 )
                {
                        *theocl = findspecies("H2O")->column;
                }

                /* O2 column density */
                else if( strcmp( chLABEL_CAPS , "O2  " )==0 )
                {
                        *theocl = findspecies("O2")->column;
                }

                /* SiO column density */
                else if( strcmp( chLABEL_CAPS , "SIO " )==0 )
                {
                        *theocl = findspecies("SiO")->column;
                }

                /* C2 column density */
                else if( strcmp( chLABEL_CAPS , "C2  " )==0 )
                {
                        *theocl = findspecies("C2")->column;
                }

                /* C3 column density */
                else if( strcmp( chLABEL_CAPS , "C3  " )==0 )
                {
                        *theocl = findspecies("C3")->column;
                }

                /* CN column density */
                else if( strcmp( chLABEL_CAPS , "CN  " )==0 )
                {
                        *theocl = findspecies("CN")->column;
                }

                /* CH column density */
                else if( strcmp( chLABEL_CAPS , "CH  " )==0 )
                {
                        *theocl = findspecies("CH")->column;
                }

                /* CH+ column density */
                else if( strcmp( chLABEL_CAPS , "CH+ " )==0 )
                {
                        *theocl = findspecies("CH+")->column;
                }

                /* ===========================================================*/
                /* end special case molecular column densities, start special cases
                 * excited state column densities */
                /* CII^* column density, population of J=3/2 upper level of split ground term */
                else if( strcmp( chLABEL_CAPS , "CII*" )==0 )
                {
                        *theocl = colden.C2Colden[1];
                }
                else if( strcmp( chLABEL_CAPS , "C11*" )==0 )
                {
                        *theocl = colden.C1Colden[0];
                }
                else if( strcmp( chLABEL_CAPS , "C12*" )==0 )
                {
                        *theocl = colden.C1Colden[1];
                }
                else if( strcmp( chLABEL_CAPS , "C13*" )==0 )
                {
                        *theocl = colden.C1Colden[2];
                }
                else if( strcmp( chLABEL_CAPS , "O11*" )==0 )
                {
                        *theocl = colden.O1Colden[0];
                }
                else if( strcmp( chLABEL_CAPS , "O12*" )==0 )
                {
                        *theocl = colden.O1Colden[1];
                }
                else if( strcmp( chLABEL_CAPS , "O13*" )==0 )
                {
                        *theocl = colden.O1Colden[2];
                }
                /* CIII excited states, upper level of 1909 */
                else if( strcmp( chLABEL_CAPS , "C30*" )==0 )
                {
                        *theocl = colden.C3Colden[1];
                }
                else if( strcmp( chLABEL_CAPS , "C31*" )==0 )
                {
                        *theocl = colden.C3Colden[2];
                }
                else if( strcmp( chLABEL_CAPS , "C32*" )==0 )
                {
                        *theocl = colden.C3Colden[3];
                }
                else if( strcmp( chLABEL_CAPS , "SI2*" )==0 )
                {
                        *theocl = colden.Si2Colden[1];
                }
                else if( strcmp( chLABEL_CAPS , "HE1*" )==0 )
                {
                        *theocl = colden.He123S;
                }
                /* special option, "H2vJ" */
                else if( strncmp(chLABEL_CAPS , "H2" , 2 ) == 0 )
                {
                        long int iVib = chLABEL_CAPS[2] - '0';
                        long int iRot = chLABEL_CAPS[3] - '0';
                        if( iVib<0 || iRot < 0 )
                        {
                                fprintf( ioQQQ, " cdColm called with insane v,J for H2=\"%4.4s\" caps=\"%4.4s\"\n", 
                                  chLabel , chLABEL_CAPS );
                                return( 1 );
                        }
                        *theocl = cdH2_colden( iVib , iRot );
                }

                /* clueless as to what was meant - bail */
                else
                {
                        fprintf( ioQQQ, " cdColm called with unknown element chLabel, org=\"%4.4s \" 0 caps=\"%4.4s\" 0\n", 
                          chLabel , chLABEL_CAPS );
                        return(1);
                }
        }
        else
        {
                /* this case, ionization stage of some element */
                /* find which element this is */
                nelem = 0;
                while( nelem < LIMELM && 
                        strncmp(chLABEL_CAPS,elementnames.chElementNameShort[nelem],4) != 0 )
                {
                        ++nelem;
                }

                /* this is true if we have one of the first 30 elements in the label,
                 * nelem is on C scale */
                if( nelem < LIMELM )
                {

                        /* sanity check - does this ionization stage exist?
                         * max2 is to pick up H2 as H 3 */
                        if( ion > MAX2(3,nelem + 2) )
                        {
                                fprintf( ioQQQ, 
                                  " cdColm asked to return ionization stage %ld for element %s but this is not physical.\n", 
                                  ion, chLabel );
                                return(1);
                        }

                        /* the column density, ion is on physics scale, but means are on C scale */
                        *theocl = mean.xIonMeans[0][nelem][ion-1];
                        /*>>chng 06 jan 23, div by factor of two
                         * special case of H2 when being tricked as H 3 - this stores 2H_2 so that
                         * the fraction of H in H0 and H+ is correct - need to remove this extra
                         * factor of two here */
                        if( nelem==ipHYDROGEN && ion==3 )
                                *theocl /= 2.;
                }
                else
                {
                        fprintf( ioQQQ, 
                          " cdColm did not understand this combination of ion %4ld and element %4.4s.\n", 
                          ion, chLabel );
                        return(1);
                }
        }
        return(0);
}


/*************************************************************************
 *
 *cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit 
 *
 ************************************************************************/

void cdErrors(FILE *ioOUT)
{
        long int nc, 
          nn, 
          npe, 
          ns, 
          nte, 
          nw ,
          nIone,
          nEdene;
        bool lgAbort_loc;

        DEBUG_ENTRY( "cdErrors()" );

        /* first check for number of warnings, cautions, etc */
        cdNwcns(&lgAbort_loc,&nw,&nc,&nn,&ns,&nte,&npe, &nIone, &nEdene );

        /* only say something is one of these problems is nonzero */
        if( nw || nc || nte || npe ||   nIone || nEdene || lgAbort_loc )
        {
                /* say the title of the model */
                fprintf( ioOUT, "%75.75s\n", input.chTitle );

                if( lgAbort_loc )
                        fprintf(ioOUT," Calculation ended with abort!\n");

                /* print warnings on the io unit */
                if( nw != 0 )
                {
                        cdWarnings(ioOUT);
                }

                /* print cautions on the io unit */
                if( nc != 0 )
                {
                        cdCautions(ioOUT);
                }

                if( nte != 0 )
                {
                        fprintf( ioOUT , "Te failures=%4ld\n", nte );
                }

                if( npe != 0 )
                {
                        fprintf( ioOUT , "Pressure failures=%4ld\n", npe );
                }

                if( nIone != 0 )
                {
                        fprintf( ioOUT , "Ionization failures=%4ld\n", nte );
                }

                if( nEdene != 0 )
                {
                        fprintf( ioOUT , "Electron density failures=%4ld\n", npe );
                }
        }
        return;
}

/*************************************************************************
 *
 *cdDepth_depth get depth structure from previous iteration 
 *
 ************************************************************************/
void cdDepth_depth( double cdDepth[] )
{
        lo