/* 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 */
/*CoolIron compute iron cooling */
/*Fe_10_11_13_cs evaluate collision stength for Fe */
/*fe14cs compute collision strengths for forbidden transitions */
/*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */
/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
/*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
/*Fe2_cooling compute cooling due to FeII emission */
/*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */
/*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */
#include "cddefines.h"
#include "physconst.h"
#include "cpu.h"
#include "dense.h"
#include "coolheavy.h"
#include "taulines.h"
#include "phycon.h"
#include "iso.h"
#include "conv.h"
#include "trace.h"
#include "hydrogenic.h"
#include "ligbar.h"
#include "cooling.h"
#include "thermal.h"
#include "lines_service.h"
#include "atoms.h"
#include "atomfeii.h"
#include "fe.h"
#define	NLFE2	6

/*lint -e777 float tests equality - several places old temp saved and compared with current */
/*Fe2Lev16 compute populations and cooling due to 16 level Fe II ion
static void Fe2Lev16(void); */

/*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */
static void Fe11Lev5(void);

/*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */
static void Fe13Lev5(void);

/*fe14cs compute collision strengths for forbidden transitions */
static void fe14cs(double te1, 
	  double *csfe14);

/*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
static void Fe7Lev8(void);

/*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion */
static void Fe3Lev14(void);

/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
static void Fe4Lev12(void);

/*Fe_10_11_13_cs evaluate collision stength for Fe */
static double Fe_10_11_13_cs(
	/* ion, one of 10, 11, or 13 on physics scale */
	int ion, 
	/* inital and final index of levels, with lowest energy 1, highest of 5 
	 * on physics scale */
	int init, 
	int final )
{
#	define N 10
	static double Fe10cs[6][6][2];
	static double Fe11cs[6][6][2];
	static double Fe13cs[6][6][2];
	int i, j , k;
	double cs;
	int index = 0;
	double temp_max, temp_min = 4;
	double temp_log = phycon.alogte;
	static int lgFirstTime = TRUE;

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

	if( lgFirstTime )
	{
		/* fit coefficients for collision strengths */
		double aFe10[N] = {10.859,-1.1541,11.593,22.333,-0.4283,7.5663,3.087,1.0937,0.8261,59.678};
		double bFe10[N] = {-1.4804,0.4956,-2.1096,-4.1288,0.1929,-1.3525,-0.5531,-0.1748,-0.1286,-11.081};
		double aFe11[N] = {5.7269,1.2885,4.0877,0.4571,1.2911,2.2339,0.3621,0.7972,0.2225,1.1021};
		double bFe11[N] = {-0.7559,-0.1671,-0.5678,-0.0653,-0.1589,-0.2924,-0.0506,-0.1038,-0.0302,-0.1062};
		double aFe13[N] = {2.9102,1.8682,-0.353,0.0622,14.229,-4.3845,0.0375,-6.9222,0.688,-0.0609};
		double bFe13[N] = {-0.4158,-0.242,0.1417,0.0023,-2.0643,1.2573,0.0286,2.0919,-0.083,0.1487};
		/* do one-time initialization 
		 * assigning array initially to zero */ 
		for (i=0; i<6; i++)
		{
			for (j=0; j<6; j++)
			{
				for( k=0; k<1; ++k )
				{
					Fe10cs[i][j][k] = 0;
					Fe11cs[i][j][k] = 0;
					Fe13cs[i][j][k] = 0;
				}
			}
		}
		/* this is done for safety */
		i = 0;
		for( j=0; j<6; ++j )
		{
			for( k=0; k<1; ++k )
			{
				SET_DBL_SNAN( Fe10cs[i][j][k] );
				SET_DBL_SNAN( Fe11cs[i][j][k] );
				SET_DBL_SNAN( Fe13cs[i][j][k] );
			}
		}
		j = 0;
		for( i=0; i<6; ++i )
		{
			for( k=0; k<1; ++k )
			{
				SET_DBL_SNAN( Fe10cs[i][j][k] );
				SET_DBL_SNAN( Fe11cs[i][j][k] );
				SET_DBL_SNAN( Fe13cs[i][j][k] );
			}
		}
		/* reading coeffs into 3D array */
		for (i=1; i<6; i++)
		{
			for (j=i+1; j<6; j++)
			{
				Fe10cs[i][j][0] = aFe10[index];
				Fe10cs[i][j][1] = bFe10[index];
				Fe11cs[i][j][0] = aFe11[index];
				Fe11cs[i][j][1] = bFe11[index];
				Fe13cs[i][j][0] = aFe13[index];
				Fe13cs[i][j][1] = bFe13[index];
				index++;
			}
		}
		lgFirstTime = FALSE;
	}

	/* Invalid entries returns '-1':the initial indices are smaller than the final indices */
	if (init >= final)
	{
		cs = -1;
	}
	/* Invalid returns '-1': the indices are greater than 5 or smaller than 0 */
	else if (init < 1 || init > 4 || final < 2 || final > 5)
	{
		cs = -1;
	}
	else
	{
		/* cs = a + b*log(T) 
		 * if temp is out of range, return boundary values */
		if (ion == 10)
		{
			temp_max = 5;
			temp_log = MAX2(temp_log, temp_min);
			temp_log = MIN2(temp_log, temp_max);
			cs = Fe10cs[init][final][0] + Fe10cs[init][final][1]*temp_log;
		}
		else if (ion == 11)
		{
			temp_max = 6.7;
			temp_log = MAX2(temp_log, temp_min);
			temp_log = MIN2(temp_log, temp_max);
			cs = Fe11cs[init][final][0] + Fe11cs[init][final][1]*temp_log;
		}
		else if (ion ==13)
		{
			temp_max = 5;
			temp_log = MAX2(temp_log, temp_min);
			temp_log = MIN2(temp_log, temp_max);
			cs = Fe13cs[init][final][0] + Fe13cs[init][final][1]*temp_log;
		}
		else
			/* this can't happen */
			TotalInsanity();
	}

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

	return cs;

#	undef N
}

/*Fe2_cooling compute cooling due to FeII emission */
static void Fe2_cooling( void )
{
	long int i , j;
	static double **AulPump,
		**CollRate,
		**AulEscp,
		**col_str ,
		**AulDest, 
		*depart,
		*pops,
		*destroy,
		*create;

	static int lgFirst=TRUE;
	int lgNegPop , 
		lgZeroPop;

	/* stat weights for Fred's 6 level model FeII atom */
	static double gFe2[NLFE2]={1.,1.,1.,1.,1.,1.};
	/* excitation energies (Kelvin) for Fred's 6 level model FeII atom */
	static double ex[NLFE2]={0.,3.32e4,5.68e4,6.95e4,1.15e5,1.31e5};

	/* these are used to only evaluated feii one time per temperature, zone, and abundance */
	static double TUsed = 0.; 
	static float AbunUsed = 0.;
	/* remember which zone last evaluated with */
	static long int nZUsed=-1,
		/* make sure at least two calls per zone */
		nCall=0;

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

	/* return if nothing to do */
	if( dense.xIonDense[ipIRON][1] == 0. )
	{
		/* zero abundance, do nothing */
		/* heating cooling and deriv from large atom */
		FeII.Fe2_large_cool = 0.;
		FeII.Fe2_large_heat = 0.;
		FeII.ddT_Fe2_large_cool = 0.;

		/* cooling and deriv from simple UV atom */
		FeII.Fe2_UVsimp_cool = 0.;
		FeII.ddT_Fe2_UVsimp_cool = 0.;

		/* now zero out intensities of all feii lines and level populations */
		FeIIIntenZero();

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

		return;
	}

	/* this can be the large 371 level FeII atom
	 * >>chng 05 dec 04, always call this but with default of 16 levels only 
	 * >>chng 00 jan 06, total rewrite mostly done 
	 * >>chng 97 jan 17, evaluate large FeII atom cooling every time temp changes
	 * >>chng 97 jan 31, added check on zone since must reeval even const temp
	 * >>chng 99 may 21, reeval when zone or temperature changes, but not when
	 * abundance changes, since we can simply rescale cooling */

	/* totally reevaluate large atom if new zone, or cooling is significant
	 * and temperture changed, we are in search phase,
	 * lgSlow option set true with atom feii slow, forces constant
	 * evaluation of atom */
	if( FeII.lgSlow || 
		nzone < 1 ||
		nzone !=nZUsed || 
		/* on new call, nCall is now set at previous zone's number of calls.
		 * it is set to zero below, so on second call, nCall is 0.  On 
		 * third call nCall is 1.  Check for <1 is to insure at least two calls */
		nCall < 1 ||
		/* check whether things have changed on later calls */
		( phycon.te != TUsed &&  fabs(FeII.Fe2_large_cool/thermal.ctot)> 0.002 &&  
		fabs(dense.xIonDense[ipIRON][1]-AbunUsed)/SDIV(AbunUsed)> 0.002 ) ||
		(phycon.te != TUsed &&  fabs(FeII.Fe2_large_cool/thermal.ctot)> 0.01) )
	{

		if( nZUsed == nzone )
		{
			/* not first call, increment, check above to make sure at least
			 * two evaluations */
			++nCall;
		}
		else
		{
			/* first call this zone set nCall to zero*/
			nCall = 0;
		}

		/* option to trace convergence and feii calls */
		if( trace.lgTrConvg>=5 )
		{
			fprintf( ioQQQ, "        CoolIron5 calling FeIILevelPops since ");
			if( phycon.te!= TUsed )
			{
				fprintf( ioQQQ, 
					"temperature changed, old new are %g %g, nCall %li ", 
					TUsed, phycon.te , nCall);
			}
			else if( nzone !=nZUsed )
			{
				fprintf( ioQQQ, 
					"new zone, nCall %li ", nCall );
			}
			else if( FeII.lgSlow )
			{
				fprintf( ioQQQ, 
					"FeII.lgSlow set  %li", nCall );
			}
			else if( conv.lgSearch )
			{
				fprintf( ioQQQ, 
					" in search phase  %li", nCall );
			}
			else if( nCall < 2 )
			{
				fprintf( ioQQQ, 
					"not second nCall %li " , nCall );
			}
			else if( phycon.te!= TUsed &&  FeII.Fe2_large_cool/thermal.ctot> 0.001 )
			{
				fprintf( ioQQQ, 
					"temp or cooling changed, new are %g %g nCall %li ", 
					phycon.te, FeII.Fe2_large_cool, nCall );
			}
			else
			{
				fprintf(ioQQQ, "????");
			}
			fprintf(ioQQQ, "\n");
		}

		/* remember parameters for current conditions */
		TUsed = phycon.te;
		AbunUsed = dense.xIonDense[ipIRON][1];
		nZUsed = nzone;

		/* this print turned on with atom feii print command */
		if( FeII.lgPrint )
		{
			fprintf(ioQQQ,
				" FeIILevelPops called zone %4li te %5f abun %10e c(fe/tot):%6f nCall %li\n", 
				nzone,phycon.te,AbunUsed,FeII.Fe2_large_cool/thermal.ctot,nCall);
		}

		/* this solves the multi-level problem, 
		 * sets FeII.Fe2_large_cool, 
				FeII.Fe2_large_heat, & 
				FeII.ddT_Fe2_large_cool 
				but does nothing with them */
		FeIILevelPops();
		{
			/*@-redef@*/
			enum{DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC && iteration>1 && nzone>=4/**/ )
			{
				/*Fe2Lev16();*/
				fprintf(ioQQQ,"DEBUG1\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
					nzone,
					phycon.te,
					dense.gas_phase[ipHYDROGEN],
					dense.eden,
					FeII.Fe2_large_cool , 
					FeII.ddT_Fe2_large_cool ,
					FeII.Fe2_large_cool/dense.eden/dense.gas_phase[ipHYDROGEN] , 
					thermal.ctot )
					;
			}
		}

		if( trace.lgTrConvg>=5 || FeII.lgPrint)
		{
			fprintf( ioQQQ, " FeIILevelPops returned cool=%.2e heat=%.2e deriv=%.2e\n ",
					FeII.Fe2_large_cool,FeII.Fe2_large_heat ,FeII.ddT_Fe2_large_cool);
		}

	}
	else if( dense.xIonDense[ipIRON][1]!=AbunUsed )
	{
		float ratio;
		/* this branch, same zone and temperature, but small change in abundance, so just
		 * rescale cooling and deriv by this change.  assumption is that very small changes
		 * in abundance occurs as ots rates damp out */
		if( trace.lgTrConvg>=5 )
		{
			fprintf( ioQQQ, 
				"       CoolIron rescaling FeIILevelPops since small change, CFe2=%.2e CTOT=%.2e\n",
				FeII.Fe2_large_cool,thermal.ctot);
		}
		ratio = dense.xIonDense[ipIRON][1]/AbunUsed;
		FeII.Fe2_large_cool *= ratio;
		FeII.ddT_Fe2_large_cool *= ratio;
		FeII.Fe2_large_heat *= ratio;
		AbunUsed = dense.xIonDense[ipIRON][1];
	}
	else
	{
		/* this is case where temp is unchanged, so heating and cooling same too */
		if( trace.lgTrConvg>=5 )
		{
			fprintf( ioQQQ, "       CoolIron NOT calling FeIILevelPops\n");
		}
	}

	/* evaluate some strong lines that would have been evaluated by the 16 level atom */
	FeIIFillLow16();

	/* now update total cooling and its deriv 
	 * all paths flow through here */
	/* FeII.Fecool = FeII.Fe2_large_cool; */
	CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_large_cool));

	/* add negative cooling to heating stack */
	thermal.heating[25][27] = MAX2(0.,FeII.Fe2_large_heat);

	/* counts as heating derivative if negative cooling */
	if( FeII.Fe2_large_cool > 0. )
	{
		/* >>chng 01 mar 16, add factor of 3 due to conv problems after changing damper */
		thermal.dCooldT += 3.*FeII.ddT_Fe2_large_cool;
	}

	if( trace.lgTrace && trace.lgCoolTr )
	{
		fprintf( ioQQQ, " Large FeII returns te, cooling, dc=%11.3e%11.3e%11.3e\n", 
		  phycon.te, FeII.Fe2_large_cool, FeII.ddT_Fe2_large_cool );
	}

	/* >>chng 05 nov 29, still do simple UV atom if only ground term is done */
	if( !FeII.lgFeIILargeOn )
	{

		/* following treatment of Fe II follows
		 * >>refer	fe2	model	Wills, B.J., Wills, D., Netzer, H. 1985, ApJ, 288, 143
		 * all elements are used, and must be set to zero if zero */

		/* set up space for simple  model of UV FeII emission */
		if( lgFirst )
		{
			/* will never do this again in this coreload */
			lgFirst = FALSE;
			/* allocate the 1D arrays*/
			if( (pops = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
				BadMalloc();
			if( (create = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
				BadMalloc();
			if( (destroy = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
				BadMalloc();
			if( (depart = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
				BadMalloc();
			/* create space for the 2D arrays */
			if( (AulPump = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
				BadMalloc();
			if( (CollRate = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
				BadMalloc();
			if( (AulDest = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
				BadMalloc();
			if( (AulEscp = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
				BadMalloc();
			if( (col_str = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
				BadMalloc();

			for( i=0; i<(NLFE2); ++i )
			{
				if( (AulPump[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
					BadMalloc();
				if( (CollRate[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
					BadMalloc();
				if( (AulDest[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
					BadMalloc();
				if( (AulEscp[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
					BadMalloc();
				if( (col_str[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
					BadMalloc();
			}
		}

		/*zero out all arrays, then check that upper diag remains zero below */
		for( i=0; i < NLFE2; i++ )
		{
			create[i] = 0.;
			destroy[i] = 0.;
			for( j=0; j < NLFE2; j++ )
			{
				/*data[j][i] = 0.;*/
				col_str[j][i] = 0.;
				AulEscp[j][i] = 0.;
				AulDest[j][i] = 0.;
				AulPump[j][i] = 0.;
			}
		}

		/* now put in real data for lines */
		AulEscp[1][0] = 1.;
		AulEscp[2][0] = ( TauLines[ipTuv3].Pesc + TauLines[ipTuv3].Pelec_esc)*TauLines[ipTuv3].Aul;
		AulDest[2][0] = TauLines[ipTuv3].Pdest*TauLines[ipTuv3].Aul;
		AulPump[0][2] = TauLines[ipTuv3].pump;
		/*ipdest[0][2] = TauLines[ipTuv3].ipCont;*/

		AulEscp[5][0] = (TauLines[ipTFe16].Pesc + TauLines[ipTFe16].Pelec_esc)*TauLines[ipTFe16].Aul;
		AulDest[5][0] = TauLines[ipTFe16].Pdest*TauLines[ipTFe16].Aul;
		/* continuum pumping of n=6 */
		AulPump[0][5] = TauLines[ipTFe16].pump;
		/* Ly-alpha pumping */
		if( dense.xIonDense[ipIRON][1] > 0. )
		{
			double
				PumpLyaFeII = iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]*dense.xIonDense[ipHYDROGEN][1]*
				EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul*
				hydro.dstfe2lya/dense.xIonDense[ipIRON][1];
			AulPump[0][5] += PumpLyaFeII;
		}

		AulEscp[2][1] = (TauLines[ipTr48].Pesc + TauLines[ipTr48].Pelec_esc)*TauLines[ipTr48].Aul;
		AulDest[2][1] = TauLines[ipTr48].Pdest*TauLines[ipTr48].Aul;
		AulPump[1][2] = TauLines[ipTr48].pump;

		AulEscp[5][1] = (TauLines[ipTFe26].Pesc + TauLines[ipTFe26].Pelec_esc)*TauLines[ipTFe26].Aul;
		AulDest[5][1] = TauLines[ipTFe26].Pdest*TauLines[ipTFe26].Aul;
		AulPump[1][5] = TauLines[ipTFe26].pump;

		AulEscp[3][2] = (TauLines[ipTFe34].Pesc + TauLines[ipTFe34].Pelec_esc)*TauLines[ipTFe34].Aul;
		AulDest[3][2] = TauLines[ipTFe34].Pdest*TauLines[ipTFe34].Aul;
		AulPump[2][3] = TauLines[ipTFe34].pump;

		AulEscp[4][2] = (TauLines[ipTFe35].Pesc + TauLines[ipTFe35].Pelec_esc)*TauLines[ipTFe35].Aul;
		AulDest[4][2] = TauLines[ipTFe35].Pdest*TauLines[ipTFe35].Aul;
		AulPump[2][4] = TauLines[ipTFe35].pump;

		AulEscp[5][3] = (TauLines[ipTFe46].Pesc + TauLines[ipTFe46].Pelec_esc)*TauLines[ipTFe46].Aul;
		AulDest[5][3] = TauLines[ipTFe46].Pdest*TauLines[ipTFe46].Aul;
		AulPump[3][5] = TauLines[ipTFe46].pump;

		AulEscp[5][4] = (TauLines[ipTFe56].Pesc + TauLines[ipTFe56].Pelec_esc)*TauLines[ipTFe56].Aul;
		AulDest[5][4] = TauLines[ipTFe56].Pdest*TauLines[ipTFe56].Aul;
		AulPump[4][5] = TauLines[ipTFe56].pump;

		/* these are collision strengths */
		col_str[1][0] = 1.;
		col_str[2][0] = 12.;
		col_str[3][0] = 1.;
		col_str[4][0] = 1.;
		col_str[5][0] = 12.;
		col_str[2][1] = 6.;
		col_str[3][1] = 1.;
		col_str[4][1] = 1.;
		col_str[5][1] = 12.;
		col_str[3][2] = 6.;
		col_str[4][2] = 12.;
		col_str[5][2] = 1.;
		col_str[4][3] = 1.;
		col_str[5][3] = 12.;
		col_str[5][4] = 6.;

		/*void atom_levelN(long,long,float,double[],double[],double[],double*,
		double*,double*,long*,float*,float*,STRING,int*);*/
		atom_levelN(NLFE2,
			dense.xIonDense[ipIRON][1],
			gFe2,
			ex,
			'K',
			pops,
			depart,
			&AulEscp ,
			&col_str,
			&AulDest,
			&AulPump,
			&CollRate,
			create,
			destroy,
			FALSE,/* say atom_levelN should evaluate coll rates from cs */
			/*&&ipdest,*/
			&FeII.Fe2_UVsimp_cool,
			&FeII.ddT_Fe2_UVsimp_cool,
			"FeII",
			&lgNegPop,
			&lgZeroPop,
			FALSE );

		/* lgNegPop positive if negative pops occured, negative if too cold */
		if( lgNegPop>0 /*negative if too cold - that is not negative and is OK */ )
		{
			fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population for simple UV FeII.\n");
		}

		/* add heating - cooling by this process */;
		CoolAdd("Fe 2",0,MAX2(0.,FeII.Fe2_UVsimp_cool));
		thermal.heating[25][27] = MAX2(0.,-FeII.Fe2_UVsimp_cool);
		thermal.dCooldT += FeII.ddT_Fe2_UVsimp_cool;

		/* LIMLEVELN is the dim of the PopLevels vector */
		ASSERT( NLFE2 <= LIMLEVELN );
		for( i=0; i<NLFE2; ++i)
		{
			atoms.PopLevels[i] = pops[i];
			atoms.DepLTELevels[i] = depart[i];
		}

		TauLines[ipTuv3].PopLo = pops[0];
		TauLines[ipTuv3].PopOpc = (pops[0] - pops[2]);
		TauLines[ipTuv3].phots = pops[2]*AulEscp[2][0];
		TauLines[ipTuv3].xIntensity = 
			TauLines[ipTuv3].phots*TauLines[ipTuv3].EnergyErg;

		TauLines[ipTr48].PopLo = pops[1];
		TauLines[ipTr48].PopOpc = (pops[1] - pops[2]);
		TauLines[ipTr48].phots = pops[2]*AulEscp[2][1];
		TauLines[ipTr48].xIntensity = 
			TauLines[ipTr48].phots*TauLines[ipTr48].EnergyErg;

		FeII.for7 = pops[1]*AulEscp[1][0]*4.65e-12;

		TauLines[ipTFe16].PopLo = pops[0];
		TauLines[ipTFe16].PopOpc = (pops[0] - pops[5]);
		TauLines[ipTFe16].phots = pops[5]*AulEscp[5][0];
		TauLines[ipTFe16].xIntensity = 
			TauLines[ipTFe16].phots*TauLines[ipTFe16].EnergyErg;

		TauLines[ipTFe26].PopLo = (pops[1]);
		TauLines[ipTFe26].PopOpc = (pops[1] - pops[5]);
		TauLines[ipTFe26].phots = pops[5]*AulEscp[5][1];
		TauLines[ipTFe26].xIntensity = 
			TauLines[ipTFe26].phots*TauLines[ipTFe26].EnergyErg;

		TauLines[ipTFe34].PopLo = pops[2];
		TauLines[ipTFe34].PopOpc = (pops[2] - pops[3]);
		TauLines[ipTFe34].phots = pops[3]*AulEscp[3][2];
		TauLines[ipTFe34].xIntensity = 
			TauLines[ipTFe34].phots*TauLines[ipTFe34].EnergyErg;

		TauLines[ipTFe35].PopLo = pops[2];
		TauLines[ipTFe35].PopOpc = (pops[2] - pops[4]);
		TauLines[ipTFe35].phots = pops[4]*AulEscp[4][2];
		TauLines[ipTFe35].xIntensity = 
			TauLines[ipTFe35].phots*TauLines[ipTFe35].EnergyErg;

		TauLines[ipTFe46].PopLo = pops[3];
		TauLines[ipTFe46].PopOpc = (pops[3] - pops[5]);
		TauLines[ipTFe46].phots = pops[5]*AulEscp[5][3];
		TauLines[ipTFe46].xIntensity = 
			TauLines[ipTFe46].phots*TauLines[ipTFe46].EnergyErg;

		TauLines[ipTFe56].PopLo = pops[4];
		/* tFe56(ipLnPopOpc) = pops(5) - pops(6)
		 * >>chng 96 apr 23 this levels mases due to Lya pump */
		TauLines[ipTFe56].PopOpc = (pops[4] - pops[5]*0.);
		TauLines[ipTFe56].phots = pops[5]*AulEscp[5][4];
		TauLines[ipTFe56].xIntensity = 
			TauLines[ipTFe56].phots*TauLines[ipTFe56].EnergyErg;

		/* Jack's funny FeII lines, data from 
		 * >>refer	fe2	energy	Johansson, S., Brage, T., Leckrone, D.S., Nave, G. &
		 * >>refercon Wahlgren, G.M. 1995, ApJ 446, 361 */
		PutCS(10.,&TauLines[ipT191]);
		atom_level2(&TauLines[ipT191]);
	}

	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && iteration>1 && nzone>=4/**/ )
		{
			fprintf(ioQQQ,"DEBUG2\t%.2e\t%.2e\t%.2e\n",
				phycon.te,
				FeII.Fe2_large_cool , 
				FeII.Fe2_UVsimp_cool );
		}
	}

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

	return;

}

void CoolIron()
{
	long int i;
	/*static long int **ipdest;*/
	double cs ,
	  cs12, cs13, cs23,
	  cs2s2p, 
	  cs2s3p;
	float p2, 
		rate;

	static int lgFe22First=TRUE;
	static long int *ipFe22Pump=NULL,
		nFe22Pump=0;
	double Fe22_pump_rate;

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

	/* cooling by FeI 24m, 34.2 m */
	/* >>chng 03 nov 15, add these lines */
	/*TODO	2	- ground term is actually a fix level system, the vectors are
	 * created, with pointers ipFe1_54m , ipFe1_111m, must add collision date, use
	 * larger model atom */
	/*>>refer	Fe1	cs	Hollenbach & McKee 89 */
	/* the 24.0 micron line */
	rate = (float)(1.2e-7 * dense.eden + 
		/* >>chng 05 jul 05, eden to cdsqte */
		/*8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
		8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
	LineConvRate2CS( &TauLines[ipFe1_24m]  , rate );

	rate = (float)(9.3e-8 * dense.eden + 
		/* >>chng 05 jul 05, eden to cdsqte */
		/*5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
		5.3e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
	LineConvRate2CS( &TauLines[ipFe1_35m]  , rate );

	rate = (float)(1.2e-7 * dense.eden + 
		/* >>chng 05 jul 05, eden to cdsqte */
		/*6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);*/
		6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]);
	TauDummy.gHi = TauLines[ipFe1_35m].gHi;
	LineConvRate2CS( &TauDummy  , rate );
	/* this says that line is a dummy, not real one */
	TauDummy.gHi = 0.;

	atom_level3(&TauLines[ipFe1_24m],&TauLines[ipFe1_35m],&TauDummy);

	/* series of FeI lines from Dima Verner's list, each 2-lev atom
	 *
	 * Fe I 3884 */
	MakeCS(&TauLines[ipFeI3884]);
	atom_level2(&TauLines[ipFeI3884]);

	/* Fe I 3729 */
	MakeCS(&TauLines[ipFeI3729]);
	atom_level2(&TauLines[ipFeI3729]);

	/* Fe I 3457 */
	MakeCS(&TauLines[ipFeI3457]);
	atom_level2(&TauLines[ipFeI3457]);

	/* Fe I 3021 */
	MakeCS(&TauLines[ipFeI3021]);
	atom_level2(&TauLines[ipFeI3021]);

	/* Fe I 2966 */
	MakeCS(&TauLines[ipFeI2966]);
	atom_level2(&TauLines[ipFeI2966]);

	/* >>chng 05 dec 03, move Fe2 FeII Fe II cooling into separate routine */
	Fe2_cooling();

	/* lump 3p and 3f together; cs=
	 * >>refer	fe3	as	Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384
	 * A from
	 * >>refer	fe3	as	Garstang, R.H., 1957, Vistas in Astronomy, 1, 268
	 * FE III 5270, is 20.9% of total 
	 * >>chng 05 feb 18, Kevin Blagrave email
	 * average wavelength is 4823 with statistical weight averaging of upper energy level,
	 * as per , change 5th number from 2.948 to 2.984, also photon energy
	 * from 3.78 to 4.12 */

	/* >>chng 05 dec 16, FeIII update by Kevin Blagrave */
	/*CoolHeavy.c5270 = atom_pop2(5.53,25.,30.,0.435,2.984e4,dense.xIonDense[ipIRON][2])*
	  4.12e-12;
	CoolAdd("Fe 3",5270,CoolHeavy.c5270);*/
	/* FeIII 1122 entire multiplet - atomic data=A's Dima, CS = guess */
	PutCS(25.,&TauLines[ipT1122]);
	atom_level2(&TauLines[ipT1122]);

	/* call 14 level atom */
	Fe3Lev14();

	/* call 12 level atom */
	Fe4Lev12();

	/* FE V 3892 + 3839, data from Shields */
	CoolHeavy.c3892 = atom_pop2(7.4,25.,5.,0.6,3.7e4,dense.xIonDense[ipIRON][4])*
	  5.11e-12;
	CoolAdd("Fe 5",3892,CoolHeavy.c3892);

	/* FE VI 5177+5146+4972+4967
	 * data from 
	 * >>refer	fe6	as	Garstang, R.H., Robb, W.D., Rountree, S.P. 1978, ApJ, 222, 384 */
	CoolHeavy.c5177 = atom_pop2(1.9,28.,18.,0.52,2.78e4,dense.xIonDense[ipIRON][5])*
	  3.84e-12;
	CoolAdd("Fe 6",5177,CoolHeavy.c5177);

	/* >>chng 04 nov 04, add multi-level fe7 */
	Fe7Lev8();

	/* Fe IX 245,242
	 * all atomic data from 
	 * >>refer	fe9	all	Flower, D.R. 1976, A&A, 56, 451 */
	/* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */
	PutCS(0.123,&TauLines[ipT245]);
	/* AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41)
	 * C245 = AtomSeqBeryllium( .087 ,.038 ,.188 , t245 ,71. )  * 8.14E-11 */
	AtomSeqBeryllium(.087,.038,.188,&TauLines[ipT245],71.);

	CoolHeavy.c242 = atoms.PopLevels[3]*8.22e-11*71.;

	/* Fe X Fe 10 Fe10 6374, A from 
	 * >>referold	fe10	as	Mason, H. 1975, MNRAS 170, 651
	 * >>referold	fe10	cs	Mohan, M., Hibbert, A., & Kingston, A.E. 1994, ApJ, 434, 389
	 * >>referold	fe10	as	Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209
	 * does not agree with Mohan et al. by factor of 4
	 * 20x larger than Mason numbers used pre 1994 */
	/*cs = 13.3-2.*MIN2(7.0,phycon.alogte); */
	/*cs = MAX2(0.1,cs); */
	/* TODO	2	following to stop insane FeX strengths
	 * >>chng 96 jul 11, removed 1 / 10 factor, so using real data, 90.01
	 * cs = cs * 0.1
	 * >>chng 96 jun 03, transferred following
	 * >>chng 97 jan 31, I give up on this mess, use cs of unity */
	/*cs = 1.;*/
	/* >>chng 00 dec 05, update Fe10 collisions strength to Tayal 2000 */
	/* >>refer	fe10	cs	Tayal, S.S., 2000, ApJ, 544, 575-580 */
	/* >>chng 04 nov 10, 172.9 was mult rather than div by temp law,
	 * had no effect due to min on cs that lie below it */
	/*cs = 172.9 /( phycon.te30 * phycon.te05 * phycon.te02 *  phycon.te005 );*/
	/* tabulated cs ends at 30,000K, do not allow cs to grow larger than largest
	 * tabulated value */
	/*cs = MIN2( cs, 3.251 );*/
	/* >>chng 05 dec 19, update As to 
	 * >>refer	fe10	As	Aggarwal, K.M. & Keenan, F.P. 2004, A&A, 427, 763 
	 * value changed from 54.9 to 68.9 */
	/* >>chng 05 dec 19, update cs to
	 *>>refer	fe10	cs	Aggarwal, K.M. & Keenan, F.P. 2005, A&A, 439, 1215 */
	cs = Fe_10_11_13_cs(
	/* ion, one of 10, 11, or 13 on physics scale */
		10, 
	/* inital and final index of levels, with lowest energy 1, highest of 5 
	 * on physics scale */
		1, 
		2 );

	PutCS(cs,&TauLines[ipFe106375]);
	atom_level2(&TauLines[ipFe106375]);
	/* c6374 = atom_pop2(cs ,4.,2.,69.5,2.26e4,dense.xIonDense(26,10))*3.12e-12
	 * dCooldT = dCooldT + c6374 * 2.26e4 * tsq1
	 * call CoolAdd( 'Fe10' , 6374 , c6374 )
	 *
	 * this is the E1 line that can pump [FeX] */
	cs = 0.85*sexp(0.045*1.259e6/phycon.te);
	cs = MAX2(0.05,cs);
	PutCS(cs,&TauLines[ipT352]);
	atom_level2(&TauLines[ipT352]);

	/* Fe XI fe11 fe 11, 7892, 6.08 mic, CS from 
	 * >>refer	fe11	as	Mason, H. 1975, MNRAS 170, 651
	 * A from 
	 * >>refer	fe11	as	Mendoza, C., & Zeippen, C.J. 1983, MNRAS, 202, 981 
	 * following reference has very extensive energy levels and As
	 * >>refer	fe11	as	Fritzsche, S., Dong, C.Z., & Traebert, E., 2000, MNRAS, 318, 263 */
	/*cs = 0.27;*/
	/* >>chng 97 jan 31, set cs to unity, above ignored resonances */
	/*cs = 1.;*/
	/* >>chng 00 dec 05, update Fe11 CS to Tayal 2000 */
	/* >>referold	fe11	cs	Tayal, S.S., 2000, ApJ, 544, 575-580 */
	/* this is the low to mid transition */
	/*cs = 0.0282 * phycon.te30*phycon.te05*phycon.te01;*/
	/* CS is about 2.0 across broad temp range in following reference:
	 * >>refer	Fe11	cs	Aggarwal, K.M., & Keenan, F.P. 2003, MNRAS, 338, 412 */
	/*cs = 2.;
	PutCS(cs,&TauLines[ipTFe07]);*/
	/* Tayal value is 10x larger than previous values */
	/* Aggarwal & Keenan is about same as Tayal */
	/*cs = 0.48; */
	/*cs = 0.50;
	PutCS(cs,&TauLines[ipTFe61]);*/
	/* Tayal value is 3x larger than previous values */
	/*cs = 0.35; tayal number */
	/*cs = 0.55;
	PutCS(cs,&TauDummy);
	atom_level3(&TauLines[ipTFe07],&TauLines[ipTFe61],&TauDummy);*/

	/* >>refer	fe11	cs	Kafatos, M., & Lynch, J.P. 1980, ApJS, 42, 611 */
	/*CoolHeavy.c1467 = atom_pop3(9.,5.,1.,.24,.028,.342,100.,1012.,9.3,5.43e4,
	  6.19e4,&p2,dense.xIonDense[ipIRON][11-1],0.,0.,0.)*1012.*1.36e-11;
	CoolHeavy.c2649 = p2*91.0*7.512e-12;*/
	/*CoolAdd("Fe11",1467,CoolHeavy.c1467);
	CoolAdd("Fe11",2469,CoolHeavy.c2649);*/

	/* >>chng 05 dec 18, use give level Fe 11 model */
	Fe11Lev5();

	/* A's for 2-1 and 3-1 from 
	 * >>refer	fe12	as	Tayal, S.S., & Henry, R.J.W. 1986, ApJ, 302, 200
	 * CS fro 
	 * >>refer	fe12	cs	Tayal, S.S., Henry, R.J.W., Pradhan, A.K. 1987, ApJ, 319, 951
	 * a(3-2) scaled from both ref */
	CoolHeavy.c2568 = atom_pop3(4.,10.,6.,0.72,0.69,2.18,8.1e4,1.84e6,1.33e6,
	  6.37e4,4.91e4,&p2,dense.xIonDense[ipIRON][12-1],0.,0.,0.)*1.33e6*6.79e-12;
	CoolAdd("Fe12",2568,CoolHeavy.c2568);
	CoolHeavy.c1242 = CoolHeavy.c2568*2.30*1.38;
	CoolAdd("Fe12",1242,CoolHeavy.c1242);
	CoolHeavy.c2170 = p2*8.09e4*8.82e-12;
	CoolAdd("Fe12",2170,CoolHeavy.c2170);

	/* [Fe XIII] fe13 fe 13 1.07, 1.08 microns */
	/* >>chng 00 dec 05, update Fe13 CS to Tayal 2000 */
	/* >>refer	fe13	cs	Tayal, S.S., 2000, ApJ, 544, 575-580
	cs = 5.08e-3 * phycon.te30* phycon.te10;
	PutCS(cs,&TauLines[ipFe1310]);
	PutCS(2.1,&TauLines[ipFe1311]);
	PutCS(.46,&TauDummy);
	atom_level3(&TauLines[ipFe1310],&TauLines[ipFe1311],&TauDummy); */

	/* Fe13 lines from 1D and 1S -
	 >>chng 01 aug 07, added these lines */
	/* >>refer	fe11	cs	Tayal, S.S., 2000, ApJ, 544, 575-580 */
	/* >>refer	fe13	as	Shirai, T., Sugar, J., Musgrove, A., & Wiese, W.L., 2000., 
	 * >>refercon	J Phys Chem Ref Data Monograph 8 */
	/* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) 
	CoolHeavy.fe13_1216 = atom_pop3(  9.,5.,1.,  5.08,0.447,0.678,  103.2 ,1010.,7.99,
	  48068.,43440.,&p2,dense.xIonDense[ipIRON][13-1],0.,0.,0.)*1010.*
	  1.64e-11;
	CoolHeavy.fe13_2302 = CoolHeavy.fe13_1216*0.528*0.00791;
	CoolHeavy.fe13_3000 = p2*103.2*7.72e-12;*/
	/*CoolAdd("Fe13",1216,CoolHeavy.fe13_1216);
	CoolAdd("Fe13",3000,CoolHeavy.fe13_3000);
	CoolAdd("Fe13",2302,CoolHeavy.fe13_2302);*/

	/* >>chng 05 dec 18, use give level Fe 13 model */
	Fe13Lev5();

	/* Fe XIV 5303, cs from 
	 * >>refer	Fe14	cs	Storey, P.J., Mason, H.E., Saraph, H.E., 1996, A&A, 309, 677 */
	fe14cs(phycon.alogte,&cs);
	/* >>chng 97 jan 31, set cs to unity, above is VERY large, >10 */
	/* >>chng 01 may 30, back to their value, as per discussion at Lexington 2000 */
	/* >>chng 05 aug 04, following line commented out, had set 5303 to constant value */
	/*cs = 3.1;*/
	/* >>chng 05 dec 22, A from 60.3 to 59.7, experimental value 59.7 +/- 0.4 in
	 * >>refer	Fe14	as	Trabert, E. 2004, A&A, 415, L39 */
	CoolHeavy.c5303 = atom_pop2(cs,2.,4.,59.7,2.71e4,dense.xIonDense[ipIRON][14-1])*
	  3.75e-12;
	thermal.dCooldT += CoolHeavy.c5303*2.71e4*thermal.tsq1;
	CoolAdd("Fe14",5303,CoolHeavy.c5303);

	/* Fe 18 974.8A;, cs from 
	 * >>refer	fe18	cs	Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */
	cs = MIN2(0.143,0.804/(phycon.te10*phycon.te05));
	PutCS(cs,&TauLines[ipFe18975]);
	atom_level2(&TauLines[ipFe18975]);

	/* O-like Fe19, 3P ground term, 7046.72A vac wl, 1328.90A */
	/* >>refer	fe19	cs	Butler, K., & Zeippen, C.J., 2001, A&A, 372, 1083 */
	cs12 = 0.0627 / phycon.te03;
	cs13 = 0.692 /(phycon.te10*phycon.te01);
	cs23 = 0.04;
	/*CoolHeavy.c7082 = atom_pop3(5.,1.,3.,0.0132,0.0404,0.0137,0.505,1.46e4,*/
	/* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
	CoolHeavy.c7082 = atom_pop3(5.,1.,3.,cs12,cs13,cs23,0.505,1.46e4,
	  41.2,1.083e5,2.03e4,&p2,dense.xIonDense[ipIRON][19-1],0.,0.,0.)*41.2*
	  2.81e-12;
	CoolHeavy.c1118 = CoolHeavy.c7082*354.4*6.335;
	CoolHeavy.c1328 = p2*0.505*1.50e-11;
	CoolAdd("Fe19",7082,CoolHeavy.c7082);
	CoolAdd("Fe19",1118,CoolHeavy.c1118);
	CoolAdd("Fe19",1328,CoolHeavy.c1328);

	CoolHeavy.c592 = atom_pop2(0.0913,9.,5.,1.64e4,2.428e5,dense.xIonDense[ipIRON][19-1])*
	  3.36e-11;
	CoolAdd("Fe19",592,CoolHeavy.c592);

	/* >>chng 01 aug 10, add this line */
	/* Fe20 721.40A, 578*/
	/* >>refer	fe20	cs	Butler, K., & Zeippen, C.J., 2001, A&A, 372, 1078 */
	/*>>refer	fe20	as	Merkelis, G., Martinson, I., Kisielius, R., & Vilkas, M.J., 1999, 
	  >>refercon	Physica Scripta, 59, 122 */
	cs = 1.17 /(phycon.te20/phycon.te01);
	PutCS(cs , &TauLines[ipTFe20_721]);
	cs = 0.248 /(phycon.te10/phycon.te01);
	PutCS(cs , &TauLines[ipTFe20_578]);
	cs = 0.301 /(phycon.te10/phycon.te02);
	PutCS(cs , &TauDummy);
	atom_level3(&TauLines[ipTFe20_721],&TauLines[ipTFe20_578],&TauDummy);

	/* >>chng 00 oct 26, much larger CS for following transition */
	/* Fe 21 1354, 2299 A, cs eval at 1e6K
	 * >>refer	Fe21	cs	Aggarwall, K.M., 1991, ApJS 77, 677 */
	PutCS(0.072,&TauLines[ipTFe13]);
	PutCS(0.269,&TauLines[ipTFe23]);
	PutCS(0.055,&TauDummy);
	atom_level3(&TauLines[ipTFe13],&TauLines[ipTFe23],&TauDummy);

	/*>>refer	Fe22	energy	Feldman, U., Curdt, W., Landi, E., Wilhelm, K., 2000, ApJ, 544, 508 */
	/*>>chng 00 oct 26, added these fe 21 lines, removed old CoolHeavy.fs846, the lowest */

	/* one time initialization if first call, and level 2 lines are on */
	if( lgFe22First && nWindLine )
	{
		lgFe22First = FALSE;
		nFe22Pump = 0;
		for( i=0; i<nWindLine; ++i )
		{
			/* don't test on nelem==ipIRON since lines on physics, not C, scale */
			if( TauLine2[i].nelem ==26 && TauLine2[i].IonStg==22  )
			{
				++nFe22Pump;
			}
		}
		if( nFe22Pump<0 )
			TotalInsanity();
		else if( nFe22Pump > 0 )
			/* create the space - can't malloc 0 bytes */
			if( (ipFe22Pump = (long *)MALLOC((unsigned)(nFe22Pump)*sizeof(long) ) )==NULL )
				BadMalloc();
		nFe22Pump = 0;
		for( i=0; i<nWindLine; ++i )
		{
			/* don't test on nelem==ipIRON since lines on physics, not C, scale */
			if( TauLine2[i].nelem ==26 && TauLine2[i].IonStg==22  )
			{
#				define	PRT_DEBUG	FALSE
#				if	PRT_DEBUG
				DumpLine( &TauLine2[i] );
#				endif
				ipFe22Pump[nFe22Pump] = i;
				++nFe22Pump;
			}
		}
	}
	else
		/* level 2 lines are not enabled */
		nFe22Pump = 0;

	/* now sum pump rates */
	Fe22_pump_rate = 0.;
	for( i=0; i<nFe22Pump; ++i )
	{
		Fe22_pump_rate += TauLine2[ipFe22Pump[i]].pump;
#		if	PRT_DEBUG
		fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n",
			i,
			TauLine2[ipFe22Pump[i]].WLAng , TauLine2[ipFe22Pump[i]].pump );
#		endif
#		undef	PRT_DEBUG
	}

	/*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
	/*>>refer	fe22	cs	Zhang, H.L., Graziani, M., & Pradhan, A.K., 1994, A&A 283, 319*/
	/*>>refer	fe22	as	Dankwort, W., & Trefftz, E., 1978, A&A 65, 93-98 */
	AtomSeqBoron(&TauLines[ipFe22_846], 
	  &TauLines[ipFe22_247], 
	  &TauLines[ipFe22_217], 
	  &TauLines[ipFe22_348], 
	  &TauLines[ipFe22_292], 
	  &TauLines[ipFe22_253], 
	  /*double cs40,
		double cs32,
		double cs42,
		double cs43, */
	  0.00670 , 0.0614 , 0.0438 , 0.122 , Fe22_pump_rate ,"Fe22");

	/* Fe 22 845.6, C.S.=guess, A from
	 * >>refer	fe22	as	Froese Fischer, C. 1983, J.Phys. B, 16, 157 
	CoolHeavy.fs846 = atom_pop2(0.2,2.,4.,1.5e4,1.701e5,dense.xIonDense[ipIRON][22-1])*
	  2.35e-11;
	CoolAdd("Fe22",846,CoolHeavy.fs846);*/

	/*TODO	2	update atomic data to Chidichimo et al 1999 AASup 137, 175*/
	/* FE 23 262.6, 287A, 1909-LIKE, 
	 * cs from 
	 * >>refer	fe23	cs	Bhatia, A.K., & Mason, H.E. 1986, A&A, 155, 413 */
	CoolHeavy.c263 = atom_pop2(0.04,1.,9.,1.6e7,5.484e5,dense.xIonDense[ipIRON][23-1])*
	  7.58e-11;
	CoolAdd("Fe23",262,CoolHeavy.c263);

	/* FE 24, 255, 192 Li seq doublet
	 * >>refer	fe24	cs	Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
	ligbar(26,&TauLines[ipT192],&TauLines[ipT11],&cs2s2p,&cs2s3p);

	PutCS(cs2s2p,&TauLines[ipT192]);
	atom_level2(&TauLines[ipT192]);

	/* funny factor (should have been 0.5) due to energy change */
	PutCS(cs2s2p*0.376,&TauLines[ipT255]);
	atom_level2(&TauLines[ipT255]);

	PutCS(cs2s3p,&TauLines[ipT11]);
	atom_level2(&TauLines[ipT11]);

	/* TODO	2		 * following not in cooling function */
	TauLines[ipT353].PopOpc = dense.xIonDense[ipIRON][11-1];
	TauLines[ipT353].PopLo = dense.xIonDense[ipIRON][11-1];
	TauLines[ipT353].PopHi = 0.;
	TauLines[ipT347].PopOpc = dense.xIonDense[ipIRON][14-1];
	TauLines[ipT347].PopLo = dense.xIonDense[ipIRON][14-1];
	TauLines[ipT347].PopHi = 0.;


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

/*fe14cs compute collision strength for forbidden transition 
 * w/in Fe XIV ground term. From 
 * >>refer	fe14	cs	Storey, P.J., Mason, H.E., Saraph, H.E., 1996, A&A, 309, 677
 * */
static void fe14cs(double te1, 
	  double *csfe14)
{
	double a, 
	  b, 
	  c, 
	  d, 
	  telog1, 
	  telog2, 
	  telog3;

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

	/* limit range in log T: */
	telog1 = te1;
	telog1 = MIN2(9.0,telog1);
	telog1 = MAX2(4.0,telog1);

	/* compute square and cube */
	telog2 = telog1*telog1;
	telog3 = telog2*telog1;

	/* compute cs:
	 * */
	if( telog1 <= 5.0 )
	{
		a = 557.05536;
		b = -324.56109;
		c = 63.437974;
		d = -4.1365147;
		*csfe14 = a + b*telog1 + c*telog2 + d*telog3;
	}
	else
	{
		a = 0.19515493;
		b = 2.9404407;
		c = 4.9578944;
		d = 0.79887506;
		*csfe14 = a + b*exp(-0.5*((telog1-c)*(telog1-c)/d));
	}

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

/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
static void Fe4Lev12(void)
{
#	define	NLFE4	12
	int lgNegPop,
		lgZeroPop;
	long int i, 
	  j, 
	  _r;
	static int lgFirst=TRUE;

	double dfe4dt;

	/*static long int **ipdest; */
	static double 
		**AulEscp ,
		**col_str ,
		**AulDest, 
		depart[NLFE4],
		pops[NLFE4], 
		destroy[NLFE4], 
		create[NLFE4] ,
		**CollRate,
		**AulPump;

	static double Fe4A[NLFE4][NLFE4], 
	  Fe4CS[NLFE4][NLFE4];

	static double gfe4[NLFE4]={6.,12.,10.,6.,8.,6.,4.,2.,8.,2.,6.,4.};

	/* excitation energies in Kelvin 
	static double ex[NLFE4]={0.,46395.8,46464.,46475.6,46483.,50725.,
	  50838.,50945.,55796.,55966.,56021.,56025.};*/
	/*>>refer	Fe3	energies	version 3 NIST Atomic Spectra Database */
	/*>>chng 05 dec 17, from Kelvin above to excitation energies in wn */
	static double excit_wn[NLFE4]={0.,32245.5,32292.8,32301.2,32305.7,35253.8,
	  35333.3,35406.6,38779.4,38896.7,38935.1,38938.2};

	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {0.,0.,0.,1.e-5,0.,1.368,.89,0.,
		  1.3e-3,1.8e-4,.056,.028};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[0][j-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {0.,0.,2.6e-8,0.,0.,0.,0.,0.,1.7e-7,
		  0.,0.,0.};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[1][j-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {0.,0.,0.,0.,3.5e-7,6.4e-10,0.,0.,
		  6.315e-4,0.,6.7e-7,0.};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[2][j-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {0.,0.,0.,0.,1.1e-6,6.8e-5,8.6e-6,
		  3.4e-10,7.6e-5,1.e-7,5.8e-4,2.8e-4};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[3][j-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {0.,0.,0.,0.,0.,1.5e-5,1.3e-9,0.,
		  7.6e-4,0.,1.1e-6,6.0e-7};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[4][j-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.,0.,0.,0.,0.,0.,1.1e-5,1.2e-13,
		  .038,9.9e-7,.022,.018};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[5][j-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {0.,0.,0.,0.,0.,0.,0.,3.7e-5,2.9e-6,
		  .034,3.5e-3,.039};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[6][j-1] = _itmp6[_r++];
			}
		}
		{ static double _itmp7[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,.058,
		  3.1e-6,1.4e-3};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[7][j-1] = _itmp7[_r++];
			}
		}
		{ static double _itmp8[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-4,
		  3.1e-14};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[8][j-1] = _itmp8[_r++];
			}
		}
		{ static double _itmp9[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.9e-19,
		  1.0e-5};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[9][j-1] = _itmp9[_r++];
			}
		}
		{ static double _itmp10[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		  0.,1.3e-7};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[10][j-1] = _itmp10[_r++];
			}
		}
		for( j=1; j <= NLFE4; j++ )
		{
			Fe4A[11][j-1] = 0.;
			}
		for( i=1; i <= NLFE4; i++ )
		{
			Fe4CS[0][i-1] = 0.;
			}
		{ static double _itmp11[] = {0.98,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		  0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[1][i-1] = _itmp11[_r++];
			}
		}
		{ static double _itmp12[] = {0.8167,3.72,0.,0.,0.,0.,0.,0.,
		  0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[2][i-1] = _itmp12[_r++];
			}
		}
		{ static double _itmp13[] = {0.49,0.0475,0.330,0.,0.,0.,0.,
		  0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[3][i-1] = _itmp13[_r++];
			}
		}
		{ static double _itmp14[] = {0.6533,0.473,2.26,1.64,0.,0.,0.,
		  0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[4][i-1] = _itmp14[_r++];
			}
		}
		{ static double _itmp15[] = {0.45,0.686,0.446,0.106,0.254,0.,
		  0.,0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[5][i-1] = _itmp15[_r++];
			}
		}
		{ static double _itmp16[] = {0.30,0.392,0.152,0.269,0.199,0.605,
		  0.,0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[6][i-1] = _itmp16[_r++];
			}
		}
		{ static double _itmp17[] = {0.15,0.0207,0.190,0.0857,0.166,
		  0.195,0.327,0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[7][i-1] = _itmp17[_r++];
			}
		}
		{ static double _itmp18[] = {0.512,1.23,0.733,0.174,0.398,0.623,
		  0.335,0.102,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[8][i-1] = _itmp18[_r++];
			}
		}
		{ static double _itmp19[] = {0.128,0.0583,0.185,0.200,0.188,
		  0.0835,0.127,0.0498,0.0787,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[9][i-1] = _itmp19[_r++];
			}
		}
		{ static double _itmp20[] = {0.384,0.578,0.534,0.363,0.417,
		  0.396,0.210,0.171,0.810,0.101,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[10][i-1] = _itmp20[_r++];
			}
		}
		{ static double _itmp21[] = {0.256,0.234,0.306,0.318,0.403,
		  0.209,0.195,0.112,0.195,0.458,0.727,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[11][i-1] = _itmp21[_r++];
			}
		}
		_aini = 0;
	}

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

	if( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* allocate the 1D arrays*/
		/* create space for the 2D arrays */
		if( (AulPump = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulDest = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulEscp = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (col_str = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE4); ++i )
		{
			if( (AulPump[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulDest[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulEscp[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (col_str[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
		}
	}

	/* bail if no Fe+3 */
	if( dense.xIonDense[ipIRON][3] <= 0. )
	{
		fe.Fe4CoolTot = 0.;
		fe.fe40401 = 0.;
		fe.fe42836 = 0.;
		fe.fe42829 = 0.;
		fe.fe42567 = 0.;
		fe.fe41207 = 0.;
		fe.fe41206 = 0.;
		fe.fe41106 = 0.;
		fe.fe41007 = 0.;
		fe.fe41008 = 0.;
		fe.fe40906 = 0.;
		CoolAdd("Fe 4",0,0.);

		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( NLFE4 <= LIMLEVELN);
		for( i=0; i < NLFE4; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->Fe4Lev12()\n", debug_fp );
#		endif
		return;
	}
	/* number of levels in model ion */

	/* these are in wavenumbers
	 * data excit_wn/ 0., 32245.5, 32293., 32301.2, 
	 *  1  32306., 35254., 35333., 35407., 38779., 38897., 38935.,
	 *  2  38938./ 
	 * excitation energies in Kelvin */

	/* A's are from Garstang, R.H., MNRAS 118, 572 (1958).
	 * each set is for a lower level indicated by second element in array,
	 * index runs over upper level
	 * A's are saved into arrays as data(up,lo) */

	/* collision strengths from Berrington and Pelan  Ast Ap S 114, 367.
	 * order is cs(low,up) */

	/* all elements are used, and must be set to zero if zero */
	for( i=0; i < NLFE4; i++ )
	{
		create[i] = 0.;
		destroy[i] = 0.;
		for( j=0; j < NLFE4; j++ )
		{
			/*data[j][i] = 1e33;*/
			col_str[j][i] = 0.;
			AulEscp[j][i] = 0.;
			AulDest[j][i] = 0.;
			AulPump[j][i] = 0.;
		}
	}

	/* fill in einstein As and collision strengths */
	for( i=0; i < NLFE4; i++ )
	{
		for( j=i + 1; j < NLFE4; j++ )
		{
			/*data[i][j] = Fe4A[i][j];*/
			AulEscp[j][i] = Fe4A[i][j];
			/*data[j][i] = Fe4CS[j][i];*/
			col_str[j][i] = Fe4CS[j][i];
		}
	}

	/* leveln itself is well-protected against zero abundances,
	 * low temperatures */

	atom_levelN(NLFE4,
		dense.xIonDense[ipIRON][3],
		gfe4,
		excit_wn,
		'w',
		pops,
		depart,
		&AulEscp ,
		&col_str ,
		&AulDest,
		&AulPump,
		&CollRate,
		create,
		destroy,
		/* say atom_levelN should evaluate coll rates from cs */
		FALSE,
		&fe.Fe4CoolTot,
		&dfe4dt,
		"FeIV",
		/* lgNegPop positive if negative pops occured, negative if too cold */
		&lgNegPop,
		&lgZeroPop,
		FALSE );

	/* LIMLEVELN is the dim of the PopLevels vector */
	ASSERT( NLFE4 <= LIMLEVELN );
	for( i=0; i<NLFE4; ++i)
	{
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}

	if( lgNegPop>0 )
	{
		fprintf( ioQQQ, " fe4levl2 found negative populations\n" );
		ShowMe();
		puts( "[Stop in fe4levl2]" );
		cdEXIT(EXIT_FAILURE);
	}

	CoolAdd("Fe 4",0,fe.Fe4CoolTot);

	thermal.dCooldT += dfe4dt;

	/* three transitions hst can observe
	 * 4 -1 3095.9A and 5 -1 3095.9A */
	fe.fe40401 = (pops[3]*Fe4A[0][3]*(excit_wn[3] - excit_wn[0]) + 
		pops[4]*Fe4A[0][4]*(excit_wn[4] - excit_wn[0]) )*T1CM*BOLTZMANN;

	fe.fe42836 = pops[5]*Fe4A[0][5]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN;

	fe.fe42829 = pops[6]*Fe4A[0][6]*(excit_wn[5] - excit_wn[0])*T1CM*BOLTZMANN;

	fe.fe42567 = (pops[10]*Fe4A[0][10]*(excit_wn[10] - excit_wn[0]) + 
		pops[11]*Fe4A[0][11]*(excit_wn[10] - excit_wn[0]))*T1CM*BOLTZMANN;

	fe.fe41207 = pops[11]*Fe4A[6][11]*(excit_wn[11] - excit_wn[6])*T1CM*BOLTZMANN;
	fe.fe41206 = pops[11]*Fe4A[5][11]*(excit_wn[11] - excit_wn[5])*T1CM*BOLTZMANN;
	fe.fe41106 = pops[10]*Fe4A[5][10]*(excit_wn[10] - excit_wn[5])*T1CM*BOLTZMANN;
	fe.fe41007 = pops[9]*Fe4A[6][9]*(excit_wn[9] - excit_wn[6])*T1CM*BOLTZMANN;
	fe.fe41008 = pops[9]*Fe4A[7][9]*(excit_wn[9] - excit_wn[7])*T1CM*BOLTZMANN;
	fe.fe40906 = pops[8]*Fe4A[5][8]*(excit_wn[8] - excit_wn[5])*T1CM*BOLTZMANN;

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

#if 0
#ifdef NLFE2
#	undef NLFE2
#endif
#define	NLFE2	16

/*Fe2Lev16 compute populations and cooling due to 16 level Fe II ion */
static void Fe2Lev16(void)
{
	int lgNegPop,
		lgZeroPop;
	long int i, 
	  j, 
	  nUsed, 
	  _r;
	/*static long int **ipdest; */

	static double 
		**AulEscp ,
		**col_str ,
		**AulDest, 
		*pops, 
		*create, 
		*destroy, 
		*depart,
		**CollRate,
		**AulPump;
	static double Fe2A[NLFE2][NLFE2], 
	  Fe2CS[NLFE2][NLFE2];

	static int lgFirst=TRUE;

	static double gfe2[NLFE2]={10.,8.,6.,4.,2.,10.,8.,6.,4.,8.,6.,4.,
	  2.,6.,4.,2.};
	/* excitation energies in Kelvin */
	static double ex[NLFE2]={0.,553.646,960.680,1241.151,1405.810,2694.300,
	  3496.489,4083.318,4485.486,11446.297,12074.546,12489.670,12728.966,
	  19387.344,19673.344,20006.633};
	static int _aini = 1;

	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {0.00e00,2.13e-03,1.00e-16,1.00e-16,
		  1.00e-16,2.80e-04,1.00e-04,1.00e-16,1.00e-16,4.74e-03,1.00e-16,
		  1.00e-16,1.00e-16,1.00e-16,1.00e-16,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[0][j-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {0.00e00,0.00e00,1.57e-03,1.00e-16,
		  1.00e-16,1.00e-16,1.80e-04,6.90e-04,1.00e-16,1.31e-03,4.80e-04,
		  1.00e-16,1.00e-16,6.64e-03,1.00e-16,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[1][j-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {0.00e00,0.00e00,0.00e00,7.18e-04,
		  1.00e-16,1.00e-16,3.40e-05,1.00e-04,2.80e-05,9.70e-04,1.98e-03,
		  1.00e-16,1.00e-16,1.00e-16,6.81e-03,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[2][j-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  1.88e-04,1.00e-16,1.00e-16,2.50e-05,4.00e-05,1.00e-16,1.17e-03,
		  2.45e-03,7.80e-04,1.00e-16,1.00e-16,6.23e-03};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[3][j-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,1.00e-16,1.00e-16,1.00e-16,9.90e-06,1.00e-16,1.00e-16,
		  1.08e-03,3.32e-03,1.00e-16,2.20e-04,1.93e-03};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[4][j-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,5.84e-03,1.00e-16,1.00e-16,5.98e-03,3.12e-03,
		  1.00e-16,1.00e-16,3.56e-02,1.00e-16,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[5][j-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,3.92e-03,1.00e-16,1.32e-03,2.49e-03,
		  4.18e-03,1.00e-16,8.83e-03,2.21e-02,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[6][j-1] = _itmp6[_r++];
			}
		}
		{ static double _itmp7[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,1.41e-03,1.00e-16,1.82e-03,
		  1.18e-03,4.75e-03,1.68e-03,1.29e-02,1.61e-02};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[7][j-1] = _itmp7[_r++];
			}
		}
		{ static double _itmp8[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,1.00e-16,1.00e-16,
		  2.12e-03,2.47e-03,1.00e-16,3.65e-03,2.13e-02};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[8][j-1] = _itmp8[_r++];
			}
		}
		{ static double _itmp9[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,2.56e-03,
		  1.00e-16,1.00e-16,2.23e-03,2.25e-03,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[9][j-1] = _itmp9[_r++];
			}
		}
		{ static double _itmp10[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  1.36e-03,1.00e-16,2.70e-04,1.00e-16,2.76e-03};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[10][j-1] = _itmp10[_r++];
			}
		}
		{ static double _itmp11[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,3.71e-04,7.70e-05,1.60e-04,2.50e-04};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[11][j-1] = _itmp11[_r++];
			}
		}
		{ static double _itmp12[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,1.00e-16,1.50e-04,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[12][j-1] = _itmp12[_r++];
			}
		}
		{ static double _itmp13[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,1.90e-04,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[13][j-1] = _itmp13[_r++];
			}
		}
		{ static double _itmp14[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,5.50e-04};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[14][j-1] = _itmp14[_r++];
			}
		}
		for( j=1; j <= NLFE2; j++ )
		{
			Fe2A[15][j-1] = 0.00e00;
			}
		for( i=1; i <= NLFE2; i++ )
		{
			Fe2CS[0][i-1] = 0.000e00;
			}
		{ static double _itmp15[] = {5.521e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[1][i-1] = _itmp15[_r++];
			}
		}
		{ static double _itmp16[] = {1.491e00,5.442e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[2][i-1] = _itmp16[_r++];
			}
		}
		{ static double _itmp17[] = {6.747e-01,1.197e00,4.319e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[3][i-1] = _itmp17[_r++];
			}
		}
		{ static double _itmp18[] = {2.839e-01,4.191e-01,7.255e-01,
		  2.640e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[4][i-1] = _itmp18[_r++];
			}
		}
		{ static double _itmp19[] = {3.595e00,1.633e00,6.373e-01,2.108e-01,
		  5.804e-02,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[5][i-1] = _itmp19[_r++];
			}
		}
		{ static double _itmp20[] = {1.507e00,1.484e00,1.113e00,6.083e-01,
		  2.231e-01,8.365e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[6][i-1] = _itmp20[_r++];
			}
		}
		{ static double _itmp21[] = {4.969e-01,1.026e00,9.581e-01,7.646e-01,
		  4.466e-01,2.227e00,7.386e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[7][i-1] = _itmp21[_r++];
			}
		}
		{ static double _itmp22[] = {1.371e-01,4.457e-01,7.373e-01,
		  7.140e-01,4.239e-01,7.005e-01,2.069e00,5.627e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[8][i-1] = _itmp22[_r++];
			}
		}
		{ static double _itmp23[] = {1.098e01,2.479e00,5.033e-01,1.468e-01,
		  5.434e-02,2.106e00,8.906e-01,2.935e-01,6.481e-02,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[9][i-1] = _itmp23[_r++];
			}
		}
		{ static double _itmp24[] = {5.599e-01,6.540e00,2.847e00,6.264e-01,
		  1.177e-01,7.248e-01,8.511e-01,6.858e-01,3.146e-01,2.230e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[10][i-1] = _itmp24[_r++];
			}
		}
		{ static double _itmp25[] = {1.909e-01,3.552e-01,3.583e00,2.445e00,
		  5.540e-01,1.477e-01,5.556e-01,4.996e-01,5.181e-01,8.177e-01,
		  1.752e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[11][i-1] = _itmp25[_r++];
			}
		}
		{ static double _itmp26[] = {6.104e-02,1.227e-01,1.956e-01,
		  1.534e00,1.649e00,2.123e-02,1.471e-01,3.672e-01,3.274e-01,
		  3.689e-01,3.784e-01,1.015e00,0.000e00,0.000e00,0.000e00,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[12][i-1] = _itmp26[_r++];
			}
		}
		{ static double _itmp27[] = {9.480e-01,6.813e-01,3.850e-01,
		  1.707e-01,5.327e-02,1.232e00,6.882e-01,3.779e-01,1.971e-01,
		  1.008e00,5.870e-01,2.563e-01,7.694e-02,0.000e00,0.000e00,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[13][i-1] = _itmp27[_r++];
			}
		}
		{ static double _itmp28[] = {5.024e-01,2.276e-01,2.681e-01,
		  2.971e-01,1.967e-01,4.114e-01,5.541e-01,4.492e-01,2.809e-01,
		  4.690e-01,2.905e-01,3.180e-01,2.262e-01,1.280e00,0.000e00,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[14][i-1] = _itmp28[_r++];
			}
		}
		{ static double _itmp29[] = {3.278e-02,2.871e-01,2.451e-01,
		  1.319e-01,4.768e-02,1.404e-01,1.990e-01,2.603e-01,2.482e-01,
		  7.626e-02,2.877e-01,2.033e-01,8.596e-02,4.837e-01,6.294e-01,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[15][i-1] = _itmp29[_r++];
			}
		}
		_aini = 0;
	}

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

	/* number of levels in model ion */

	/* statistical weights */
	/* excitation energies in Kelvin
	 * excitation energies in Kelvin */

	/* each set is for a lower level indicated by second element in array,
	 * index runs over upper level
	 * A's are saved into arrays as AulEscp[up][lo] */

	/* collision strengths from 
	 * >>refer	fe2	cs	Zhang, H.L., & Pradhan, A., 1995, A&A, 293, 953
	 * order is cs(low,up), for 10,000 K */
	if( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* allocate the 1D arrays*/
		if( (pops = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
			BadMalloc();
		if( (create = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
			BadMalloc();
		if( (destroy = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
			BadMalloc();
		if( (depart = (double *)MALLOC( sizeof(double)*(NLFE2) )) == NULL )
			BadMalloc();
		/* create space for the 2D arrays */
		if( (AulPump = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulDest = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulEscp = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (col_str = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE2); ++i )
		{
			if( (AulPump[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulDest[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulEscp[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (col_str[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
		}
	}

	/* bail if no ions */
	if( dense.xIonDense[ipIRON][1] <= 0. )
	{
		/* these are only predicted by the very large feii atom */
		FeII.fe25to6 = 0.;
		FeII.fe27to7 = 0.;
		FeII.fe28to8 = 0.;
		FeII.fe29to9 = 0.;
 		FeII.fe32to6 = 0.;
 		FeII.fe33to7 = 0.;
 		FeII.fe37to7 = 0.;
 		FeII.fe39to8 = 0.;
 		FeII.fe40to9 = 0.;
 		FeII.fe37to6 = 0.;
 		FeII.fe39to7 = 0.;
 		FeII.fe40to8 = 0.;
 		FeII.fe41to9 = 0.;
 		FeII.fe39to6 = 0.;
 		FeII.fe40to7 = 0.;
 		FeII.fe41to8 = 0.;
		FeII.fe42to6 = 0.;
		FeII.fe43to7 = 0.;
		FeII.fe42to7 = 0.;
		FeII.fe80to28 =0.;
		FeII.fe36to2 = 0.;
		FeII.fe36to3 = 0.;
		FeII.fe32to1 = 0.;
		FeII.fe33to2 = 0.;
		FeII.fe36to5 = 0.;
		FeII.fe32to2 = 0.;
		FeII.fe33to3 = 0.;
		FeII.fe30to3 = 0.;
		FeII.fe33to6 = 0.;
		FeII.fe24to2 = 0.;
		FeII.fe32to7 = 0.;
		FeII.fe35to8 = 0.;
		FeII.fe34to8 = 0.;
		FeII.fe27to6 = 0.;
		FeII.fe28to7 = 0.;
		FeII.fe30to8 = 0.;
		FeII.fe24to6 = 0.;
		FeII.fe29to8 = 0.;
		FeII.fe24to7 = 0.;
		FeII.fe22to7 = 0.;
		FeII.fe38to11 =0.;
		FeII.fe19to8 = 0.;
		FeII.fe17to6 = 0.;
		FeII.fe18to7 = 0.;
		FeII.fe18to8 = 0.;
		FeII.fe17to7 = 0.;
		FeII.fe21406 = 0.;
		FeII.fe21507 = 0.;
		FeII.fe21508 = 0.;
		FeII.fe21609 = 0.;
		FeII.fe21308 = 0.;
		FeII.fe21207 = 0.;
		FeII.fe21106 = 0.;
		FeII.fe21006 = 0.;
		FeII.fe21204 = 0.;
		FeII.fe21103 = 0.;
		FeII.fe21104 = 0.;
		FeII.fe21001 = 0.;
		FeII.fe21002 = 0.;
		FeII.fe20201 = 0.;
		FeII.fe20302 = 0.;
		FeII.fe20706 = 0.;
		FeII.fe20807 = 0.;
		FeII.fe20908 = 0.;
		FeII.fe21007 = 0.;
		FeII.fe21107 = 0.;
		FeII.fe21108 = 0.;
		FeII.fe21110 = 0.;
		FeII.fe21208 = 0.;
		FeII.fe21209 = 0.;
		FeII.fe21211 = 0.;

		CoolAdd("Fe 2",2,0.);

		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( NLFE2 <= LIMLEVELN);
		for( i=0; i < NLFE2; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->Fe2Lev16()\n", debug_fp );
#		endif
		return;
	}

	/* all elements are used, and must be set to zero if zero */
	for( i=0; i < NLFE2; i++ )
	{
		create[i] = 0.;
		destroy[i] = 0.;
		for( j=0; j < NLFE2; j++ )
		{
			/*data[j][i] = 1e33;*/
			col_str[j][i] = 0.;
			AulEscp[j][i] = 0.;
			AulDest[j][i] = 0.;
			AulPump[j][i] = 0.;
		}
	}

	/* fill in einstein As and collision strengths */
	for( i=0; i < NLFE2; i++ )
	{
		for( j=i + 1; j < NLFE2; j++ )
		{
			AulEscp[j][i] = Fe2A[i][j];
			col_str[j][i] = Fe2CS[j][i];
		}
	}

	/* atom_levelN itself is well-protected against zero abundances, and low temperatures
	 * >>chng 96 oct 16, reduce number of levels for low temp since
	 * can still be important coolant when highest level not populated */

	if( phycon.te < 500. )
	{
		nUsed = 5;
	}

	else if( phycon.te < 1e3 )
	{
		nUsed = 9;
	}

	else
	{
		nUsed = NLFE2;
	}

	/*void atom_levelN(long,long,float,double[],double[],double[],double*,
	double*,double*,long*,float*,float*,STRING,int*);*/
	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(nUsed,
		dense.xIonDense[ipIRON][1],
		gfe2,
		ex,
		'K',
		pops,
		depart,
		&AulEscp ,
		&col_str ,
		&AulDest,
		&AulPump,
		&CollRate,
		create,
		destroy,
		FALSE,/* say atom_levelN should evaluate coll rates from cs */
		/*&ipdest,*/
		&FeII.Fe2_16levl_cool,
		&FeII.ddT_Fe2_16levl_cool,
		"FeII",
		&lgNegPop,
		&lgZeroPop,
		FALSE );

	/* LIMLEVELN is the dim of the PopLevels vector */
	ASSERT( NLFE2 <= LIMLEVELN );
	for( i=0; i<nUsed; ++i)
	{
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}
	for( i=nUsed; i<NLFE2; ++i)
	{
		pops[i] = 0.;
		depart[i] = 1.;
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}

	if( lgNegPop>0 )
	{
		fprintf( ioQQQ, " fe2levl6 found negative populations\n" );
		ShowMe();
		puts( "[Stop in fe2levl6]" );
		cdEXIT(EXIT_FAILURE);
	}

#	if 0
	/* >>chng 05 dec 03, move to following call to this routine */
	/* cooling due to lowest 16 levels of FeII - evaluated in call to atomn */
	CoolAdd("Fe 2",2,FeII.Fe2_16levl_cool);

	thermal.dCooldT += dfe2dt;
#	endif

	/* all 10 transitions within ground term */
	FeII.fe21308 = pops[12]*Fe2A[7][12]*(ex[12] - ex[7])*BOLTZMANN;
	FeII.fe21207 = pops[11]*Fe2A[6][11]*(ex[11] - ex[6])*BOLTZMANN;
	FeII.fe21106 = pops[10]*Fe2A[5][10]*(ex[10] - ex[5])*BOLTZMANN;
	FeII.fe21006 = pops[9]*Fe2A[5][9]*(ex[9] - ex[5])*BOLTZMANN;
	FeII.fe21204 = pops[11]*Fe2A[3][11]*(ex[11] - ex[3])*BOLTZMANN;
	FeII.fe21103 = pops[10]*Fe2A[2][10]*(ex[10] - ex[2])*BOLTZMANN;
	FeII.fe21104 = pops[10]*Fe2A[3][10]*(ex[10] - ex[3])*BOLTZMANN;
	FeII.fe21001 = pops[9]*Fe2A[0][9]*(ex[9] - ex[0])*BOLTZMANN;
	FeII.fe21002 = pops[9]*Fe2A[1][9]*(ex[9] - ex[1])*BOLTZMANN;
	FeII.fe20201 = pops[1]*Fe2A[0][1]*(ex[1] - ex[0])*BOLTZMANN;
	FeII.fe20302 = pops[2]*Fe2A[1][2]*(ex[2] - ex[1])*BOLTZMANN;
	FeII.fe20706 = pops[6]*Fe2A[5][6]*(ex[6] - ex[5])*BOLTZMANN;
	FeII.fe20807 = pops[7]*Fe2A[6][7]*(ex[7] - ex[6])*BOLTZMANN;
	FeII.fe20908 = pops[8]*Fe2A[7][8]*(ex[8] - ex[7])*BOLTZMANN;
	FeII.fe21007 = pops[9]*Fe2A[6][9]*(ex[9] - ex[6])*BOLTZMANN;
	FeII.fe21107 = pops[10]*Fe2A[6][10]*(ex[10] - ex[6])*BOLTZMANN;
	FeII.fe21108 = pops[10]*Fe2A[7][10]*(ex[10] - ex[7])*BOLTZMANN;
	FeII.fe21110 = pops[10]*Fe2A[9][10]*(ex[10] - ex[9])*BOLTZMANN;
	FeII.fe21208 = pops[11]*Fe2A[7][11]*(ex[11] - ex[7])*BOLTZMANN;
	FeII.fe21209 = pops[11]*Fe2A[8][11]*(ex[11] - ex[8])*BOLTZMANN;
	FeII.fe21211 = pops[11]*Fe2A[10][11]*(ex[11] - ex[10])*BOLTZMANN;

	/* these are the highest lines in the 16 level atom, and so most uncertain */
	FeII.fe21406 = pops[13]*Fe2A[5][13]*(ex[13] - ex[5])*BOLTZMANN;
	FeII.fe21507 = pops[14]*Fe2A[6][14]*(ex[14] - ex[6])*BOLTZMANN;
	FeII.fe21508 = pops[14]*Fe2A[7][14]*(ex[14] - ex[7])*BOLTZMANN;
	FeII.fe21609 = pops[15]*Fe2A[8][15]*(ex[15] - ex[8])*BOLTZMANN;

	/* these are only predicted by the very large feii atom */
	FeII.fe25to6 = 0.;
	FeII.fe27to7 = 0.;
	FeII.fe28to8 = 0.;
	FeII.fe29to9 = 0.;
 	FeII.fe32to6 = 0.;
 	FeII.fe33to7 = 0.;
 	FeII.fe37to7 = 0.;
 	FeII.fe39to8 = 0.;
 	FeII.fe40to9 = 0.;
 	FeII.fe37to6 = 0.;
 	FeII.fe39to7 = 0.;
 	FeII.fe40to8 = 0.;
 	FeII.fe41to9 = 0.;
 	FeII.fe39to6 = 0.;
 	FeII.fe40to7 = 0.;
 	FeII.fe41to8 = 0.;
	FeII.fe42to6 = 0.;
	FeII.fe43to7 = 0.;
	FeII.fe42to7 = 0.;
	FeII.fe80to28 =0.;
	FeII.fe36to2 = 0.;
	FeII.fe36to3 = 0.;
	FeII.fe32to1 = 0.;
	FeII.fe33to2 = 0.;
	FeII.fe36to5 = 0.;
	FeII.fe32to2 = 0.;
	FeII.fe33to3 = 0.;
	FeII.fe30to3 = 0.;
	FeII.fe33to6 = 0.;
	FeII.fe24to2 = 0.;
	FeII.fe32to7 = 0.;
	FeII.fe35to8 = 0.;
	FeII.fe34to8 = 0.;
	FeII.fe27to6 = 0.;
	FeII.fe28to7 = 0.;
	FeII.fe30to8 = 0.;
	FeII.fe24to6 = 0.;
	FeII.fe29to8 = 0.;
	FeII.fe24to7 = 0.;
	FeII.fe22to7 = 0.;
	FeII.fe38to11 =0.;
	FeII.fe19to8 = 0.;
	FeII.fe17to6 = 0.;
	FeII.fe18to7 = 0.;
	FeII.fe18to8 = 0.;
	FeII.fe17to7 = 0.;

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

/*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
static void Fe7Lev8(void)
{
	int lgNegPop,
		lgZeroPop;
	double scale;
	long int i, 
	  j;
	static int lgFirst=TRUE;
	static long int ipPump=-1;

	double dfe7dt,
		FUV_pump;

	long int ihi , ilo;
	/*static long int **ipdest; */
	static double 
	  *depart,
	  *pops, 
	  *destroy, 
	  *create ,
	  **AulDest, 
	  **CollRate,
	  **AulPump,
	  **AulNet, 
	  **col_str;
	/*
	following are FeVII lines predicted in limit_laser_2000
	Fe 7  5721A -24.399   0.1028 
	Fe 7  4989A -24.607   0.0637 
	Fe 7  4894A -24.838   0.0374
	Fe 7  4699A -25.693   0.0052
	Fe 7  6087A -24.216   0.1566
	Fe 7  5159A -24.680   0.0538
	Fe 7  4943A -25.048   0.0231
	Fe 7  3587A -24.285   0.1336
	Fe 7  5277A -25.053   0.0228
	Fe 7  3759A -24.139   0.1870
	Fe 7 3.384m -25.621   0.0062
	Fe 7 2.629m -25.357   0.0113
	Fe 7 9.510m -24.467   0.0878
	Fe 7 7.810m -24.944   0.0293
	*/

	/* statistical weights for the n levels */
	static double gfe7[NLFE7]={5.,7.,9.,5.,1.,3.,5.,9.};

	/*refer	Fe7	energies	Ekbert, J.O. 1981, Phys. Scr 23, 7 */
	/*static double ex[NLFE7]={0.,1047.,2327.,17475.,20037.,20428.,21275.,28915.};*/
	/* excitation energies in wavenumbers 
	 * >>chng 05 dec 15, rest of code had assumed that there were energies in Kelvin
	 * rather than wavenumbers.  Bug caught by Kevin Blagrave 
	 * modified atom_levelN to accept either Kelvin or wavenumbers */
	static double excit_wn[NLFE7]={0. , 1051.5 , 2331.5 , 17475.5 , 20040.3 , 20430.1 , 21278.6 , 28927.3 };

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

	if( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* allocate the 1D arrays*/
		/* create space for the 2D arrays */
		if( (depart = ((double *)MALLOC((NLFE7)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (pops = ((double *)MALLOC((NLFE7)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (destroy = ((double *)MALLOC((NLFE7)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (create = ((double *)MALLOC((NLFE7)*sizeof(double)))) == NULL )
			BadMalloc();
		/* now the 2-d arrays */
		if( (fe.Fe7_wl = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (fe.Fe7_emiss = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulNet = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (col_str = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulPump = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulDest = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE7); ++i )
		{
			if( (fe.Fe7_wl[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (fe.Fe7_emiss[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulNet[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (col_str[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulPump[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulDest[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
		}

		/* set some to constant values after zeroing out */
		for( i=0; i<NLFE7; ++i )
		{
			create[i] = 0.;
			destroy[i] = 0.;
			for( j=0; j<NLFE7; ++j )
			{
				AulNet[i][j] = 0.;
				col_str[i][j] = 0.;
				CollRate[i][j] = 0.;
				AulDest[i][j] = 0.;
				AulPump[i][j] = 0.;
				fe.Fe7_wl[i][j] = 0.;
				fe.Fe7_emiss[i][j] = 0.;
			}
		}
		SET_DBL_SNAN( fe.Fe7_wl[2][1] );
		SET_DBL_SNAN( fe.Fe7_emiss[2][1] );
		SET_DBL_SNAN( fe.Fe7_wl[1][0] );
		SET_DBL_SNAN( fe.Fe7_emiss[1][0] );
		/* two levels within ground term must never be addressed in this array -
		 * they are fully transferred */
		for( ilo=0; ilo<NLFE7-1; ++ilo )
		{
			/* must not do 1-0 or 2-1, which are transferred lines */
			for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
			{
				fe.Fe7_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / RefIndex( excit_wn[ihi]-excit_wn[ilo] );
			}
		}

		/* assume FeVII is optically thin - just use As as net escape */
		/*>>refer	Fe7	As	Berrington, K.A., Nakazaki, S., & Norrington, P.H. 2000, A&AS, 142, 313 */
		AulNet[1][0] = 0.0325;
		AulNet[2][0] = 0.167e-8;
		/* following corrected from 3.25 to 0.372 as per Keith Berrington email */
		AulNet[3][0] = 0.372;
		AulNet[4][0] = 0.135;
		AulNet[5][0] = 0.502e-1;
		/* following corrected from 0.174 to 0.0150 as per Keith Berrington email */
		AulNet[6][0] = 0.0150;
		AulNet[7][0] = 0.959e-3;

		AulNet[2][1] = 0.466e-1;
		AulNet[3][1] = 0.603;
		AulNet[5][1] = 0.762e-1;
		AulNet[6][1] = 0.697e-1;
		AulNet[7][1] = 0.343;

		AulNet[3][2] = 0.139e-2;
		AulNet[6][2] = 0.735e-1;
		AulNet[7][2] = 0.503;

		AulNet[4][3] = 0.472e-6;
		AulNet[5][3] = 0.572e-1;
		/* following corrected from 0.191 to 0.182 as per Keith Berrington email */
		AulNet[6][3] = 0.182;
		AulNet[7][3] = 0.414e-2;

		AulNet[5][4] = 0.115e-2;
		AulNet[6][4] = 0.139e-7;

		AulNet[6][5] = 0.743e-2;

		AulNet[7][6] = 0.454e-4;

		/* collision strengths at log T 4.5 */
		/*TODO	2	put in temperature dependence */
		/*>>refer	Fe7	cs	Berrington, K.A., Nakazaki, S., & Norrington, P.H. 2000, A&AS, 142, 313 */
#		if 0
		if( fudge(-1) )
		{
			fprintf(ioQQQ,"DEBUG fudge call cool_iron\n");
			scale = fudge(0);
		}
		else
#		endif
		scale = 1.;

		col_str[1][0] = 3.35;
		col_str[2][0] = 1.17;
		col_str[3][0] = 0.959;
		col_str[4][0] = 0.299;
		col_str[5][0] = 0.633;
		col_str[6][0] = 0.549;
		col_str[7][0] = 1.24*scale;

		col_str[2][1] = 4.11;
		col_str[3][1] = 1.29;
		col_str[4][1] = 0.235;
		col_str[5][1] = 0.833;
		col_str[6][1] = 1.06;
		col_str[7][1] = 1.74*scale;

		col_str[3][2] = 1.60;
		col_str[4][2] = 0.187;
		col_str[5][2] = 0.690;
		col_str[6][2] = 1.94;
		col_str[7][2] = 2.25*scale;

		col_str[4][3] = 0.172;
		col_str[5][3] = 0.531;
		col_str[6][3] = 1.06;
		col_str[7][3] = 2.02;

		col_str[5][4] = 0.370;
		col_str[6][4] = 0.324;
		col_str[7][4] = 0.164;

		col_str[6][5] = 1.17;
		col_str[7][5] = 0.495;

		col_str[7][6] = 0.903;

		/* check whether level 2 lines are on, and if so, find the driving lines that
		 * can pump the upper levels of this atom */
		if( nWindLine > 0 )
		{
			ipPump = -1;
			for( i=0; i<nWindLine; ++i )
			{
				/* don't test on nelem==ipIRON since lines on physics, not C, scale */
				if( TauLine2[i].nelem ==26 && TauLine2[i].IonStg==7 &&
					/* >>chng 05 jul 10, move line to wavelength that agrees with nist tables
					 * here and in level2 data file 
					 * NB must add a few wn to second number to get hit -
					 * logic is that this is lowest E1 transition */
					(TauLine2[i].EnergyWN-4.27360E+05) < 0. )
				{
					ipPump = i;
					break;
				}
			}
			if( ipPump<0 )
			{
				fprintf(ioQQQ,"PROBLEM Fe7Lev8 cannot identify the FUV driving line.\n");
				TotalInsanity();
			}
		}
		else
		{
			ipPump = 0;
		}
	}

	/* bail if no ions */
	if( dense.xIonDense[ipIRON][6] <= 0. )
	{
		CoolAdd("Fe 7",0,0.);

		fe.Fe7CoolTot = 0.;
		for( ilo=0; ilo<NLFE7-1; ++ilo )
		{
			/* must not do 1-0 or 2-1, which are transferred lines */
			for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
			{
				fe.Fe7_emiss[ihi][ilo] = 0.;
			}
		}
		TauLines[ipFe0795].PopLo = 0.;
		TauLines[ipFe0795].PopOpc = 0.;
		TauLines[ipFe0795].PopHi = 0.;
		TauLines[ipFe0795].xIntensity = 0.;
		TauLines[ipFe0795].cool = 0.;
		TauLines[ipFe0795].phots = 0.;
		TauLines[ipFe0795].ColOvTot = 0.;
		TauLines[ipFe0795].heat = 0.;
		CoolAdd( "Fe 7", TauLines[ipFe0795].WLAng , 0.);
		TauLines[ipFe0778].PopLo = 0.;
		TauLines[ipFe0778].PopOpc = 0.;
		TauLines[ipFe0778].PopHi = 0.;
		TauLines[ipFe0778].xIntensity = 0.;
		TauLines[ipFe0778].cool = 0.;
		TauLines[ipFe0778].phots = 0.;
		TauLines[ipFe0778].ColOvTot = 0.;
		TauLines[ipFe0778].heat = 0.;
		CoolAdd( "Fe 7", TauLines[ipFe0778].WLAng , 0.);
		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( NLFE7 <= LIMLEVELN);
		for( i=0; i < NLFE7; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->Fe7Lev8()\n", debug_fp );
#		endif
		return;
	}

	/* do pump rate for FUV excitation of 3P (levels 5-7 on physics scale, not C scale) */
	if( ipPump )
	{
		FUV_pump = TauLine2[ipPump].pump * 0.3 /(0.3+TauLine2[ipPump].Pesc);
	}
	else
	{
		FUV_pump = 0.;
	}
	/*if( fudge(-1) )
		FUV_pump *= fudge(0);*/
	/*fprintf(ioQQQ,"DEBUG pump %li %.2e %.2e %.2e\n", 
		nzone, FUV_pump , dense.xIonDense[ipIRON][6],TauLine2[ipPump].TauIn );*/

	/* this represents photoexcitation of 3P from ground level
	 * >>chng 04 nov 22, upper array elements were on physics not c scale, off by one, TE */
	AulPump[0][4] = FUV_pump;
	AulPump[1][4] = FUV_pump;
	AulPump[2][4] = FUV_pump;
	AulPump[0][5] = FUV_pump;
	AulPump[1][5] = FUV_pump;
	AulPump[2][5] = FUV_pump;
	AulPump[0][6] = FUV_pump;
	AulPump[1][6] = FUV_pump;
	AulPump[2][6] = FUV_pump;
	/* these were set in the previous call to atom_levelN as the previous pump times
	 * ratio of staticial weights, so this is the downward pump rate */
	AulPump[4][0] = 0;
	AulPump[4][1] = 0;
	AulPump[4][2] = 0;
	AulPump[5][0] = 0;
	AulPump[5][1] = 0;
	AulPump[5][2] = 0;
	AulPump[6][0] = 0;
	AulPump[6][1] = 0;
	AulPump[6][2] = 0;

	/* within ground term update to rt results */
	AulNet[1][0] = TauLines[ipFe0795].Aul*(TauLines[ipFe0795].Pesc + TauLines[ipFe0795].Pelec_esc);
	AulDest[1][0] = TauLines[ipFe0795].Aul*TauLines[ipFe0795].Pdest;
	AulPump[0][1] = TauLines[ipFe0795].pump;
	AulPump[1][0] = 0.;

	AulNet[2][1] = TauLines[ipFe0778].Aul*(TauLines[ipFe0778].Pesc + TauLines[ipFe0778].Pelec_esc);
	AulDest[2][1] = TauLines[ipFe0778].Aul*TauLines[ipFe0778].Pdest;
	AulPump[1][2] = TauLines[ipFe0778].pump;
	AulPump[2][1] = 0.;

	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(
		/* number of levels */
		NLFE7,
		/* the abundance of the ion, cm-3 */
		dense.xIonDense[ipIRON][6],
		/* the statistical weights */
		gfe7,
		/* the excitation energies in wavenumbers */
		excit_wn,
		/* units are wavenumbers */
		'w',
		/* the derived populations - cm-3 */
		pops,
		/* the derived departure coefficients */
		depart,
		/* the net emission rate, Aul * escp prob */
		&AulNet ,
		/* the collision strengths */
		&col_str ,
		/* A * dest prob */
		&AulDest,
		/* pumping rate */
		&AulPump,
		/* collision rate, s-1, must defined if no collision strengths */
		&CollRate,
		/* creation vector */
		create,
		/* destruction vector */
		destroy,
		/* say atom_levelN should evaluate coll rates from cs */
		FALSE,
		&fe.Fe7CoolTot,
		&dfe7dt,
		"Fe 7",
		&lgNegPop,
		&lgZeroPop,
		FALSE );

	/* LIMLEVELN is the dim of the PopLevels vector */
	ASSERT( NLFE7 <= LIMLEVELN );
	for( i=0; i<NLFE7; ++i)
	{
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}

	if( lgZeroPop )
	{
		/* this case, too cool to excite upper levels of atom, but will want to do ir lines
		 * in ground term */
		PutCS(col_str[1][0],&TauLines[ipFe0795]);
		PutCS(col_str[2][1],&TauLines[ipFe0778]);
		PutCS(col_str[2][0],&TauDummy);
		atom_level3(&TauLines[ipFe0795],&TauLines[ipFe0778],&TauDummy);
		atoms.PopLevels[0] = TauLines[ipFe0795].PopLo;
		atoms.PopLevels[1] = TauLines[ipFe0795].PopHi;
		atoms.PopLevels[2] = TauLines[ipFe0778].PopHi;
		for( ilo=0; ilo<NLFE7-1; ++ilo )
		{
			/* must not do 1-0 or 2-1, which are transferred lines */
			for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
			{
				fe.Fe7_emiss[ihi][ilo] = 0.;
			}
		}
	}
	else
	{
		/* this case non-zero pops, but must set up vars within emline structure */
		TauLines[ipFe0795].PopLo = atoms.PopLevels[0];
		TauLines[ipFe0795].PopHi = atoms.PopLevels[1];
		TauLines[ipFe0795].PopOpc = (pops[0] - pops[1]*gfe7[0]/gfe7[1]);
		TauLines[ipFe0795].xIntensity = TauLines[ipFe0795].phots*TauLines[ipFe0795].EnergyErg;
		TauLines[ipFe0795].phots = TauLines[ipFe0795].Aul*(TauLines[ipFe0795].Pesc + TauLines[ipFe0795].Pelec_esc)*pops[1];
		TauLines[ipFe0795].ColOvTot = (float)(CollRate[0][1]/(CollRate[0][1]+TauLines[ipFe0795].pump) );
		TauLines[ipFe0795].cool = 0.;
		TauLines[ipFe0795].heat = 0.;

		TauLines[ipFe0778].PopLo = atoms.PopLevels[1];
		TauLines[ipFe0778].PopHi = atoms.PopLevels[2];
		TauLines[ipFe0778].PopOpc = (pops[1] - pops[2]*gfe7[1]/gfe7[2]);
		TauLines[ipFe0778].xIntensity = TauLines[ipFe0778].phots*TauLines[ipFe0778].EnergyErg;
		TauLines[ipFe0778].phots = TauLines[ipFe0778].Aul*(TauLines[ipFe0778].Pesc + TauLines[ipFe0778].Pelec_esc)*pops[2];
		TauLines[ipFe0778].ColOvTot = (float)(CollRate[1][2]/(CollRate[1][2]+TauLines[ipFe0778].pump) );
		TauLines[ipFe0778].heat = 0.;
		TauLines[ipFe0778].cool = 0.;
	}

	if( lgNegPop>0 )
	{
		fprintf( ioQQQ, "PROBLEM Fe7Lev8 found negative populations\n" );
		ShowMe();
		puts( "[Stop in fe7lev8]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* add cooling then its deriv */
	CoolAdd("Fe 7",0,fe.Fe7CoolTot);
	/* deriv of cooling */
	thermal.dCooldT += dfe7dt;

	/* find emission in each line */
	for( ilo=0; ilo<NLFE7-1; ++ilo )
	{
		/* must not do 1-0 or 2-1, which are transferred lines */
		for( ihi=MAX2(3,ilo+1); ihi<NLFE7; ++ihi )
		{
			/* emission in these lines */
			fe.Fe7_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*ERG1CM;
		}
	}

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

/*Fe3Lev14 compute populations and cooling due to 14 level Fe III ion 
 * >>chng 05 dec 17, code provided by Kevin Blagrave */
static void Fe3Lev14(void)
{
	int lgNegPop,
		lgZeroPop;
	long int i,
		j;
	static int lgFirst=TRUE;

	double dfe3dt;

	long int ihi , ilo;
	static double
		*depart,
		*pops,
		*destroy,
		*create ,
		**AulDest,
		**CollRate,
		**AulPump,
		**AulNet,
		**col_str;

	/* statistical weights for the n levels */
	static double gfe3[NLFE3]={9.,7.,5.,3.,1.,5.,13.,11.,9.,3.,1.,9.,7.,5.};

	/*refer Fe3	energies	NIST version 3 Atomic Spectra Database */
	/* from smallest to largest energy (not in multiplet groupings) */

	/* energy in wavenumbers */
	static double excit_wn[NLFE3]={
		0.0    ,   436.2,   738.9,   932.4,  1027.3,
		19404.8, 20051.1, 20300.8, 20481.9, 20688.4,
		21208.5, 21462.2, 21699.9, 21857.2 };

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

	if( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* allocate the 1D arrays*/
		/* create space for the 2D arrays */
		if( (depart = ((double *)MALLOC((NLFE3)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (pops = ((double *)MALLOC((NLFE3)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (destroy = ((double *)MALLOC((NLFE3)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (create = ((double *)MALLOC((NLFE3)*sizeof(double)))) == NULL )
			BadMalloc();
		/* now the 2-d arrays */
		if( (fe.Fe3_wl = ((double **)MALLOC((NLFE3)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (fe.Fe3_emiss = ((double **)MALLOC((NLFE3)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulNet = ((double **)MALLOC((NLFE3)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (col_str = ((double **)MALLOC((NLFE3)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulPump = ((double **)MALLOC((NLFE3)*sizeof(double *)))) == NULL )
			BadMalloc();  
		if( (CollRate = ((double **)MALLOC((NLFE3)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulDest = ((double **)MALLOC((NLFE3)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE3); ++i )
		{
			if( (fe.Fe3_wl[i] = ((double *)MALLOC((NLFE3)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (fe.Fe3_emiss[i] = ((double *)MALLOC((NLFE3)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulNet[i] = ((double *)MALLOC((NLFE3)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (col_str[i] = ((double *)MALLOC((NLFE3)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulPump[i] = ((double *)MALLOC((NLFE3)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE3)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulDest[i] = ((double *)MALLOC((NLFE3)*sizeof(double )))) == NULL )
				BadMalloc();
		}

		/* set some to constant values after zeroing out */
		for( i=0; i<NLFE3; ++i )
		{
			create[i] = 0.;
			destroy[i] = 0.;
			for( j=0; j<NLFE3; ++j )
			{
				AulNet[i][j] = 0.;
				col_str[i][j] = 0.;
				CollRate[i][j] = 0.;
				AulDest[i][j] = 0.;
				AulPump[i][j] = 0.;
				fe.Fe3_wl[i][j] = 0.;
				fe.Fe3_emiss[i][j] = 0.;
			}
		}
		/* calculates wavelengths of transitions */
		/* dividing by RefIndex converts the vacuum wavelength to air wavelength */
		for( ihi=1; ihi<NLFE3; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				fe.Fe3_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / RefIndex( (excit_wn[ihi]-excit_wn[ilo]) );
			}
		}

		/* assume FeIII is optically thin - just use As as net escape */
		/*>>refer	Fe3	as	Quinet, P., 1996, A&AS, 116, 573      */
		AulNet[1][0] = 2.8e-3;
		AulNet[7][0] = 4.9e-6;
		AulNet[8][0] = 5.7e-3;
		AulNet[11][0] = 4.5e-1;
		AulNet[12][0] = 4.2e-2;

		AulNet[2][1] = 1.8e-3;
		AulNet[5][1] = 4.2e-1;
		AulNet[8][1] = 1.0e-3;
		AulNet[11][1] = 8.4e-2;
		AulNet[12][1] = 2.5e-1;
		AulNet[13][1] = 2.7e-2;

		AulNet[3][2] = 7.0e-4;
		AulNet[5][2] = 5.1e-5;
		AulNet[9][2] = 5.4e-1;
		AulNet[12][2] = 8.5e-2;
		AulNet[13][2] = 9.8e-2;

		AulNet[4][3] = 1.4e-4;
		AulNet[5][3] = 3.9e-2;
		AulNet[9][3] = 4.1e-5;
		AulNet[10][3] = 7.0e-1;
		AulNet[13][3] = 4.7e-2;

		AulNet[9][4] = 9.3e-2;

		AulNet[9][5] = 4.7e-2;
		AulNet[12][5] = 2.5e-6;
		AulNet[13][5] = 1.7e-5;

		AulNet[7][6] = 2.7e-4;

		AulNet[8][7] = 1.2e-4;
		AulNet[11][7] = 6.6e-4;

		AulNet[11][8] = 1.6e-3;
		AulNet[12][8] = 7.8e-4;

		AulNet[10][9] = 8.4e-3;
		AulNet[13][9] = 2.8e-7;

		AulNet[12][11] = 3.0e-4;

		AulNet[13][12] = 1.4e-4;

		/* collision strengths at log T 4 */
		/*TODO	2	put in temperature dependence */
		/*>>refer	Fe3	cs	Zhang, H.  1996, 119, 523 */
		col_str[1][0] = 2.92;
		col_str[2][0] = 1.24;
		col_str[3][0] = 0.595;
		col_str[4][0] = 0.180;
		col_str[5][0] = 0.580;
		col_str[6][0] = 1.34; 
		col_str[7][0] = 0.489;
		col_str[8][0] = 0.0926;
		col_str[9][0] = 0.165;
		col_str[10][0] = 0.0213;
		col_str[11][0] = 1.07;
		col_str[12][0] = 0.435;
		col_str[13][0] = 0.157;

		col_str[2][1] = 2.06;
		col_str[3][1] = 0.799;
		col_str[4][1] = 0.225;
		col_str[5][1] = 0.335;
		col_str[6][1] = 0.555;
		col_str[7][1] = 0.609; 
		col_str[8][1] = 0.367;
		col_str[9][1] = 0.195;
		col_str[10][1] = 0.0698;
		col_str[11][1] = 0.538;
		col_str[12][1] = 0.484;
		col_str[13][1] = 0.285;

		col_str[3][2] = 1.29;
		col_str[4][2] = 0.312;
		col_str[5][2] = 0.173;
		col_str[6][2] = 0.178;
		col_str[7][2] = 0.430;
		col_str[8][2] = 0.486;
		col_str[9][2] = 0.179;
		col_str[10][2] = 0.0741;
		col_str[11][2] = 0.249;
		col_str[12][2] = 0.362;
		col_str[13][2] = 0.324;

		col_str[4][3] = 0.493;
		col_str[5][3] = 0.0767;
		col_str[6][3] = 0.0348;
		col_str[7][3] = 0.223;
		col_str[8][3] = 0.401;
		col_str[9][3] = 0.126;
		col_str[10][3] = 0.0528;
		col_str[11][3] = 0.101;
		col_str[12][3] = 0.207;
		col_str[13][3] = 0.253;

		col_str[5][4] = 0.0211;
		col_str[6][4] = 0.00122; 
		col_str[7][4] = 0.0653;
		col_str[8][4] = 0.154;
		col_str[9][4] = 0.0453;
		col_str[10][4] = 0.0189;
		col_str[11][4] = 0.0265;
		col_str[12][4] = 0.0654;
		col_str[13][4] = 0.0950;

		col_str[6][5] = 0.403;
		col_str[7][5] = 0.213;
		col_str[8][5] = 0.0939;
		col_str[9][5] = 1.10;
		col_str[10][5] = 0.282;
		col_str[11][5] = 0.942;
		col_str[12][5] = 0.768;
		col_str[13][5] = 0.579;

		col_str[7][6] = 2.84; /* 10-9 */
		col_str[8][6] = 0.379; /* 11-9 */
		col_str[9][6] = 0.0876;  /* 7-9 */
		col_str[10][6] = 0.00807; /* 8-9 */
		col_str[11][6] = 1.85; /* 12-9 */
		col_str[12][6] = 0.667; /* 13-9 */
		col_str[13][6] = 0.0905; /* 14-9 */

		col_str[8][7] = 3.07; /* 11-10 */
		col_str[9][7] = 0.167;   /* 7-10 */
		col_str[10][7] = 0.0526;  /* 8-10 */
		col_str[11][7] = 0.814; /* 12-10 */
		col_str[12][7] = 0.837; /* 13-10 */
		col_str[13][7] = 0.626; /* 14-10 */

		col_str[9][8] = 0.181; /* 7-11 */
		col_str[10][8] = 0.0854; /* 8-11 */
		col_str[11][8] = 0.180; /* 12-11 */
		col_str[12][8] = 0.778; /* 13-11 */
		col_str[13][8] = 0.941; /* 14-11 */

		col_str[10][9] = 0.377; /* 8-7 */
		col_str[11][9] = 0.603; /* 12-7 */
		col_str[12][9] = 0.472; /* 13-7 */
		col_str[13][9] = 0.302; /* 14-7 */

		col_str[11][10] = 0.216; /* 12-8 */
		col_str[12][10] = 0.137; /* 13-8 */
		col_str[13][10] = 0.106; /* 14-8 */

		col_str[12][11] = 1.25;
		col_str[13][11] = 0.292;

		col_str[13][12] = 1.10;
	}

	/* bail if no ions */
	if( dense.xIonDense[ipIRON][2] <= 0. )
	{
		CoolAdd("Fe 3",0,0.);

		fe.Fe3CoolTot = 0.;   
		for( ihi=1; ihi<NLFE3; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				fe.Fe3_emiss[ihi][ilo] = 0.;
			}
		}
		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( NLFE3 <= LIMLEVELN);
		for( i=0; i < NLFE3; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}

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

	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(
		/* number of levels */
		NLFE3,
		/* the abundance of the ion, cm-3 */
		dense.xIonDense[ipIRON][2],
		/* the statistical weights */
		gfe3,
		/* the excitation energies */
		excit_wn,
		'w',
		/* the derived populations - cm-3 */
		pops,
		/* the derived departure coefficients */
		depart,
		/* the net emission rate, Aul * escp prob */
		&AulNet ,
		/* the collision strengths */
		&col_str ,
		/* A * dest prob */
		&AulDest,               
		/* pumping rate */
		&AulPump,
		/* collision rate, s-1, must defined if no collision strengths */
		&CollRate,
		/* creation vector */
		create,
		/* destruction vector */
		destroy,
		/* say atom_levelN should evaluate coll rates from cs */
		FALSE,   
		&fe.Fe3CoolTot,
		&dfe3dt,
		"Fe 3",
		&lgNegPop,
		&lgZeroPop,
		FALSE );

	/* LIMLEVELN is the dim of the PopLevels vector */
	ASSERT( NLFE3 <= LIMLEVELN );
	for( i=0; i<NLFE3; ++i)  
	{
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}

	if( lgNegPop>0 )  
	{
		fprintf( ioQQQ, " Fe3Lev14 found negative populations\n" );
		ShowMe();
		puts( "[Stop in fe3lev14]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* add cooling then its deriv */
	CoolAdd("Fe 3",0,fe.Fe3CoolTot);
	/* deriv of cooling */
	thermal.dCooldT += dfe3dt;

	/* find emission in each line */
	for( ihi=1; ihi<NLFE3; ++ihi )
	{
		for( ilo=0; ilo<ihi; ++ilo )
		{
			/* emission in these lines */
			fe.Fe3_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN;
		}
	}

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

/*Fe11Lev5 compute populations and cooling due to 5 level Fe 11 ion */
static void Fe11Lev5(void)
{
	int lgNegPop,
		lgZeroPop;
	long int i,
		j;
	static int lgFirst=TRUE;

	double dCool_dT;

	long int ihi , ilo;
	static double
		*depart,
		*pops,
		*destroy,
		*create ,
		**AulDest,
		**CollRate,
		**AulPump,
		**AulNet,
		**col_str ,
		TeUsed=-1.;

	/* statistical weights for the n levels */
	static double stat_wght[NLFE11]={5.,3.,1.,5.,1.};

	/*refer Fe11	energies	NIST version 3 Atomic Spectra Database */
	/* from smallest to largest energy (not in multiplet groupings) */

	/* energy in wavenumbers */
	static double excit_wn[NLFE11]={
		0.0 , 12667.9 , 14312. , 37743.6 , 80814.7 };

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

	if( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* allocate the 1D arrays*/
		/* create space for the 2D arrays */
		if( (depart = ((double *)MALLOC((NLFE11)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (pops = ((double *)MALLOC((NLFE11)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (destroy = ((double *)MALLOC((NLFE11)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (create = ((double *)MALLOC((NLFE11)*sizeof(double)))) == NULL )
			BadMalloc();
		/* now the 2-d arrays */
		if( (fe.Fe11_wl = ((double **)MALLOC((NLFE11)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (fe.Fe11_emiss = ((double **)MALLOC((NLFE11)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulNet = ((double **)MALLOC((NLFE11)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (col_str = ((double **)MALLOC((NLFE11)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulPump = ((double **)MALLOC((NLFE11)*sizeof(double *)))) == NULL )
			BadMalloc();  
		if( (CollRate = ((double **)MALLOC((NLFE11)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulDest = ((double **)MALLOC((NLFE11)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE11); ++i )
		{
			if( (fe.Fe11_wl[i] = ((double *)MALLOC((NLFE11)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (fe.Fe11_emiss[i] = ((double *)MALLOC((NLFE11)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulNet[i] = ((double *)MALLOC((NLFE11)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (col_str[i] = ((double *)MALLOC((NLFE11)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulPump[i] = ((double *)MALLOC((NLFE11)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE11)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulDest[i] = ((double *)MALLOC((NLFE11)*sizeof(double )))) == NULL )
				BadMalloc();
		}

		/* set some to constant values after zeroing out */
		for( i=0; i<NLFE11; ++i )
		{
			create[i] = 0.;
			destroy[i] = 0.;
			for( j=0; j<NLFE11; ++j )
			{
				AulNet[i][j] = 0.;
				col_str[i][j] = 0.;
				CollRate[i][j] = 0.;
				AulDest[i][j] = 0.;
				AulPump[i][j] = 0.;
				fe.Fe11_wl[i][j] = 0.;
				fe.Fe11_emiss[i][j] = 0.;
			}
		}
		/* calculates wavelengths of transitions */
		/* dividing by RefIndex converts the vacuum wavelength to air wavelength */
		for( ihi=1; ihi<NLFE11; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				fe.Fe11_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / RefIndex( (excit_wn[ihi]-excit_wn[ilo]) );
			}
		}

		/*>>refer	Fe11	as	Mendoza C. & Zeippen, C. J. 1983, MNRAS, 202, 981 */
		AulNet[4][0] = 1.7;
		AulNet[4][1] = 990.;
		AulNet[4][3] = 8.3;
		AulNet[3][0] = 92.3;
		AulNet[3][1] = 9.44;
		AulNet[3][2] = 1.4e-3;
		AulNet[2][0] = 9.9e-3;
		AulNet[1][0] = 43.7;
		AulNet[2][1] = 0.226;

	}

	/* bail if no ions */
	if( dense.xIonDense[ipIRON][10] <= 0. )
	{
		CoolAdd("Fe11",0,0.);

		fe.Fe11CoolTot = 0.;   
		for( ihi=1; ihi<NLFE11; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				fe.Fe11_emiss[ihi][ilo] = 0.;
			}
		}
		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( NLFE11 <= LIMLEVELN);
		for( i=0; i < NLFE11; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}

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

	/* evaluate collision strengths if Te changed by too much */
	if( fabs(phycon.te / TeUsed - 1. ) > 0.05 )
	{
		TeUsed = phycon.te;

		/* collision strengths at current temperature */
		for( ihi=1; ihi<NLFE11; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				col_str[ihi][ilo] = Fe_10_11_13_cs( 11 , ilo+1 , ihi+1 );
				ASSERT( col_str[ihi][ilo] > 0. );
			}
		}
	}

	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(
		/* number of levels */
		NLFE11,
		/* the abundance of the ion, cm-3 */
		dense.xIonDense[ipIRON][10],
		/* the statistical weights */
		stat_wght,
		/* the excitation energies */
		excit_wn,
		'w',
		/* the derived populations - cm-3 */
		pops,
		/* the derived departure coefficients */
		depart,
		/* the net emission rate, Aul * escp prob */
		&AulNet ,
		/* the collision strengths */
		&col_str ,
		/* A * dest prob */
		&AulDest,               
		/* pumping rate */
		&AulPump,
		/* collision rate, s-1, must defined if no collision strengths */
		&CollRate,
		/* creation vector */
		create,
		/* destruction vector */
		destroy,
		/* say atom_levelN should evaluate coll rates from cs */
		FALSE,   
		&fe.Fe11CoolTot,
		&dCool_dT,
		"Fe11",
		&lgNegPop,
		&lgZeroPop,
		FALSE );

	/* LIMLEVELN is the dim of the PopLevels vector */
	ASSERT( NLFE11 <= LIMLEVELN );
	for( i=0; i<NLFE11; ++i)  
	{
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}

	if( lgNegPop>0 )  
	{
		fprintf( ioQQQ, " Fe11Lev5 found negative populations\n" );
		ShowMe();
		puts( "[Stop in Fe11Lev5]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* add cooling then its deriv */
	CoolAdd("Fe11",0,fe.Fe11CoolTot);
	/* deriv of cooling */
	thermal.dCooldT += dCool_dT;

	/* find emission in each line */
	for( ihi=1; ihi<NLFE11; ++ihi )
	{
		for( ilo=0; ilo<ihi; ++ilo )
		{
			/* emission in these lines */
			fe.Fe11_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN;
		}
	}

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

/*Fe13Lev5 compute populations and cooling due to 5 level Fe 13 ion */
static void Fe13Lev5(void)
{
	int lgNegPop,
		lgZeroPop;
	long int i,
		j;
	static int lgFirst=TRUE;

	double dCool_dT;

	long int ihi , ilo;
	static double
		*depart,
		*pops,
		*destroy,
		*create ,
		**AulDest,
		**CollRate,
		**AulPump,
		**AulNet,
		**col_str ,
		TeUsed=-1.;

	/* statistical weights for the n levels */
	static double stat_wght[NLFE13]={1.,3.,5.,5.,1.};

	/*refer Fe13	energies	NIST version 3 Atomic Spectra Database */
	/* energy in wavenumbers */
	static double excit_wn[NLFE13]={
		0.0 , 9302.5 , 18561.0 , 48068. , 91508. };

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

	if( lgFirst )
	{
		/* will never do this again */
		lgFirst = FALSE;
		/* allocate the 1D arrays*/
		/* create space for the 2D arrays */
		if( (depart = ((double *)MALLOC((NLFE13)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (pops = ((double *)MALLOC((NLFE13)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (destroy = ((double *)MALLOC((NLFE13)*sizeof(double)))) == NULL )
			BadMalloc();
		if( (create = ((double *)MALLOC((NLFE13)*sizeof(double)))) == NULL )
			BadMalloc();
		/* now the 2-d arrays */
		if( (fe.Fe13_wl = ((double **)MALLOC((NLFE13)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (fe.Fe13_emiss = ((double **)MALLOC((NLFE13)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulNet = ((double **)MALLOC((NLFE13)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (col_str = ((double **)MALLOC((NLFE13)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulPump = ((double **)MALLOC((NLFE13)*sizeof(double *)))) == NULL )
			BadMalloc();  
		if( (CollRate = ((double **)MALLOC((NLFE13)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (AulDest = ((double **)MALLOC((NLFE13)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE13); ++i )
		{
			if( (fe.Fe13_wl[i] = ((double *)MALLOC((NLFE13)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (fe.Fe13_emiss[i] = ((double *)MALLOC((NLFE13)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulNet[i] = ((double *)MALLOC((NLFE13)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (col_str[i] = ((double *)MALLOC((NLFE13)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulPump[i] = ((double *)MALLOC((NLFE13)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE13)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (AulDest[i] = ((double *)MALLOC((NLFE13)*sizeof(double )))) == NULL )
				BadMalloc();
		}

		/* set some to constant values after zeroing out */
		for( i=0; i<NLFE13; ++i )
		{
			create[i] = 0.;
			destroy[i] = 0.;
			for( j=0; j<NLFE13; ++j )
			{
				AulNet[i][j] = 0.;
				col_str[i][j] = 0.;
				CollRate[i][j] = 0.;
				AulDest[i][j] = 0.;
				AulPump[i][j] = 0.;
				fe.Fe13_wl[i][j] = 0.;
				fe.Fe13_emiss[i][j] = 0.;
			}
		}
		/* calculates wavelengths of transitions */
		/* dividing by RefIndex converts the vacuum wavelength to air wavelength */
		for( ihi=1; ihi<NLFE13; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				fe.Fe13_wl[ihi][ilo] = 1e8/(excit_wn[ihi]-excit_wn[ilo]) / RefIndex( (excit_wn[ihi]-excit_wn[ilo]) );
			}
		}

		/*>>refer	Fe13	as	Huang, K.-N. 1985, At. Data Nucl. Data Tables 32, 503 */
		AulNet[4][1] = 1010.;
		AulNet[4][2] = 3.8;
		AulNet[4][3] = 8.1;
		AulNet[3][1] = 45.7;
		AulNet[3][2] = 57.6;
		AulNet[2][0] = 0.0063;
		AulNet[1][0] = 14.0;
		AulNet[2][1] = 9.87;

	}

	/* bail if no ions */
	if( dense.xIonDense[ipIRON][12] <= 0. )
	{
		CoolAdd("Fe13",0,0.);

		fe.Fe13CoolTot = 0.;   
		for( ihi=1; ihi<NLFE13; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				fe.Fe13_emiss[ihi][ilo] = 0.;
			}
		}
		/* level populations */
		/* LIMLEVELN is the dimension of the atoms vectors */
		ASSERT( NLFE13 <= LIMLEVELN);
		for( i=0; i < NLFE13; i++ )
		{
			atoms.PopLevels[i] = 0.;
			atoms.DepLTELevels[i] = 1.;
		}

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

	/* evaluate collision strengths if Te changed by too much */
	if( fabs(phycon.te / TeUsed - 1. ) > 0.05 )
	{
		TeUsed = phycon.te;

		/* collision strengths at current temperature */
		for( ihi=1; ihi<NLFE13; ++ihi )
		{
			for( ilo=0; ilo<ihi; ++ilo )
			{
				col_str[ihi][ilo] = Fe_10_11_13_cs( 13 , ilo+1 , ihi+1 );
				ASSERT( col_str[ihi][ilo] > 0. );
			}
		}
	}

	/*fprintf(ioQQQ,"DEBUG Fe13 N %.2e %.2e %.2e %.2e %.2e %.2e\n",
		col_str[1][0] , col_str[2][1] , col_str[2][0] ,
		AulNet[1][0] , AulNet[2][1] , AulNet[2][0]);*/

	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(
		/* number of levels */
		NLFE13,
		/* the abundance of the ion, cm-3 */
		dense.xIonDense[ipIRON][12],
		/* the statistical weights */
		stat_wght,
		/* the excitation energies */
		excit_wn,
		'w',
		/* the derived populations - cm-3 */
		pops,
		/* the derived departure coefficients */
		depart,
		/* the net emission rate, Aul * escp prob */
		&AulNet ,
		/* the collision strengths */
		&col_str ,
		/* A * dest prob */
		&AulDest,               
		/* pumping rate */
		&AulPump,
		/* collision rate, s-1, must defined if no collision strengths */
		&CollRate,
		/* creation vector */
		create,
		/* destruction vector */
		destroy,
		/* say atom_levelN should evaluate coll rates from cs */
		FALSE,   
		&fe.Fe13CoolTot,
		&dCool_dT,
		"Fe13",
		&lgNegPop,
		&lgZeroPop,
		FALSE );

	/* LIMLEVELN is the dim of the PopLevels vector */
	ASSERT( NLFE13 <= LIMLEVELN );
	for( i=0; i<NLFE13; ++i)  
	{
		atoms.PopLevels[i] = pops[i];
		atoms.DepLTELevels[i] = depart[i];
	}

	if( lgNegPop>0 )  
	{
		fprintf( ioQQQ, " Fe13Lev5 found negative populations\n" );
		ShowMe();
		puts( "[Stop in Fe13Lev5]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* add cooling then its deriv */
	CoolAdd("Fe13",0,fe.Fe13CoolTot);
	/* deriv of cooling */
	thermal.dCooldT += dCool_dT;

	/* find emission in each line */
	for( ihi=1; ihi<NLFE13; ++ihi )
	{
		for( ilo=0; ilo<ihi; ++ilo )
		{
			/* emission in these lines */
			fe.Fe13_emiss[ihi][ilo] = pops[ihi]*AulNet[ihi][ilo]*(excit_wn[ihi] - excit_wn[ilo])*T1CM*BOLTZMANN;
		}
	}

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

/*lint +e777 float tests equality - several places old temp saved and compared with current */

