/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*atmdat_readin read in some data files, but only if this is very first call, 
 * called by Cloudy */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "mewecoef.h"
#include "heavy.h"
#include "mole.h"
#include "yield.h"
#include "path.h"
#include "trace.h"
#include "lines.h"
#include "lines_service.h"
#include "ionbal.h"
#include "struc.h"
#include "geometry.h"
#include "dynamics.h"
#include "elementnames.h"
#include "hyperfine.h"
#include "atmdat.h"
/* */
/* this was needed to get array to crash out of bounds if not set.
 * std limits on limits.h did not work with visual studio! */
#define INTBIG 2000000000

/* these are the individual pointers to the level 1 lines, they are set to
 * very large negative.  
 * NB NB NB!!
 * these occur two times in the code!!
 * They are declared in taulines.h, and defined here  */
long 
   ipT1656=INTBIG,    ipT9830=INTBIG,  /*ipH21cm=INTBIG, ipHe3cm=INTBIG,*/
   ipT8727=INTBIG,    ipT1335=INTBIG,  ipT1909=INTBIG,
	ipT977=INTBIG,    ipT1550=INTBIG,  ipT1548=INTBIG, ipT386=INTBIG, 
	ipT310=INTBIG,    ipc31175=INTBIG, ipT291=INTBIG,  ipT280=INTBIG,
	ipT274=INTBIG,    ipT270=INTBIG,   ipT312=INTBIG,  ipT610=INTBIG, 
	ipT370=INTBIG,    ipT157=INTBIG,   ipT1085=INTBIG, 
	ipT990=INTBIG,    ipT1486=INTBIG,  ipT765=INTBIG,  ipT1243=INTBIG, 
	ipT1239=INTBIG,   ipT374g=INTBIG,  ipT374x=INTBIG, ipT1200=INTBIG,
	ipT2140=INTBIG,   ipT671=INTBIG,   ipT315=INTBIG,  ipT324=INTBIG, 
	ipT333=INTBIG,    ipT209=INTBIG,   ipT122=INTBIG,  ipT205=INTBIG,
	ipT57=INTBIG,     ipT6300=INTBIG,  ipT6363=INTBIG, ipT5577=INTBIG, 
	ipT834=INTBIG,    ipT1661=INTBIG,  ipT1666=INTBIG, ipT835=INTBIG,
	ipT789=INTBIG,    ipT630=INTBIG,   ipT1304=INTBIG, ipSi10_606=INTBIG,
	ipT1039=INTBIG,   ipT8446=INTBIG,  ipT4368=INTBIG, ipTOI13=INTBIG,
	ipTOI11=INTBIG,   ipTOI29=INTBIG,  ipTOI46=INTBIG, ipTO1025=INTBIG, 
	ipT304=INTBIG,    ipT1214=INTBIG,  ipT150=INTBIG,  ipT146=INTBIG,
	ipT63=INTBIG,     ipTO88=INTBIG,   ipT52=INTBIG,   ipT26=INTBIG, 
	ipT1032=INTBIG,   ipT1037=INTBIG,  ipF0229=INTBIG, ipF0267=INTBIG,
	ipF444=INTBIG,    ipF425=INTBIG,   ipT770=INTBIG, ipT780=INTBIG, 
	ipxNe0676=INTBIG, ipT895=INTBIG,   ipT88=INTBIG, ipTNe13=INTBIG,
	ipTNe36=INTBIG,   ipTNe16=INTBIG,  ipTNe14=INTBIG, ipTNe24=INTBIG, 
	ipT5895=INTBIG,   ipfsNa373=INTBIG,ipfsNa490=INTBIG, ipfsNa421=INTBIG,
	ipxNa6143=INTBIG, ipxNa6862=INTBIG,ipxNa0746=INTBIG, ipMgI2853=INTBIG, 
	ipMgI2026=INTBIG, ipT2796=INTBIG,  ipT2804=INTBIG,
	ipT705=INTBIG,    ipT4561=INTBIG,  ipxMg51325=INTBIG, ipxMg52417=INTBIG, 
	ipxMg52855=INTBIG,ipxMg71190=INTBIG, ipxMg72261=INTBIG,
	ipxMg72569=INTBIG,ipxMg08303=INTBIG, ipTMg610=INTBIG, ipTMg625=INTBIG, 
	ipT58=INTBIG,     ipTMg4=INTBIG, ipTMg14=INTBIG, ipTMg6=INTBIG,
	ipfsMg790=INTBIG, ipfsMg755=INTBIG, ipAlI3957=INTBIG, ipAlI3090=INTBIG, 
	ipT1855=INTBIG,   ipT1863=INTBIG, ipT2670=INTBIG,
	ipAl529=INTBIG,   ipAl6366=INTBIG, ipAl6912=INTBIG, ipAl8575=INTBIG, 
	ipAl8370=INTBIG,  ipAl09204=INTBIG, ipT639=INTBIG,
	ipTAl550=INTBIG,  ipTAl568=INTBIG, ipTAl48=INTBIG, ipSii2518=INTBIG, 
	ipSii2215=INTBIG, ipT1808=INTBIG,
	ipT1207=INTBIG,   ipT1895=INTBIG, ipT1394=INTBIG, ipT1403=INTBIG, 
	ipT1527=INTBIG,   ipT1305=INTBIG, ipT1260=INTBIG, ipSi619=INTBIG,
	ipSi10143=INTBIG, ipTSi499=INTBIG, ipTSi521=INTBIG, ipTSi41=INTBIG, 
	ipTSi35=INTBIG,   ipTSi25=INTBIG, ipTSi65=INTBIG,
	ipTSi3=INTBIG,    ipTSi4=INTBIG, ipP0260=INTBIG, ipP0233=INTBIG, 
	ipP0318=INTBIG,   ipP713=INTBIG, ipP848=INTBIG, ipP817=INTBIG,
	ipP1027=INTBIG,   ipP1018=INTBIG, ipT1256=INTBIG, ipT1194=INTBIG, 
	ipTS1720=INTBIG,  ipT1198=INTBIG, ipT786=INTBIG,
	ipT933=INTBIG,    ipT944=INTBIG, ipfsS810=INTBIG, ipfsS912=INTBIG, 
	ipfsS938=INTBIG,  ipfsS1119=INTBIG, ipfsS1114=INTBIG, ipfsS1207=INTBIG,
	ipTSu418=INTBIG,  ipTSu446=INTBIG, ipTSu30=INTBIG, ipTS19=INTBIG, 
	ipTS34=INTBIG,    ipTS11=INTBIG, ipfsCl214=INTBIG, ipfsCl233=INTBIG,
	ipCl04203=INTBIG, ipCl04117=INTBIG, ipCl973=INTBIG, ipCl1030=INTBIG, 
	ipCl1092=INTBIG,  ipT354=INTBIG, ipT389=INTBIG, ipT25=INTBIG,
	ipTAr7=INTBIG,    ipTAr9=INTBIG, ipTAr22=INTBIG, ipTAr13=INTBIG, 
	ipTAr8=INTBIG,    ipAr06453=INTBIG, ipAr1055=INTBIG, ipAr1126=INTBIG,
	ipAr1178=INTBIG,  ipKI7745=INTBIG, ipxK03462=INTBIG, ipxK04598=INTBIG, 
	ipxK04154=INTBIG, ipxK06882=INTBIG, ipxK06557=INTBIG,
	ipxK07319=INTBIG, ipxK11425=INTBIG, ipCaI4228=INTBIG, ipT3934=INTBIG, 
	ipT3969=INTBIG,   ipT8498=INTBIG, ipT8542=INTBIG,
	ipT8662=INTBIG,   ipT7291=INTBIG, ipT7324=INTBIG, ipTCa302=INTBIG, 
	ipTCa345=INTBIG,  ipTCa19=INTBIG, ipTCa3=INTBIG, ipTCa12=INTBIG,
	ipTCa4=INTBIG,    ipCa0741=INTBIG, ipCa0761=INTBIG, ipCa08232=INTBIG, 
	ipCa12333=INTBIG, ipSc05231=INTBIG, ipSc13264=INTBIG,
	ipTi06172=INTBIG, ipTi14212=INTBIG, ipVa07130=INTBIG, ipVa15172=INTBIG, 
	ipCr08101=INTBIG, ipCr16141=INTBIG, ipxMn0979=INTBIG,
	ipxMn1712=INTBIG, ipFeI3884=INTBIG, ipFeI3729=INTBIG, ipFeI3457=INTBIG, 
	ipFeI3021=INTBIG, ipFeI2966=INTBIG, ipTuv3=INTBIG,
	ipTr48=INTBIG,    ipTFe16=INTBIG, ipTFe26=INTBIG, ipTFe34=INTBIG, 
	ipTFe35=INTBIG,   ipTFe46=INTBIG, ipTFe56=INTBIG, ipT1122=INTBIG,
	ipFe0795=INTBIG,  ipFe0778=INTBIG, ipT245=INTBIG, ipT352=INTBIG, 
	ipFe106375=INTBIG,ipT353=INTBIG, /*ipFe1310=INTBIG, ipFe1311=INTBIG,*/
	ipT347=INTBIG,    ipT192=INTBIG, ipT255=INTBIG, ipT11=INTBIG, 
	ipT191=INTBIG,    /*ipTFe07=INTBIG, ipTFe61=INTBIG,*/ ipFe18975=INTBIG, ipTFe23=INTBIG,
	ipTFe13=INTBIG,   ipCo11527=INTBIG, ipxNi1242=INTBIG ;
long 	ipS4_1405=INTBIG,ipS4_1398=INTBIG,ipS4_1424=INTBIG,ipS4_1417=INTBIG,ipS4_1407=INTBIG,
		ipO4_1400=INTBIG,ipO4_1397=INTBIG,ipO4_1407=INTBIG,ipO4_1405=INTBIG,ipO4_1401=INTBIG,
		ipN3_1749=INTBIG,ipN3_1747=INTBIG,ipN3_1754=INTBIG,ipN3_1752=INTBIG,ipN3_1751=INTBIG,
		ipC2_2325=INTBIG,ipC2_2324=INTBIG,ipC2_2329=INTBIG,ipC2_2328=INTBIG,ipC2_2327=INTBIG,
		ipSi2_2334=INTBIG,ipSi2_2329=INTBIG,ipSi2_2350=INTBIG,ipSi2_2344=INTBIG,ipSi2_2336=INTBIG,
		ipFe22_247=INTBIG,ipFe22_217=INTBIG,ipFe22_348=INTBIG,ipFe22_292=INTBIG,
		ipFe22_253=INTBIG,ipFe22_846=INTBIG,ipTFe20_721=INTBIG,ipTFe20_578=INTBIG , ipZn04363=INTBIG,
		ipS12_520=INTBIG,ipS1_25m=INTBIG,ipS1_56m=INTBIG,ipCl1_11m=INTBIG,
		ipFe1_24m=INTBIG,ipFe1_35m=INTBIG,ipFe1_54m=INTBIG,ipFe1_111m=INTBIG,
		ipNi1_7m=INTBIG ,ipNi1_11m=INTBIG,ipSi1_130m=INTBIG,ipSi1_68m=INTBIG;

/* above are the individual pointers to the level 1 lines, they are set to
 * very large negative.  
 * NB NB NB!!
 * these occur two times in the code!!
 * They are declared in taulines.h, and defined here  */

/* definition for whether level 2 lines are enabled, will be set to -1 
 * with no level2 command */
/*long nWindLine = NWINDDIM;*/
/*float TauLine2[NWINDDIM][NTA];*/
/*float **TauLine2;*/
#include "atomfeii.h"

void atmdat_readin(void)
{
	long int i, 
	  iCase ,
	  ion, 
	  ipDens ,
	  ipISO ,
	  ipTemp ,
	  j, 
	  J,
	  ig0, 
	  ig1, 
	  imax, 
	  nDataSet,
	  nelem, 
	  nelec, 
	  ns ,
	  magic1, 
	  magic2,
	  mol ,
	  nUTA_Romas,
#	define	NIONROMAS	26
	  nFeIonRomas[NIONROMAS];

	char cha;
	char chS2[3];

	long ipZ;
	int lgEOL;

	float temp[15];

	FILE *ioDATA, *ioROMAS , *ioBEHAR ;
	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;

	static int lgFirstCall = TRUE;

#	ifdef DEBUG_FUN
	fputs( "<+>atmdat_readin()\n", debug_fp );
#	endif

	/* do nothing if not first call */
	if( !lgFirstCall )
	{ 
		/* do not do anything, but make sure that number of zones has not increased */
		int lgTooBig = FALSE;
		for( j=0; j < ITRDIM; j++ )
		{
			if( geometry.nend[j]>=struc.nzlim )
				lgTooBig = TRUE;
		}
		if( lgTooBig )
		{
			fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
			fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
			fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
			fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
			fprintf(ioQQQ," Sorry.\n");
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}

#		ifdef DEBUG_FUN
		fputs( " <->atmdat_readin()\n", debug_fp );
#		endif
		return;
	}

	lgFirstCall = FALSE; /* do not reevaluate again */

	/*Badnell_rec_init This code is written by Terry Yun, 2005 *
	 * It reads dielectronic recombination rate coefficient fits into 3D arrays */
	Badnell_rec_init();

	/* create space for the structure variables */
	/* nzlim will be limit, and number allocated */
	/* >>chng 01 jul 28, define this var, do all following mallocs */
	for( j=0; j < ITRDIM; j++ )
	{
		struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] );
	}

	/* sloppy, but add one extra for safely */
	++struc.nzlim ;

	if( (struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )))) == NULL )
		BadMalloc();

	if( (struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )))) == NULL )
		BadMalloc();

	if( (struc.testr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.volstr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.drad_x_fillfac = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.histr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.hiistr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.ednstr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.o3str = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.pressure = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.PresRadCurr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.GasPressure = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.hden = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.DenParticles = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.DenMass = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.drad = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.depth = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.xLyman_depth = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (struc.Hmolec = ((float**)MALLOC( (size_t)(N_H_MOLEC)*sizeof(float* )))) == NULL )
		BadMalloc();

	if( (struc.COmolec = ((float**)MALLOC( (size_t)(NUM_COMOLE_CALC)*sizeof(float* )))) == NULL )
		BadMalloc();

	for (mol=0;mol<N_H_MOLEC;mol++)
	{
		if( (struc.Hmolec[mol] = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
			BadMalloc();
	}

	for (mol=0;mol<NUM_COMOLE_CALC;mol++)
	{
		if( (struc.COmolec[mol] = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
			BadMalloc();
	}

	/* create space for gas phase abundances array, first create space for the elements */
	if( (struc.gas_phase = (float **)MALLOC(sizeof(float *)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();
	/* last the zones */
	for( ipZ=0; ipZ< LIMELM ;++ipZ )
	{
		if( (struc.gas_phase[ipZ] = 
			(float*)MALLOC(sizeof(float)*(unsigned)(struc.nzlim) ))==NULL )
			BadMalloc();
	}

	/* create space for struc.xIonDense array, first create space for the zones */
	if( (struc.xIonDense = (float ***)MALLOC(sizeof(float **)*(unsigned)(LIMELM+3) ) )==NULL )
		BadMalloc();

	for( ipZ=0; ipZ< (LIMELM+3) ;++ipZ )
	{
		/* space for the ionizaiton stage */
		if( (struc.xIonDense[ipZ] = 
			(float**)MALLOC(sizeof(float*)*(unsigned)(LIMELM+1) ))==NULL )
			BadMalloc();
		/* now create diagonal of space for zones */
		for( ion=0; ion < (LIMELM+1); ++ion )
		{
			if( (struc.xIonDense[ipZ][ion] = 
				(float*)MALLOC(sizeof(float)*(unsigned)(struc.nzlim) ))==NULL )
				BadMalloc();
		}
	}

	/* some structure variables */
	for( i=0; i < struc.nzlim; i++ )
	{
		struc.testr[i] = 0.;
		struc.volstr[i] = 0.;
		struc.drad_x_fillfac[i] = 0.;
		struc.histr[i] = 0.;
		struc.hiistr[i] = 0.;
		struc.ednstr[i] = 0.;
		struc.o3str[i] = 0.;
		struc.heatstr[i] = 0.;
		struc.coolstr[i] = 0.;
		struc.pressure[i] = 0.;
		struc.PresRadCurr[i] = 0.;
		struc.GasPressure[i] = 0.;
		struc.DenParticles[i] = 0.;
		struc.depth[i] = 0.;
		for (mol=0;mol<N_H_MOLEC;mol++)
		{
			struc.Hmolec[mol][i] = 0.;
		}
		for (mol=0;mol<NUM_COMOLE_CALC;mol++)
		{
			struc.COmolec[mol][i] = 0.;
		}
	}

	/* allocate space for some arrays used by dynamics routines, and zero out vars */
	DynaCreateArrays( );

	/*************************************************************
	 *                                                           *
	 * get the level 1 lines, with real collision data set       *
	 *                                                           *
	 *************************************************************/

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "level1.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "level1.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," atmdat_readin opening level1.dat:");

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " atmdat_readin could not open level1.dat\n" );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " atmdat_readin could not open level1.dat\n");
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* first line is a version number and does not count */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* count how many lines are in the file, ignoring all lines
	 * starting with '#' */
	nLevel1 = 0;
	while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		/* we want to count the lines that do not start with #
		 * since these contain data */
		if( chLine[0] != '#')
			++nLevel1;
	}

	/* now MALLOC the TauLines array */
	if( (TauLines = ((EmLine *)MALLOC( (size_t)(nLevel1+1)*sizeof(EmLine )))) == NULL )
	{
		fprintf(ioQQQ," atmdat_readin could not MALLOC TauLines\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now rewind the file so we can read it a second time*/
	if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
	{
		fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}
	i = 1;
	/* level 1 magic number */
	nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of level1.dat  */
	if( ( nelem != 5 ) || ( nelec != 12 ) || ( ion != 19 ) )
	{
		fprintf( ioQQQ, 
			" atmdat_readin: the version of level1.dat is not the current version.\n" );
		fprintf( ioQQQ, 
			" Please obtain the current version from the Cloudy web site.\n" );
		fprintf( ioQQQ, 
			" I expected to find the number 04 11 5 and got %li %li %li instead.\n" ,
			nelem , nelec , ion );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this starts at 1 not 0 since zero is reserved for the
	 * dummy line - we counted number of lines above */
	i = 1;
	while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		if( i >= (nLevel1+1) )
		{
			fprintf(ioQQQ," version conflict.\n");
		}
		/* only look at lines without '#' in first col */
		if( chLine[0] != '#')
		{
			/* >>chng 00 apr 24, prevents crash in PresTotCurrent on call from ContSetIntensity, by PvH */
			EmLineJunk( &TauLines[i] );

			/*sscanf( chLine ,
			"%2s%2li%9li %11f %2li %2li %10g %10g %2li\n",
			chS2,&i2,&i3,&f1,&i4,&i5,&f2,&f3,&i6 );*/

			/* first two cols of line has chem element symbol, 
			 * try to match it with the list of names 
			 * there are no H lines in the level 1 set so zero
			 * is an invalid result */

			/* copy first two char into chS2 and null terminate */
			strncpy( chS2 , chLine , 2);
			chS2[2] = 0;

			ipZ = 0;
			for( j=0; j<LIMELM; ++j)
			{
				if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 )
				{
					/* ipZ is the actual atomic number starting from 1 */
					ipZ = j + 1;
					break;
				}
			}

			/* this happens if no valid chemical element in first two cols on this line */
			if( ipZ == 0 )
			{
				fprintf( ioQQQ, 
					" atmdat_readin could not identify chem symbol on this level 1line:\n");
				fprintf( ioQQQ, "%s\n", chLine );
				fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* now stuff them into the array, will convert over to proper units later */
			TauLines[i].nelem = (int)ipZ;
			/* start the scan early on the line -- the element label will not be
			 * picked up, first number will be the ion stage after the label, as in c  4 */
			j = 1;
			TauLines[i].IonStg = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 1.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/*TauLines[i].IonStg = i2;*/
			TauLines[i].WLAng = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 2.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/*TauLines[i].WLAng = (float)i3;*/
			TauLines[i].EnergyWN = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 3.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/*TauLines[i].EnergyWN = (float)f1;*/
			TauLines[i].gLo = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 4.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/*TauLines[i].gLo = (float)i4;*/
			TauLines[i].gHi = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 5.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/*TauLines[i].gHi = (float)i5;*/

			/* gf is log if negative */
			TauLines[i].gf = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 6.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/*TauLines[i].gf = (float)f2;*/
			if( TauLines[i].gf < 0. )
				TauLines[i].gf = (float)pow(10.f,TauLines[i].gf);

			/* Aul is log if negative */
			TauLines[i].Aul = FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 7.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			if( TauLines[i].Aul < 0. )
				TauLines[i].Aul = pow(10.,TauLines[i].Aul);
			/*TauLines[i].Aul = f3;*/

			TauLines[i].iRedisFun = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " There should have been a number on this level1 line 8.   Sorry.\n" );
				fprintf( ioQQQ, "string==%s==\n" ,chLine );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/*TauLines[i].iRedisFun = i6;*/

			/* this is new in c94.01 and returns nothing (0.) for most lines */
			TauLines[i].cs = FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);

			/* finally increment i */
			++i;
		}
	}

	/* now close the file */
	fclose(ioDATA);

	/*************************************************************
	 *                                                           *
	 * get the level 2 line, opacity project, data set           *
	 *                                                           *
	 *************************************************************/

	/* nWindLine is initialized to the dimension of the vector when it is
	 * initialized in the definition at the start of this file.  
	 * it is set to -1 with the "no level2" command, which
	 * stops us from trying to establish this vector */
	if( nWindLine > 0 )
	{
		/* begin level2 level 2 wind block */
		/* open file with level 2 line data */

		/* create the TauLine2 emline array */
		if( (TauLine2 = ((EmLine *)MALLOC( (size_t)nWindLine*sizeof(EmLine )))) == NULL )
			BadMalloc();
		if( (cs1_flag_lev2 = ((float *)MALLOC( (size_t)nWindLine*sizeof(float )))) == NULL )
			BadMalloc();

		/* first initialize entire array to dangerously large negative numbers */
		for( i=0; i< nWindLine; ++i )
		{
			/* >>chng 99 jul 16, from setting all t[] to flt_max, to call for
			 * following, each member of structure set to own type of impossible value */
			EmLineJunk( &TauLine2[i] );
		}

		/* check on path if file not here and path set */
		/* path was parsed in getset */
		if( lgDataPathSet == TRUE )
		{
			/*path set, so look only there */
			strcpy( chFilename , chDataPath );
			strcat( chFilename , "level2.dat" );
		}
		else
		{
			/* path not set, check local space only */
			strcpy( chFilename , "level2.dat" );
		}

		if( trace.lgTrace )
			fprintf( ioQQQ," atmdat_readin opening level2.dat:");

		if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
		{
			fprintf( ioQQQ, " atmdat_readin could not open level2.dat\n" );
			if( lgDataPathSet == TRUE )
			{
				fprintf( ioQQQ, " atmdat_readin could not open level2.dat, even tried path.\n");
				fprintf( ioQQQ, " path is *%s*\n",chDataPath );
				fprintf( ioQQQ, " final path is *%s*\n",chFilename );
			}
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* get magic number off first line */
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " level2.dat error getting magic number\n" );
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* scan off three numbers, should be the yr mn dy of creation date */
		i = 1;
		nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* the following is the set of numbers that appear at the start of level2.dat 01 05 29 */
		if( lgEOL || ( nelem != 01 ) || ( nelec != 05 ) || ( ion != 29 ) )
		{
			fprintf( ioQQQ, 
				" atmdat_readin: the version of level2.dat is not the current version.\n" );
			fprintf( ioQQQ, 
				" I expected to find the number 01 05 29 and got %li %li %li instead.\n" ,
				nelem , nelec , ion );
			fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
			fprintf( ioQQQ, "Please obtain the correct version.\n");
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now get the actual data */
		i = 0;
		while( i < nWindLine )
		{
			float tt[7];
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			{
				fprintf( ioQQQ, " level2.dat error getting line  %li\n", i );
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* option to skip any line starting with # */
			if( chLine[0]!='#' )
			{
				/*printf("string is %s\n",chLine ) ;*/
				sscanf( chLine , "%f %f %f %f %f %f %f " , 
					&tt[0] ,
					&tt[1] ,
					&tt[2] ,
					&tt[3] ,
					&tt[4] ,
					&tt[5] ,
					&tt[6] );
				/* these are readjusted into their final form in the structure 
				 * in routine lines_setup*/
				TauLine2[i].nelem = (int)tt[0];
				TauLine2[i].IonStg = (int)tt[1];
				TauLine2[i].gLo = tt[2];
				TauLine2[i].gHi = tt[3];
				TauLine2[i].gf = tt[4];
				TauLine2[i].EnergyWN = tt[5];
				cs1_flag_lev2[i] = tt[6];
				++i;
			}
		}
		/* get magic number off last line */
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}
		sscanf( chLine , "%ld" , &magic2 );
		if( 999 != magic2 )
		{
			fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n", 
			  magic2 );
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}
		fclose( ioDATA );
		if( trace.lgTrace )
			fprintf( ioQQQ," OK \n");

		/* end of level2 block*/
	}

	/*************************************************************
	 *                                                           *
	 * get the UTA line set, both date in behar_uta.dat and other*
	 *                                                           *
	 *************************************************************/
	/*# >>refer	Fe	inner shell	Behar, E., Sako, M, Kahn, S.M., 
	# >>refercon	2001, ApJ, 563, 497-504
	# the fits were based on the paper
	# >>refer	Fe	inner shell	Behar, E., & Netzer, H., 2002, ApJ, 570, 165-170 */

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "behar_uta.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "behar_uta.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," atmdat_readin opening behar_uta.dat:");

	if( ( ioBEHAR = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " atmdat_readin could not open behar_uta.dat\n" );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " atmdat_readin could not open behar_uta.dat\n");
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now do exactly same thing for Romas Kissielius data */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "kisielius_uta.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "kisielius_uta.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," atmdat_readin opening kisielius_uta.dat:");

	if( ( ioROMAS = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " atmdat_readin could not open kisielius_uta.dat\n" );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " atmdat_readin could not open kisielius_uta.dat\n");
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* first line is a version number and does not count */
	if( fgets( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
	{
		fprintf( ioQQQ, " atmdat_readin could not read first line of kisielius_uta.dat.\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* count how many lines are in the file, ignoring all lines
	 * starting with '#' */
	nUTA = 0;

	/* first fill in Behar & Netzer interpolation - first pass just count number of lines */
	for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
	{
		for( nelem=ipISO; nelem<LIMELM; ++nelem )
		{
			++nUTA;
		}
	}

	while( fgets( chLine , (int)sizeof(chLine) , ioBEHAR ) != NULL )
	{
		/* we want to count the lines that do not start with #
		 * since these contain data */
		if( chLine[0] != '#')
			++nUTA;
	}

	/* now do same for Romas data */
	for( ion=0; ion<NIONROMAS; ++ion )
	{
		nFeIonRomas[ion] = 0;
	}
	nUTA_Romas = 0;
	while( fgets( chLine , (int)sizeof(chLine) , ioROMAS ) != NULL )
	{
		/* we want to count the lines that do not start with #
		 * since these contain data */
		if( chLine[0] != '#')
		{
			++nUTA_Romas;
			i = 1;
			ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( ion<0 || ion>=ipIRON )
			{
				fprintf(ioQQQ," impossible ion for romas data\n");
			}
			ASSERT( ion < NIONROMAS );
			++nFeIonRomas[ion];
		}
	}

	/* now malloc the UTA lines array with enough space for both
	 * data sets */
	if( (UTALines = ((EmLine *)MALLOC( (size_t)(nUTA+nUTA_Romas)*sizeof(EmLine )))) == NULL )
		BadMalloc();

	/* this will count number of lines, must equal nUTA at end */
	i = 0;
	/* first fill in Behar & Netzer interpolation - now do real thing,
	 * these are Nick Abel's fits to the BN data */
	for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
	{
		for( nelem=ipISO; nelem<LIMELM; ++nelem )
		{
			double f1;
			double sum_spon_auto;
#			define N 7

			double alam[N]={ 1.859035 , 4.692138 ,  10.324543 , 19.294364 ,
				40.401292  ,  44.442754 , 72.959719 };
			double blam[N]={-9.7855968,-11.159739, -12.914931 , -14.987272,
				-18.853446 , -19.271781 ,-23.828388 };
			double clam[N]={ 8.2874628, 8.3043002,  8.3164038 ,  8.3269937,
				8.4312895  , 8.3730893  , 8.493802 };
			/* >>chng 05 mar 07, by NA, update fits to better form */

			double aA[N] = {1.9067 , 1.8715 ,  2.1033, 2.2815 ,
				1.9511 , 1.9966 ,  2.0852};
			double bA[N] = {8.4538 , 8.5528, 7.616, 6.9247 ,
				7.9479, 7.9777 , 7.8349 };
		
			double aDep[N] = {29.54 , 12.07 , 24.23 , 7.35 , 7.37 , 11.14 , 11.99 };
			double bDep[N] = {-14.853 , -7.606 , -7.844 ,  9.987 , 10.503 , 8.541 , 8.865 };

			/* NB - these are plus one is because 1 will be subtracted when used */
			/* derive element and ion */
			UTALines[i].nelem = nelem+1;
			ion = nelem + 1 - ipISO;
			UTALines[i].IonStg = ion;

			/* these were the statistical weights */
			/*TODO	2	- check these */
			UTALines[i].gHi = 4.f;
			UTALines[i].gLo = 2.f;

			f1 = alam[ipISO-2]*1e-4 + blam[ipISO-2]*1e-4*(nelem+1) + 
				clam[ipISO-2]*1e-4*POW2(nelem+1.);
			f1 = (1./f1);
			UTALines[i].WLAng = (float)f1;
			UTALines[i].EnergyWN = (float)(1e8/f1);

			f1 = aA[ipISO-2]*log((double)(nelem+1)) + bA[ipISO-2] ;

			UTALines[i].Aul = pow(10., f1 );
			/* generate gf from A */
			UTALines[i].gf = 
				(float)(UTALines[i].Aul*UTALines[i].gHi*
				1.4992e-8*POW2(1e4/UTALines[i].EnergyWN));

			UTALines[i].iRedisFun = 1;

			/* find sum of spontaneous relax and autoionization As */
			f1 = log((double)(nelem+1));
			f1 = POW2( f1 );
			sum_spon_auto = aDep[ipISO-2]*1e-2*f1 + bDep[ipISO-2]*1e-1 ;

			if( ipISO == ipBE_LIKE )
			{
				sum_spon_auto = pow( E , sum_spon_auto );
				sum_spon_auto = pow(10., sum_spon_auto )*1e13;
			}
			else if( ipISO <= ipF_LIKE )
			{
				sum_spon_auto = pow(10., sum_spon_auto )*1e13;
			}
			else
				TotalInsanity();

			/* last number is negative branching ratio for autoionization, 
			 * which we will store as negative of cs so not misinterpreted as cs */
			/* >>chng 03 feb 18, change to negative of br */

			/*ASSERT( UTALines[i].Aul / sum_spon_auto <= 1. && UTALines[i].Aul / sum_spon_auto >= 0. );*/
			UTALines[i].cs = MIN2(0., UTALines[i].Aul / sum_spon_auto - 1. );

			/* option to print UTAs */
#			if 0
			fprintf(ioQQQ," %2s%2li\t%.4f\t%.3e\t%.3e\t%.3e\n",
				elementnames.chElementSym[nelem], 
				ion, 
				UTALines[i].WLAng ,
				UTALines[i].Aul ,
				sum_spon_auto , 
				-UTALines[i].cs);
#			endif

			/* finally increment i */
			++i;
#			undef N
		}
	}
	/*exit(1);*/
	/* now rewind the file so we can read it a second time*/
	if( fseek( ioBEHAR , 0 , SEEK_SET ) != 0 )
	{
		fprintf( ioQQQ, " atmdat_readin could not rewind behar_uta.dat.\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* now rewind the file so we can read it a second time*/
	if( fseek( ioROMAS , 0 , SEEK_SET ) != 0 )
	{
		fprintf( ioQQQ, " atmdat_readin could not rewind behar_uta.dat.\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
	{
		fprintf( ioQQQ, " atmdat_readin could not read first line of behar_uta.dat.\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}
	j = 1;
	/* behar_uta.dat magic number */
	nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
	nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
	ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of behar_uta.dat 2002 8 19 */
	/* >>chng 02 oct 22, had incorrect stat weights for all stages of ionization except atom.
	 * no longer used stored values when stat weight is zero */
	if( ( nelem != 2002 ) || ( nelec != 10 ) || ( ion != 22 ) )
	{
		fprintf( ioQQQ, 
			" atmdat_readin: the version of behar_uta.dat is not the current version.\n" );
		fprintf( ioQQQ, 
			" I expected to find the number 2002 10 22 and got %li %li %li instead.\n" ,
			nelem , nelec , ion );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read in the data - we already know that there are nUTA of them */
	for( nDataSet=0; nDataSet<2; ++nDataSet )
	{
		if( nDataSet==0 )
		{
			ioDATA = ioBEHAR;
		}
		else if( nDataSet==1 )
		{
			ioDATA = ioROMAS;
		}
		else
			TotalInsanity();

		while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
		{
			/* only look at lines without '#' in first col */
			if( chLine[0] != '#')
			{
				long int i1, i2, i3;
				double f1, f2, oscill;
				double frac_relax;
				EmLineJunk( &UTALines[i] );

				sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
					&i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );

				/* all these are iron, first number was ion stage with 0 the atom */
				UTALines[i].nelem = ipIRON+1;

				/* now do stage of ionization */
				ASSERT( i1 >= 0 && i1 < ipIRON );
				/* NB - the plus one is because 1 will be subtracted when used,
				* in the original data file i1 is 0 for the atom */
				UTALines[i].IonStg = i1 + 1;

				/* these were the statistical weights */
				if( i2>0 && i3>0 )
				{
					UTALines[i].gHi = (float)i2;
					UTALines[i].gLo = (float)i3;
				}
				else
				{
					/* nearly all stages of ionization do not include the stat weights,
					* so use unity */
					UTALines[i].gHi = 1.f;
					UTALines[i].gLo = 1.f;
				}

				/*TauLines[i].IonStg = i2;*/
				UTALines[i].WLAng = (float)f1;

				UTALines[i].EnergyWN = (float)(1e8/f1);

				/* this is true spontaneous rate for doubly excited state to inner shell,
				* and is used for pumping, and also relaxing back to inner shell */
				UTALines[i].Aul = f2;
				/* generate gf from A */
				/* >>chng 02 oct 22, use absorption oscillator strength in data file
				* since don't have stat weights for most lines */
				UTALines[i].gf = UTALines[i].gLo * (float)oscill;

				UTALines[i].iRedisFun = 1;

				/* last number is branching ratio for autoionizing, 
				* which we will store in negative in cs so not misinterpreted as real cs*/
				/* >>chng 03 feb 18, cs to negative of br */
				ASSERT( frac_relax >=0. &&  frac_relax <=1. );
				UTALines[i].cs = frac_relax-1.;
	#			if 0
				fprintf(ioQQQ,"data set %li line %li %2s%2li\t%.4f\t%.3e\t%.3e\n",
					nDataSet,
					i,
					elementnames.chElementSym[UTALines[i].nelem-1], 
					UTALines[i].IonStg, 
					UTALines[i].WLAng ,
					UTALines[i].Aul ,
					-UTALines[i].cs);
	#			endif

				/* finally increment i */
				++i;
			}
		}
	}

	/* these have to agree */
	ASSERT( i == (nUTA+nUTA_Romas) );

	/* this is option to use the data Romas computed.  */
	if( ionbal.lgInnerShell_Romas )
	{
		for( i=0; i<nUTA; ++i )
		{
			if( UTALines[i].nelem-1 == ipIRON )
			{
				/* nFeIonRomas counts how many of the Romas lines
				 * fall in this stage of ionization - if zero then
				 * none, so keep Behar data, if positive, then
				 * there are new data, and zero out Behar data */
				ASSERT( UTALines[i].IonStg-1< NIONROMAS );
				if( nFeIonRomas[UTALines[i].IonStg-1] )
				{
					UTALines[i].Aul = -1.;
				}
			}
		}
		nUTA += nUTA_Romas;
	}

	/* now close the files */
	fclose(ioBEHAR);
	fclose(ioROMAS);

	/* here there are two options - we can add the Romas data to the
	 * first block of UTA data, but zero out any overlapping data,
	 * or simply ignore the Romas data */

	/* allocate space for the carbon monoxide CO rotation levels */
	ASSERT( nCORotate > 0);
	if( (C12O16Rotate = ((EmLine *)MALLOC( (size_t)nCORotate*sizeof(EmLine )))) == NULL )
	{
		fprintf(ioQQQ," atmdat_readin could not MALLOC 1 C12O16Rotate\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( (C13O16Rotate = ((EmLine *)MALLOC( (size_t)nCORotate*sizeof(EmLine )))) == NULL )
	{
		fprintf(ioQQQ," atmdat_readin could not MALLOC 1 C13O16Rotate\n");
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this says we cannot change the number of levels again, since space is allocated */
	lgCORotateMalloc = TRUE;

	/* next initialize entire array to dangerously large negative numbers */
	for( J=0; J< nCORotate; ++J )
	{
		EmLineJunk( &C12O16Rotate[J] );
		EmLineJunk( &C13O16Rotate[J] );
	}

	/* read in data for the set of hyperfine structure lines, and allocate
	 * space for the EmLine HFLines[nHFLines] structure */
	HyperfineCreate();

	/* initialize the large blcok of level 1 real lines, and OP level 2 lines */
	lines_setup();

	/*************************************************************
	 *                                                           *
	 * get the Auger electron yield data set                     *
	 *                                                           *
	 *************************************************************/

	/* mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat ========*/

	/* open file with Mewe coefficients */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "mewe_gbar.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "mewe_gbar.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," atmdat_readin opening mewe_gbar.dat:");

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " could not open mewe_gbar.dat\n" );
			if( lgDataPathSet == TRUE )
			{
				fprintf( ioQQQ, " atmdat_readin could not open mewe_gbar.dat, even tried path.\n");
				fprintf( ioQQQ, " path is *%s*\n",chDataPath );
				fprintf( ioQQQ, " final path is *%s*\n",chFilename );
			}
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* get magic number off first line */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " mewe_gbar.dat error getting magic number\n" );
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* check the number */
	sscanf( chLine , "%ld" , &magic1 );
	if( magic1 != 9101 )
	{
		fprintf( ioQQQ, " mewe_gbar.dat starts with wrong magic number=%ld \n", 
		  magic1 );
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now get the actual data, indices are correct for c, in Fort went from 2 to 210 */
	for( i=1; i < 210; i++ )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " mewe_gbar.dat error getting line  %li\n", i );
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}
		/*printf("%s\n",chLine);*/
		sscanf( chLine , "%f %f %f %f " , 
			&MeweCoef.g[i][0] ,
			&MeweCoef.g[i][1] ,
			&MeweCoef.g[i][2] ,
			&MeweCoef.g[i][3] );
	}

	/* get magic number off last line */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" );
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	sscanf( chLine , "%ld" , &magic2 );

	if( magic1 != magic2 )
	{
		fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n", 
		  magic2 );
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	fclose( ioDATA );

	if( trace.lgTrace )
		fprintf( ioQQQ," OK \n");

	 /* preset two arrays that will hold auger data 
	  * set to very large values so that code will blow
	  * if not set properly below */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		for( ion=0; ion < LIMELM; ion++ )
		{
			Heavy.nsShells[nelem][ion] = 100000000;
			for( ns=0; ns < 7; ns++ )
			{
				yield.n_elec_eject[nelem][ion][ns] = 1000000000;
				for( nelec=0; nelec < 10; nelec++ )
				{
					yield.frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX;
				}
			}
		}
	}

	/* hydrogen and helium will not be done below, so set yields here*/
	yield.n_elec_eject[ipHYDROGEN][0][0] = 1;
	yield.n_elec_eject[ipHELIUM][0][0] = 1;
	yield.n_elec_eject[ipHELIUM][1][0] = 1;

	yield.frac_elec_eject[ipHYDROGEN][0][0][0] = 1;
	yield.frac_elec_eject[ipHELIUM][0][0][0] = 1;
	yield.frac_elec_eject[ipHELIUM][1][0][0] = 1;

	/* open file auger.dat, yield data file that came from
	 * >>refer	all	auger	Kaastra, J.S., & Mewe, R. 1993, A&AS, 97, 443-482 */
	/* check on path is file not here and path set */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "mewe_nelectron.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "mewe_nelectron.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," atmdat_readin opening mewe_nelectron.dat:");

	/* open the file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " could not open auger.dat\n" );
		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " atmdat_readin could not open auger.dat, even on path..\n");
			fprintf( ioQQQ, " path is *%s*\n",chDataPath );
			fprintf( ioQQQ, " final path is *%s*\n",chFilename );
		}
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read in all auger yields
	 * will do elements from li on up,
	 * skip any line starting with #
	 * this loop goes from lithium to Zn */
	for( nelem=2; nelem < LIMELM; nelem++ )
	{
		/* nelem is on the shifted C scale, so 2 is Li */
		for( ion=0; ion <= nelem; ion++ )
		{
			/* number of bound electrons, = atomic number for neutral */
			nelec = nelem - ion + 1;
			/* one of dima's routines to determine the number of electrons
			 * for this species, nelem +1 to shift to physical number */
			/* subroutine atmdat_outer_shell(iz,in,imax,ig0,ig1)
			 * iz - atomic number from 1 to 30 (integer) 
			 * in - number of electrons from 1 to iz (integer)
			 * Output: imax - number of the outer shell
			 */
			atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1);
			/* loop over all the electrons */
		   /* nsShells(nelem,ion) is outer shell number for ion with nelec electrons
		    * on physics scale, with K shell being 1 */
			Heavy.nsShells[nelem][ion] = imax;
			for( ns=0; ns < imax; ns++ )
			{
				char ch1 = '#';
				/* the * is the old comment char, accept it, but really want # */
				while( ch1 == '#' || ch1 == '*' )
				{
					if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
					{
						fprintf( ioQQQ, " auger.dat error getting line  %li\n", ns );
						puts( "[Stop in atmdat_readin]" );
						cdEXIT(EXIT_FAILURE);
					}
					ch1 = chLine[0];
				}
				/*printf("string is %s\n",chLine ) ;*/
				sscanf( chLine , "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f" , 
					&temp[0] ,&temp[1] ,&temp[2] ,&temp[3] ,&temp[4] ,
					&temp[5] ,&temp[6] ,&temp[7] ,&temp[8] ,&temp[9] ,
					&temp[10] ,&temp[11] ,&temp[12] ,&temp[13] ,&temp[14] );
				yield.n_elec_eject[nelem][ion][ns] = (long int)temp[3];

				/* can pop off up to 10 auger electrons, these are the probabilities*/
				for( j=0; j<10; j++)
				{
					yield.frac_elec_eject[nelem][ion][ns][j] = temp[j+4];
				}
			}
		}
		/* activate this print statement to get yield for k-shell of all atoms */
		/*fprintf(ioQQQ,"yyyield\t%li\t%.4e\n", nelem+1, yield.frac_elec_eject[nelem][0][0][0] );*/
	}

	fclose( ioDATA );

	/* validate the data, this is a nested set of asserts that should be
	 * optimized away in fast mode */
	for( nelem=2; nelem < LIMELM; nelem++ )
	{
		for( ion=0; ion <= nelem; ion++ )
		{
			ASSERT( Heavy.nsShells[nelem][ion] > 0 );
			ASSERT( Heavy.nsShells[nelem][ion] <= 10 );
			for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
			{
				ASSERT( yield.n_elec_eject[nelem][ion][ns] >= 0 );
				ASSERT( yield.n_elec_eject[nelem][ion][ns] < 11 );
				for( i=0; i < yield.n_elec_eject[nelem][ion][ns]; i++ )
				{
					ASSERT( yield.frac_elec_eject[nelem][ion][ns][i] >= 0. );
				}
			}
		}
	}

	/* this is set true with no auger command */
	if( yield.lgKillAuger )
	{
		int nstage , nshell;
		if( trace.lgTrace )
			fprintf( ioQQQ," Auger data set zeroed. \n");
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			for( nstage=0; nstage < nelem; nstage++ )
			{
				for( nshell=0; nshell < 7; nshell++ )
				{
					/* only one electron comes out */
					yield.n_elec_eject[nelem][nstage][nshell] = 1;
					yield.frac_elec_eject[nelem][nstage][nshell][0] = 1.;
				}
			}
		}
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," OK \n");

	/* open file mewe_fluor.dat, yield data file that came from
	 * >>refer	all	auger	Kaastra, J.S., & Mewe, R. 1993, 97, 443-482 */
	/* check on path is file not here and path set */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "mewe_fluor.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "mewe_fluor.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," atmdat_readin opening mewe_fluor.dat:");

	/* open the file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " could not open mewe_fluor.dat\n" );
		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " atmdat_readin could not open mewe_fluor.dat, even on path..\n");
			fprintf( ioQQQ, " path is *%s*\n",chDataPath );
			fprintf( ioQQQ, " final path is *%s*\n",chFilename );
		}
		puts( "[Stop in atmdat_readin]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read in all auger yields
	 * will do elements from li on up,
	 * skip any line starting with # */
	do 
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " mewe_fluor.dat error getting line  %i\n", 0 );
			puts( "[Stop in atmdat_readin]" );
			cdEXIT(EXIT_FAILURE);
		}
	} while( chLine[0] == '#' );

	lgEOL = FALSE;

	yield.nfl_lines = 0;
	do
	{
#		define NKM	10
		int nDima[NKM] = {0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
		int nAuger;
		/*printf("string is %s\n",chLine ) ;*/
		sscanf( chLine , "%f %f %f %f %f %f %f" , 
			&temp[0] ,&temp[1] ,&temp[2] ,&temp[3] ,&temp[4] ,
			&temp[5] ,&temp[6]  );

		/* the atomic number, C is 5 */
		yield.nfl_nelem[yield.nfl_lines] = (int)temp[0]-1;
		ASSERT(yield.nfl_nelem[yield.nfl_lines]>=0 && yield.nfl_nelem[yield.nfl_lines]<LIMELM);

		/* the ion stage for target, atom is 0 */
		yield.nfl_ion[yield.nfl_lines] = (int)temp[1]-1;
		ASSERT(yield.nfl_ion[yield.nfl_lines]>=0 && yield.nfl_ion[yield.nfl_lines]<=yield.nfl_nelem[yield.nfl_lines]+1);

		/* the target's shell */
		yield.nfl_nshell[yield.nfl_lines] = nDima[(long)temp[2]-1];
		ASSERT(yield.nfl_nshell[yield.nfl_lines]>=0 && yield.nfl_nshell[yield.nfl_lines]<
			/* nsShells is shell number, where K is 1 */
			(Heavy.nsShells[yield.nfl_nelem[yield.nfl_lines]][yield.nfl_ion[yield.nfl_lines]]-1));

		/* this is the number of Auger electrons ejected */
		nAuger = (int)temp[3];
		/* so this is the spectrum of the photons that are emitted */
		yield.nfl_ion_emit[yield.nfl_lines] = yield.nfl_ion[yield.nfl_lines] + nAuger + 1;
		/* must be gt 0 since at least photoelectron comes off */
		ASSERT(yield.nfl_ion_emit[yield.nfl_lines]>0 && yield.nfl_ion_emit[yield.nfl_lines]<=yield.nfl_nelem[yield.nfl_lines]+1);

		/* this is the type of line as defined in their paper */
		yield.nfl_nLine[yield.nfl_lines] = (int)temp[4];

		/* energy in Ryd */
		yield.fl_energy[yield.nfl_lines] = temp[5] / (float)EVRYD;
		ASSERT( yield.fl_energy[yield.nfl_lines] > 0. );

		/* fluor yield */
		yield.fl_yield[yield.nfl_lines] = temp[6];
		/* NB cannot assert <=1 since data file has yields around 1.3 - 1.4 */
		ASSERT(yield.fl_yield[yield.nfl_lines]>=0 );

		++yield.nfl_lines;
		if( yield.nfl_lines > MEWE_FLUOR )
			TotalInsanity();
		do
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				lgEOL = TRUE;
		} while( chLine[0]=='#'  && !lgEOL );
	} while(  !lgEOL );

	fclose( ioDATA );

	/****************************************************************
	 *                                                              *
	 * get the Hummer and Storey model case A, B results, these are *
	 * the two data files e1b.dat and e2b.dat, for H and He         *
	 *                                                              *
	 ****************************************************************/

	/* now get Hummer and Storey case b data, loop over H, He */
	/* >>chng 01 aug 08, generalized this to both case A and B, and all elements
	 * up to HS_NZ, now 8 */
	for( ipZ=0; ipZ<HS_NZ; ++ipZ )
	{
		/* don't do the minor elements, Li-B */
		if( ipZ>1 && ipZ<5 ) continue;

		for( iCase=0; iCase<2; ++iCase )
		{
			/* open file with case B data */
			/* check on path is file not here and path set */
			if( lgDataPathSet == TRUE )
			{
				/*path set, so look only there */
				strcpy( chFilename , chDataPath );
				/* >>chng 01 aug 08, add HS_ to start of file names to indicate Hummer Storey origin */
				strcat( chFilename , "HS_e" );
			}
			else
			{
				/* path not set, check local space only */
				strcpy( chFilename , "e" );
			}

			/* create remainder of file name for this charge
			 * first character after e is charge number for this element,
			 * stuff this into temporary work array */
			chLine[0] = (char)(ipZ + '1'); 

			/* terminate it*/
			chLine[1] = '\0' ; 

			/* add it to end of file name */
			strcat( chFilename , chLine );

			if( iCase == 0 )
			{
				/* add case a to end of name */
				strcat( chFilename , "a" );
			}
			else if (iCase == 1 )
			{
				/* add case b to end of name */
				strcat( chFilename , "b" );
			}

			/* finally end the name with .dat*/
			strcat( chFilename , ".dat") ;

			if( trace.lgTrace )
				fprintf( ioQQQ," atmdat_readin opening Hummer Storey emission file %s:",chFilename);

			/* open the file */
			if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
			{
				fprintf( ioQQQ, " could not open HS case b %s\n",chFilename);
				if( lgDataPathSet == TRUE )
				{
					fprintf( ioQQQ, " atmdat_readin could not open even on path..\n");
					fprintf( ioQQQ, " path is *%s*\n",chDataPath );
					fprintf( ioQQQ, " final path is *%s*\n",chFilename );
				}
				puts( "[Stop in atmdat_readin]" );
				cdEXIT(EXIT_FAILURE);
			}
			if( trace.lgTrace )
			{
				fprintf( ioQQQ," OK\n");
			}

			/* read in the number of temperatures and densities*/
			i = fscanf( ioDATA, "%li %li ", 
				&atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] );

			/* check that ntemp and nDensity are below NHSDIM, 
			 * set to 15 in atmdat_HS_caseB.h */
			assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM );
			assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM) ;

			/* loop reading in line emissivities for all temperatures*/
			for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ )
			{
				for( ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ )
				{
					long int junk, junk2 , ne;
					fscanf( ioDATA, " %lf %li %lf %c %li %ld ", 
						&atmdat.Density[iCase][ipZ][ipDens], &junk , 
						&atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 , 
						&atmdat.ncut[iCase][ipZ] );

					ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2;
					ASSERT( ne<=NLINEHS );
					for( j=0; j < ne; j++ )
					{
						fscanf( ioDATA, "%lf ", &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] );
					}
				}
			}

			/*this is end of read-in loop */
			fclose(ioDATA); 
#			if 0
			/* print list of densities and temperatures */
			for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ )
			{
				fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]) ;
			}
			fprintf(ioQQQ,"\n") ;
			for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ )
			{
				fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]) ;
			}
			fprintf(ioQQQ,"\n") ;
#			endif
		}
	}

	/* Verner's model FeII atom
	 * do not call if Netzer model used, or it Fe(2) is zero
	 * exception is when code is searching for first soln */
	/* read the atomic data from files */
	/* convert line arrays to internal form needed by code */
	FeIICreate(  );

#	ifdef DEBUG_FUN
	fputs( " <->atmdat_readin()\n", debug_fp );
#	endif
	return;
}

