/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*CoolIron compute iron cooling */
/*fe14cs compute collision strengths for forbidden transitions */
/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
/*Fe2Lev16 compute populations and cooling due to 16 level Fe II ion */
/*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
#include "cddefines.h"
#include "physconst.h"
#include "dense.h"
#include "fe.h"
#include "coolheavy.h"
#include "taulines.h"
#include "phycon.h"
#include "iso.h"
#include "converge.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"
#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);

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

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

void CoolIron()
{
	int lgNegPop;
	long int i, 
	  j;
	/*static long int **ipdest;*/
	static double **pump,
	  **CollRate,
	  **data, 
	  **dest, 
	  *depart,
	  *pops,
	  *destroy,
	  *create;
	double cs ,
	  cs12, cs13, cs23,
	  dFe2dt, 
	  tused,
	  cs2s2p, 
	  cs2s3p,
	  PumpLyaFeII;
	float p2, 
		rate;

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

	static int lgFirst=TRUE;

#	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 + 
		8.0e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);
	LineConvRate2CS( &TauLines[ipFe1_24m]  , rate );

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

	rate = (float)((1.2e-7 * dense.eden + 
		6.9e-10*pow((phycon.te/100.), 0.17 )*dense.xIonDense[ipHYDROGEN][0]) / dense.eden);
	TauDummy.gHi = TauLines[ipFe1_35m].gHi;
	LineConvRate2CS( &TauDummy  , rate );
	TauDummy.gHi = 0.;

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

	/* set up FeII model atom */
	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( (pump = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (dest = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (data = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();

		for( i=0; i<(NLFE2); ++i )
		{
			if( (pump[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (dest[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (data[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
		}
	}

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


	if( FeII.lgFeIION )
	{
		/* this is the large 371 level FeII atom
		 * >>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 */

		/* only redo large feii when redoing everything */
		if( dense.xIonDense[ipIRON][1] == 0. )
		{
			/* zero abundance, do nothing */
			FeII.dfcool = 0.;
			FeII.fcool = 0.;
			FeII.feheat = 0.;
			/* now zero out intensities of all feii lines */
			FeIIIntenZero();
		}
		else
		{
			/* these are used to only evaluated feii one time per temperature, zone, and abundance */
			static float TUsed = 0. , 
				AbunUsed = 0.;
			/* remember which zone last evaluated with */
			static long int nZUsed=-1,
				/* make sure at least two calls per zone */
				nCall=0;

			/* 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.fcool/thermal.ctot)> 0.002 &&  
				fabs(dense.xIonDense[ipIRON][1]-AbunUsed)/SDIV(AbunUsed)> 0.002 ) ||
				(phycon.te != TUsed &&  fabs(FeII.fcool/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.fcool/thermal.ctot> 0.001 )
					{
						fprintf( ioQQQ, 
							"temp or cooling changed, new are %g %g nCall %li ", 
							phycon.te, FeII.fcool, nCall );
					}
					else
					{
						fprintf(ioQQQ, "????");
					}
				}

				/* 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.fcool/thermal.ctot,nCall);
				}

				FeIILevelPops();

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

			}
			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.fcool,thermal.ctot);
				}
				ratio = dense.xIonDense[ipIRON][1]/AbunUsed;
				FeII.fcool *= ratio;
				FeII.dfcool *= ratio;
				FeII.feheat *= 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 heating/cooling */
			/* total cooling */
			CoolHeavy.Fecool = FeII.fcool;

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

			CoolAdd("Fe 2",0,MAX2(0.,FeII.fcool));

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

			if( trace.lgTrace && trace.lgCoolTr )
			{
				fprintf( ioQQQ, " Large FeII returns te, cooling, dc=%11.3e%11.3e%11.3e\n", 
				  phycon.te, CoolHeavy.Fecool, FeII.dfcool );
			}
		}
	}
	else
	{
		FeII.feheat = 0.;

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

		/*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.;
				dest[j][i] = 0.;
				pump[j][i] = 0.;
			}
		}

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

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

		/*ipdest[0][5] = TauLines[ipTFe16].ipCont;*/

		data[1][2] = (TauLines[ipTr48].Pesc + TauLines[ipTr48].Pelec_esc)*TauLines[ipTr48].Aul;
		dest[1][2] = TauLines[ipTr48].Pdest*TauLines[ipTr48].Aul;
		pump[1][2] = TauLines[ipTr48].pump;
		/*ipdest[1][2] = TauLines[ipTr48].ipCont;*/

		data[1][5] = (TauLines[ipTFe26].Pesc + TauLines[ipTFe26].Pelec_esc)*TauLines[ipTFe26].Aul;
		dest[1][5] = TauLines[ipTFe26].Pdest*TauLines[ipTFe26].Aul;
		pump[1][5] = TauLines[ipTFe26].pump;
		/*ipdest[1][5] = TauLines[ipTFe26].ipCont;*/

		data[2][3] = (TauLines[ipTFe34].Pesc + TauLines[ipTFe34].Pelec_esc)*TauLines[ipTFe34].Aul;
		dest[2][3] = TauLines[ipTFe34].Pdest*TauLines[ipTFe34].Aul;
		pump[2][3] = TauLines[ipTFe34].pump;
		/*ipdest[2][3] = TauLines[ipTFe34].ipCont;*/

		data[2][4] = (TauLines[ipTFe35].Pesc + TauLines[ipTFe35].Pelec_esc)*TauLines[ipTFe35].Aul;
		dest[2][4] = TauLines[ipTFe35].Pdest*TauLines[ipTFe35].Aul;
		pump[2][4] = TauLines[ipTFe35].pump;
		/*ipdest[2][4] = TauLines[ipTFe35].ipCont;*/

		data[3][5] = (TauLines[ipTFe46].Pesc + TauLines[ipTFe46].Pelec_esc)*TauLines[ipTFe46].Aul;
		dest[3][5] = TauLines[ipTFe46].Pdest*TauLines[ipTFe46].Aul;
		pump[3][5] = TauLines[ipTFe46].pump;
		/*ipdest[3][5] = TauLines[ipTFe46].ipCont;*/

		data[4][5] = (TauLines[ipTFe56].Pesc + TauLines[ipTFe56].Pelec_esc)*TauLines[ipTFe56].Aul;
		dest[4][5] = TauLines[ipTFe56].Pdest*TauLines[ipTFe56].Aul;
		pump[4][5] = TauLines[ipTFe56].pump;
		/*ipdest[4][5] = TauLines[ipTFe56].ipCont;*/

		/* these are collision strengths */
		data[1][0] = 1.;
		data[2][0] = 12.;
		data[3][0] = 1.;
		data[4][0] = 1.;
		data[5][0] = 12.;
		data[2][1] = 6.;
		data[3][1] = 1.;
		data[4][1] = 1.;
		data[5][1] = 12.;
		data[3][2] = 6.;
		data[4][2] = 12.;
		data[5][2] = 1.;
		data[4][3] = 1.;
		data[5][3] = 12.;
		data[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,
		  pops,
		  depart,
		  &data,
		  &dest,
		  &pump,
		  &CollRate,
		  create,
		  destroy,
		  FALSE,/* say atom_levelN should evaluate coll rates from cs */
		  /*&&ipdest,*/
		  &CoolHeavy.Fecool,
		  &dFe2dt,
		  "FeII",
		  &lgNegPop,
		  FALSE );
		/* 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];
		}

		/* lgNegPop positive if negative pops occured, negative if too cold */
		if( lgNegPop>0 )
		{
			fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population for FeII.\n");
		}

		CoolAdd("Fe 2",0,MAX2(0.,CoolHeavy.Fecool));
		thermal.heating[25][27] = MAX2(0.,-CoolHeavy.Fecool);
		thermal.dCooldT += dFe2dt;

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

		TauLines[ipTr48].PopLo = pops[1];
		TauLines[ipTr48].PopOpc = (pops[1] - pops[2]);
		/*CoolHeavy.rmt48 = pops[2]*data[1][2]*TauLines[ipTr48].EnergyErg;*/
		TauLines[ipTr48].phots = pops[2]*data[1][2];
		TauLines[ipTr48].xIntensity = 
			TauLines[ipTr48].phots*TauLines[ipTr48].EnergyErg;

		CoolHeavy.for7 = pops[1]*data[0][1]*4.65e-12;

		TauLines[ipTFe16].PopLo = pops[0];
		TauLines[ipTFe16].PopOpc = (pops[0] - pops[5]);
		TauLines[ipTFe16].phots = pops[5]*data[0][5];
		TauLines[ipTFe16].xIntensity = 
			TauLines[ipTFe16].phots*TauLines[ipTFe16].EnergyErg;
		/*CoolHeavy.cFe16 = TauLines[ipTFe16].phots*TauLines[ipTFe16].EnergyErg;*/

		TauLines[ipTFe26].PopLo = (pops[1]);
		TauLines[ipTFe26].PopOpc = (pops[1] - pops[5]);
		TauLines[ipTFe26].phots = pops[5]*data[1][5];
		TauLines[ipTFe26].xIntensity = 
			TauLines[ipTFe26].phots*TauLines[ipTFe26].EnergyErg;
		/*CoolHeavy.cFe26 = TauLines[ipTFe26].phots*TauLines[ipTFe26].EnergyErg;*/

		TauLines[ipTFe34].PopLo = pops[2];
		TauLines[ipTFe34].PopOpc = (pops[2] - pops[3]);
		TauLines[ipTFe34].phots = pops[3]*data[2][3];
		TauLines[ipTFe34].xIntensity = 
			TauLines[ipTFe34].phots*TauLines[ipTFe34].EnergyErg;
		/*CoolHeavy.cFe34 = TauLines[ipTFe34].phots*TauLines[ipTFe34].EnergyErg;*/

		TauLines[ipTFe35].PopLo = pops[2];
		TauLines[ipTFe35].PopOpc = (pops[2] - pops[4]);
		TauLines[ipTFe35].phots = pops[4]*data[2][4];
		TauLines[ipTFe35].xIntensity = 
			TauLines[ipTFe35].phots*TauLines[ipTFe35].EnergyErg;
		/*CoolHeavy.cFe35 = TauLines[ipTFe35].phots*TauLines[ipTFe35].EnergyErg;*/

		TauLines[ipTFe46].PopLo = pops[3];
		TauLines[ipTFe46].PopOpc = (pops[3] - pops[5]);
		TauLines[ipTFe46].phots = pops[5]*data[3][5];
		TauLines[ipTFe46].xIntensity = 
			TauLines[ipTFe46].phots*TauLines[ipTFe46].EnergyErg;
		/*CoolHeavy.cFe46 = 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]*data[4][5];
		TauLines[ipTFe56].xIntensity = 
			TauLines[ipTFe56].phots*TauLines[ipTFe56].EnergyErg;
		/*CoolHeavy.cFe56 = 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]);

		/* do 16 level ground + first two terms */
		Fe2Lev16();
	}

	/* lump 3p and 3f together; cs=
	 * >>refer	fe3	vs	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 */
	CoolHeavy.c5270 = atom_pop2(5.53,25.,30.,0.435,2.948e4,dense.xIonDense[ipIRON][2])*
	  3.78e-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 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 97 jan 06 added fir lines, cs from 
	 * >>refer	fe7	as	Keenan, F.P., & Norrington, P.H. 1987 A&A, 181, 370
	 * [FeVII] 9.51, 7.81 micron lines */
	if( phycon.te < 1.5e4 )
	{
		cs = MIN2(1.49,374.8/(phycon.sqrte*phycon.te10));
		PutCS(cs,&TauLines[ipFe0795]);
		cs = MIN2(1.57,188.5/(phycon.sqrte*phycon.te02));
		PutCS(cs,&TauLines[ipFe0778]);
		cs = MIN2(0.389,51.32/(phycon.sqrte*phycon.te03));
		PutCS(cs,&TauDummy);
	}
	else
	{
		cs = MIN2(1.17,7.929/(phycon.te20/phycon.te001));
		PutCS(cs,&TauLines[ipFe0795]);
		cs = MIN2(1.27,5.638/(phycon.te10*phycon.te05*phycon.te005));
		PutCS(cs,&TauLines[ipFe0778]);
		cs = MIN2(0.314,1.491/(phycon.te10*phycon.te05*phycon.te01*
		  phycon.te001*phycon.te001));
		PutCS(cs,&TauDummy);
	}
	atom_level3(&TauLines[ipFe0795],&TauLines[ipFe0778],&TauDummy);

	/* FE VII fe07 fe 07 fe 7 fe7 6087, A from 
	 * >>refer	fe7	as	Nussbaumer, H., & Storey, P.J. 1982, A&A, 113, 21
	 * cs from 
	 * >>refer	fe7	cs	Keenan, F.P., & Norrington, P.H. 1987 A&A, 181, 370 */
	tused = MAX2(4.0,phycon.alogte);
	tused = MIN2(5.1,tused);
	cs = 2.2684 - tused*(0.60893 - 0.052004*tused);
	/* >>chng 97 jan 27, should have been total collision strength of multiplet
	 * thanks for Tino Oliva for finding this problem */
	cs *= 2.96;
	CoolHeavy.c6087 = atom_pop2(cs,21.,5.,0.94,2.32e4,dense.xIonDense[ipIRON][6])*
	  3.27e-12*0.577/0.94;
	CoolHeavy.Fe5722 = CoolHeavy.c6087*0.629*6087./5722.;
	thermal.dCooldT += CoolHeavy.c6087*2.37e4*thermal.tsq1;
	CoolAdd("Fe 7",6087,CoolHeavy.c6087);
	CoolAdd("Fe 7",5722,CoolHeavy.Fe5722);

	/* 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 
	 * >>refer	fe10	as	Mason, H. 1975, MNRAS 170, 651
	 * cs 
	 * >>refer	fe10	cs	Mohan, M., Hibbert, A., & Kingston, A.E. 1994, ApJ, 434, 389
	 * not used newer cs from 
	 * >>refer	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	cs	fe10	Tayal, S.S., 2000, ApJ, 544, 575-580 */
	cs = 172.9 * phycon.te30 * phycon.te05 * phycon.te02 *  phycon.te005;
	/* above data ends at 40,000K, do not allow cs to grow larger than largest
	 * tabulated value */
	cs = MIN2( cs, 3.251 );
	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;
	CoolAdd("Fe11",1467,CoolHeavy.c1467);
	CoolHeavy.c2649 = p2*91.0*7.512e-12;
	CoolAdd("Fe11",2469,CoolHeavy.c2649);

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

	/* 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 */
	cs = 3.1;
	CoolHeavy.c5303 = atom_pop2(cs,2.,4.,60.3,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	cs	fe19	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	cs	fe20	Butler, K., & Zeippen, C.J., 2001, A&A, 372, 1078 */
	/*>>refer	as	fe20	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);

	/*>>chng 00 oct 26, added these fe 21 lines, removed old CoolHeavy.fs846, the lowest */

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

	/* 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 strengths for forbidden transitions 
 * 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;
	long int i, 
	  j, 
	  _r;
	static int lgFirst=TRUE;

	double dfe4dt;

	/*static long int **ipdest; */
	static double **data, 
	  **dest, 
	  depart[NLFE4],
	  pops[NLFE4], 
	  destroy[NLFE4], 
	  create[NLFE4] ,
	  **CollRate,
	  **pump;

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

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

	static double ex[NLFE4]={0.,46395.8,46464.,46475.6,46483.,50725.,
	  50838.,50945.,55796.,55966.,56021.,56025.};

	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( (pump = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (dest = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (data = ((double **)MALLOC((NLFE4)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE4); ++i )
		{
			if( (pump[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (dest[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (data[i] = ((double *)MALLOC((NLFE4)*sizeof(double )))) == NULL )
				BadMalloc();
		}
	}

	/* bail if no ions */
	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 ex/ 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;
			dest[j][i] = 0.;
			pump[j][i] = 0.;
		}
		data[i][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];
			data[j][i] = Fe4CS[j][i];
		}
	}

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

	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(NLFE4,
		dense.xIonDense[ipIRON][3],
		gfe4,
		ex,
		pops,
		depart,
		&data,
		&dest,
		&pump,
		&CollRate,
		create,
		destroy,
		FALSE,/* say atom_levelN should evaluate coll rates from cs */
		/*&ipdest,*/
		&fe.Fe4CoolTot,
		&dfe4dt,
		"FeIV",
		&lgNegPop,
		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] + pops[4]*Fe4A[0][4])*
	  6.42e-12;
	fe.fe42836 = pops[5]*Fe4A[0][5]*7.02e-12;
	fe.fe42829 = pops[6]*Fe4A[0][6]*7.03e-12;
	fe.fe42567 = (pops[10]*Fe4A[0][10] + pops[11]*Fe4A[0][11])*
	  7.752e-12;
	fe.fe41207 = pops[11]*Fe4A[6][11]*(ex[11] - ex[6])*BOLTZMANN;
	fe.fe41206 = pops[11]*Fe4A[5][11]*(ex[11] - ex[5])*BOLTZMANN;
	fe.fe41106 = pops[10]*Fe4A[5][10]*(ex[10] - ex[5])*BOLTZMANN;
	fe.fe41007 = pops[9]*Fe4A[6][9]*(ex[9] - ex[6])*BOLTZMANN;
	fe.fe41008 = pops[9]*Fe4A[7][9]*(ex[9] - ex[7])*BOLTZMANN;
	fe.fe40906 = pops[8]*Fe4A[5][8]*(ex[8] - ex[5])*BOLTZMANN;

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

#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;
	long int i, 
	  j, 
	  nUsed, 
	  _r;
	/*static long int **ipdest; */

	double dfe2dt;
	static double **data, 
	  **dest, 
	  *pops, 
	  *create, 
	  *destroy, 
	  *depart,
	  **CollRate,
	  **pump;
	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.};
	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 data(up,lo)
	 * A's are saved into arrays as data(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( (pump = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (dest = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (data = ((double **)MALLOC((NLFE2)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE2); ++i )
		{
			if( (pump[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (dest[i] = ((double *)MALLOC((NLFE2)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (data[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.Fe2L16Tot = 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;
			dest[j][i] = 0.;
			pump[j][i] = 0.;
		}
		data[i][i] = 0.;
	}

	/* fill in einstein As and collision strengths */
	for( i=0; i < NLFE2; i++ )
	{
		for( j=i + 1; j < NLFE2; j++ )
		{
			data[i][j] = Fe2A[i][j];
			data[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,
		pops,
		depart,
		&data,
		&dest,
		&pump,
		&CollRate,
		create,
		destroy,
		FALSE,/* say atom_levelN should evaluate coll rates from cs */
		/*&ipdest,*/
		&FeII.Fe2L16Tot,
		&dfe2dt,
		"FeII",
		&lgNegPop,
		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);
	}

	CoolAdd("Fe 2",2,FeII.Fe2L16Tot);

	thermal.dCooldT += dfe2dt;

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

#if 0
/*Fe7Lev8 compute populations and cooling due to 8 level Fe VII ion */
static void Fe7Lev8(void)
{
#	define	NLFE7	12
	int lgNegPop;
	long int i, 
	  j, 
	  _r;
	static int lgFirst=TRUE;

	double dfe4dt;

	/*static long int **ipdest; */
	static double **data, 
	  **dest, 
	  depart[NLFE7],
	  pops[NLFE7], 
	  destroy[NLFE7], 
	  create[NLFE7] ,
	  **CollRate,
	  **pump;

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

	static double gfe4[NLFE7]={5.,7.,9.,5.,1.,3.,5.,9.};

	static double ex[NLFE7]={0.,1047.,2327.,17475.,20037.,20428.,21275.,28915.};
#if 0
	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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; j++ )
		{
			Fe4A[10][j-1] = _itmp10[_r++];
			}
		}
		for( j=1; j <= NLFE7; j++ )
		{
			Fe4A[11][j-1] = 0.;
			}
		for( i=1; i <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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 <= NLFE7; 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( (pump = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (dest = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (data = ((double **)MALLOC((NLFE7)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(NLFE7); ++i )
		{
			if( (pump[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (dest[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (data[i] = ((double *)MALLOC((NLFE7)*sizeof(double )))) == NULL )
				BadMalloc();
		}
	}

	/* bail if no ions */
	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( NLFE7 <= LIMLEVELN);
		for( i=0; i < NLFE7; 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 ex/ 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 < NLFE7; i++ )
	{
		create[i] = 0.;
		destroy[i] = 0.;
		for( j=0; j < NLFE7; j++ )
		{
			data[j][i] = 1e33;
			dest[j][i] = 0.;
			pump[j][i] = 0.;
		}
		data[i][i] = 0.;
	}

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

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

	/* lgNegPop positive if negative pops occured, negative if too cold */
	atom_levelN(NLFE7,
		dense.xIonDense[ipIRON][3],
		gfe4,
		ex,
		pops,
		depart,
		&data,
		&dest,
		&pump,
		&CollRate,
		create,
		destroy,
		FALSE,/* say atom_levelN should evaluate coll rates from cs */
		/*&ipdest,*/
		&fe.Fe4CoolTot,
		&dfe4dt,
		"FeIV",
		&lgNegPop,
		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( 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] + pops[4]*Fe4A[0][4])*
	  6.42e-12;
	fe.fe42836 = pops[5]*Fe4A[0][5]*7.02e-12;
	fe.fe42829 = pops[6]*Fe4A[0][6]*7.03e-12;
	fe.fe42567 = (pops[10]*Fe4A[0][10] + pops[11]*Fe4A[0][11])*
	  7.752e-12;
	fe.fe41207 = pops[11]*Fe4A[6][11]*(ex[11] - ex[6])*BOLTZMANN;
	fe.fe41206 = pops[11]*Fe4A[5][11]*(ex[11] - ex[5])*BOLTZMANN;
	fe.fe41106 = pops[10]*Fe4A[5][10]*(ex[10] - ex[5])*BOLTZMANN;
	fe.fe41007 = pops[9]*Fe4A[6][9]*(ex[9] - ex[6])*BOLTZMANN;
	fe.fe41008 = pops[9]*Fe4A[7][9]*(ex[9] - ex[7])*BOLTZMANN;
	fe.fe40906 = pops[8]*Fe4A[5][8]*(ex[8] - ex[5])*BOLTZMANN;
#endif

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



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

