/*PrtZone print out individual zone results */
#include <string.h>
#include "cddefines.h"
/* following is quantum heating punch io that will need to be removed -
 * now it punches at all times */
#include "physconst.h"
#include "nhe3lvl.h"
#include "nhe1lvl.h"
#include "grainvar.h"
#include "nh.h"
#include "pressure.h"
#include "wind.h"
#include "converge.h"
#include "pbn.h"
#include "trace.h"
#include "dcneg.h"
#include "called.h"
#include "pheat.h"
#include "heating.h"
#include "he1bn.h"
#include "he3n.h"
#include "he3lines.h"
#include "secondaries.h"
#include "opac.h"
#include "he3bn.h"
#include "phe1lv.h"
#include "photn2.h"
#include "colden.h"
#include "filfac.h"
#include "hmi.h"
#include "rfield.h"
#include "heat.h"
#include "compton.h"
#include "hehp.h"
#include "h.h"
#include "he.h"
#include "radius.h"
#include "zonecnt.h"
#include "phycon.h"
#include "hmidep.h"
#include "hevmolec.h"
#include "elmton.h"
#include "lyaext.h"
#include "fluct.h"
#include "ionfracs.h"
#include "hydrogenic.h"
#include "abuntabllg.h"
#include "ionrange.h"
#include "printe82.h"
#include "elementnames.h"
#include "prtlinepres.h"
#include "pop371.h"
#include "prthydro.h"
#include "prtzone.h"

void PrtZone(void)
{
	char chField7[8];
	char chLet, 
	  chQHMark;
	long int i, 
	  ishift, 
	  nelem ,
	  nd, 
	  nq;
	double cdif, 
	  coninc, 
	  dmidle, 
	  fac, 
	  hatmic, 
	  qp1[NDUST], 
	  qt1[NDUST], 
	  rmidle;

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

	if( dcneg.lgdCNeg )
	{
		chLet = 'u';
	}
	else
	{
		chLet = ' ';
	}

	/* middle of zone for printing */
	rmidle = radius.Radius - radius.drad*0.5*radius.dRadSign;
	dmidle = radius.depth - radius.drad*0.5;

	/* option to print single line when quiet but tracing convergence
	 * with "trace convergence" command */
	if( called.lgTalk || trace.lgTrConvg )
	{
		/* print either ####123 or ###1234 */
		if( ZoneCnt.nzone <= 999 )
		{
			sprintf( chField7, "####%3ld", ZoneCnt.nzone );
		}
		else
		{
			sprintf( chField7, "###%4ld", ZoneCnt.nzone );
		}

		fprintf(ioQQQ, " %7.7s %cTe:",chField7, chLet);
		PrintE93(ioQQQ,phycon.te);
		fprintf(ioQQQ," Hden:");
		PrintE93(ioQQQ,phycon.hden);
		fprintf(ioQQQ," Ne:");
		PrintE93(ioQQQ,phycon.eden);
		fprintf(ioQQQ," R:");
		PrintE93(ioQQQ,rmidle );
		fprintf(ioQQQ," R-R0:");
		PrintE93(ioQQQ,dmidle);
		fprintf(ioQQQ," dR:");
		PrintE93(ioQQQ,radius.drad);
		fprintf(ioQQQ," NTR:%3ld Htot:",conv.nPres2Ioniz);
		PrintE93(ioQQQ,heat.htot);
		fprintf(ioQQQ," T912:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",opac.tauabs[0][nh.ipHn[0][IP1S]-1] ));
		fprintf(ioQQQ,"###\n");
		  
		if( trace.lgTrConvg )
		{
			fprintf( ioQQQ, " H:%10.2e%10.2e He:%10.2e%10.2e%10.2e trip%10.2e\n", 
			  h.hi/(h.hi + h.hii), 
			  h.hii/(h.hi + h.hii), 
			  xIonFracs[1][ipHE]/xIonFracs[0][ipHE], 
			  xIonFracs[2][ipHE]/xIonFracs[0][ipHE], 
			  xIonFracs[3][ipHE]/xIonFracs[0][ipHE], 
			  he.hei3/xIonFracs[0][ipHE] );
		}
	}

	/* now return if not talking */
	if( !called.lgTalk || trace.lgTrConvg )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->PrtZone()\n", debug_fp );
#		endif
		return;
	}

	/* lgDenFluc set to true in zero, only false when variable abundances
	 * lgAbTaON set true when element table used */
	if( !fluct.lgDenFluc || AbunTablLG.lgAbTaON )
	{
		fprintf( ioQQQ, " Abun:" );
		for( i=1; i <= LIMELM; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%8.1e", xIonFracs[0][i-1] ));
		}
		fprintf( ioQQQ, "\n" );
	}

	/*-------------------------------------------------
	 * print wind parameters if windy model */
	fac = wind.windv*phycon.hden*radius.r1r0sq;
	if( wind.windv != 0. )
	{
		if( wind.AccelTot == 0. )
		{
			fac = 1.;
		}
		else
		{
			fac = wind.AccelTot;
		}
		/*fprintf( ioQQQ, " WIND; V:%9.3ekm/s G:%10.3e Accel:%10.3e Fr(cont):%6.3f Fr(line):%6.3f Fr(dP):%6.3f\n", 
		  wind.windv/1e5, -wind.agrav, wind.AccelTot, wind.AccelCont/
		  fac, wind.AccelLine/fac, wind.AccelPres/fac );*/
		fprintf( ioQQQ, " WIND; V:");
		fprintf(ioQQQ,PrintEfmt("%10.3e",wind.windv/1e5));
		fprintf( ioQQQ, "km/s G:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",-wind.agrav));
		fprintf( ioQQQ, " Accel:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",wind.AccelTot));
		fprintf( ioQQQ, " Fr(cont):%6.3f Fr(line):%6.3f Fr(dP):%6.3f\n", 
		  wind.AccelCont/ fac, wind.AccelLine/fac, wind.AccelPres/fac );
	}

	/* print line with radiation pressure if significant */
	if( wind.windv != 0. || pressure.pbeta > .05 )
	{
		PrtLinePres();
	}

	/*---------------------------------------------------- */
	hatmic = (h.hi + h.hii)/(h.hi + h.hii + hmi.hminus + 2.*(hmi.htwo + 
	  hmi.h2plus) + 3*hmi.h3plus);

	fprintf( ioQQQ, " Hydrogen     ");
	fprintf(ioQQQ,PrintEfmt("%9.2e",h.hi/(h.hi + h.hii)));
	fprintf(ioQQQ,PrintEfmt("%9.2e",h.hii/(h.hi + h.hii)));
	fprintf( ioQQQ, " H+o/Hden");
	fprintf(ioQQQ,PrintEfmt("%9.2e",hatmic ));
	fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.hminus/phycon.hden ));
	fprintf( ioQQQ, " H-    H2");
	fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.htwo/phycon.hden));
	fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.h2plus/phycon.hden));
	fprintf( ioQQQ, " H2+ HeH+");
	fprintf(ioQQQ,PrintEfmt("%9.2e",hehpCom.hehp/phycon.hden));
	fprintf( ioQQQ, " Ho+ ColD");
	fprintf(ioQQQ,PrintEfmt("%9.2e",coldenCom.colden[IPCHI-1]));
	fprintf(ioQQQ,PrintEfmt("%9.2e",coldenCom.colden[IPCHII-1]));
	fprintf( ioQQQ, "\n");

	/* print departure coef if desired */
	if( pbn.lgPrtBN )
	{
		fprintf( ioQQQ, " Hydrogen     " );
		fprintf(ioQQQ,PrintEfmt("%9.2e",  hydro.hbn[0][IP1S]));
		fprintf(ioQQQ,PrintEfmt("%9.2e",  1.));
		fprintf( ioQQQ, " H+o/Hden");
		fprintf(ioQQQ,PrintEfmt("%9.2e", (h.hi + h.hii)/phycon.hden));
		fprintf(ioQQQ,PrintEfmt("%9.2e", hmidepCom.hmidep));
		fprintf( ioQQQ, " H-    H2");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hmidepCom.h2dep));
		fprintf( ioQQQ, "      H2+");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hmidepCom.h2pdep));
		fprintf( ioQQQ, "      H3+");
		fprintf(ioQQQ,PrintEfmt("%9.2e",hmidepCom.h3pdep));
		fprintf( ioQQQ, "\n" );
	}

	if( pheat.lgPrintHeating )
	{
		fprintf( ioQQQ, "              ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[0][0]/heat.htot));
		fprintf( ioQQQ,"         ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[15][0]/heat.htot));
		fprintf( ioQQQ,"         ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[16][0]/heat.htot));
		fprintf( ioQQQ,"\n");
	}

	/* temperature correspoding to radiation fields */
	coninc = 0.;
	cdif = 0.;
	for( i=0; i < rfield.nflux; i++ )
	{
		coninc += rfield.flux[i]*(rfield.anu[i]*EN1RYD);
		cdif += (rfield.ConOutNoInter[i] + rfield.outlin[i] + 
			rfield.outcon[i])* (rfield.anu[i]*EN1RYD);
	}
	coninc = pow(coninc/2.997925e10/7.56464e-15,0.25);
	cdif = pow(cdif/2.997925e10/7.56464e-15,0.25);


	/* print info on hydrogen level populations */
	PrtHydro(0);

	if( pheat.lgPrintHeating )
	{
		fprintf( ioQQQ, "              ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[1][0]/heat.htot ));
		fprintf( ioQQQ, "         ");
		fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
		fprintf( ioQQQ, " BoundCom");
		fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[18][0]/ heat.htot));
		fprintf( ioQQQ, "   Extra:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",HeatingCom.heating[20][0]/heat.htot));
		fprintf( ioQQQ, "   Pairs:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[21][0]/ heat.htot ));
		fprintf( ioQQQ,"  H-lines\n");
	}

	/* Helium */
	if( elmton.lgElmtOn[1] )
	{
		fprintf( ioQQQ, " Helium       " );
		for( i=1; i < 4; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", xIonFracs[i][ipHE]/xIonFracs[0][ipHE]) );
		}
		/*fprintf( ioQQQ, "%10.3e%10.3e He I 2S3 %10.3e Comp H, C%10.3e\n", 
		  Compton.cmheat, Compton.cmcool*phycon.te, filfac.FillFac, 
		  photn2Com.photn2 );*/
		fprintf( ioQQQ, " He I2SP3");
		fprintf(ioQQQ,PrintEfmt("%9.2e", he.hei3/xIonFracs[0][ipHE] ));
		fprintf(ioQQQ,PrintEfmt("%9.2e", he3lines.p2p/xIonFracs[0][ipHE] ));
		fprintf( ioQQQ, " Comp H,C");
		fprintf(ioQQQ,PrintEfmt("%9.2e", Compton.cmheat ));
		fprintf(ioQQQ,PrintEfmt("%9.2e",  Compton.cmcool*phycon.te));
		fprintf( ioQQQ , " Fill Fac");
		fprintf(ioQQQ,PrintEfmt("%9.2e", filfac.FillFac));
		fprintf( ioQQQ , " Gam1/tot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", photn2Com.photn2));
		fprintf( ioQQQ, "\n");

		/* option to print departure coef */
		if( pbn.lgPrtBN )
		{
			fprintf( ioQQQ, " Helium       " );
			fprintf(ioQQQ,PrintEfmt("%9.2e", he1bnCOM.he1bn[0]));
			fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.hbn[1][IP1S]));
			fprintf(ioQQQ,PrintEfmt("%9.2e", 1.));
			fprintf( ioQQQ, " He I 2S3");
			fprintf(ioQQQ,PrintEfmt("%9.2e", he3bnCom.he3bn[0] ));

			/*fprintf( ioQQQ, "%10.3e%10.3e He I 2S3 %10.3e Comp H, C%10.3e\n", 
			  Compton.cmheat, Compton.cmcool*phycon.te, filfac.FillFac, 
			  photn2Com.photn2 );*/
			fprintf( ioQQQ, " Comp H,C");
			fprintf(ioQQQ,PrintEfmt("%9.2e", Compton.cmheat ));
			fprintf(ioQQQ,PrintEfmt("%9.2e", Compton.cmcool*phycon.te ));
			fprintf( ioQQQ , " Fill Fac");
			fprintf(ioQQQ,PrintEfmt("%9.2e", filfac.FillFac ));
			fprintf( ioQQQ , " Gam1/tot");
			fprintf(ioQQQ,PrintEfmt("%9.2e", photn2Com.photn2));
			fprintf( ioQQQ, "\n");
		}

		/* print heating from He (and others) if desired
		 * entry "lines" is induced line heating
		 * 1,12 ffheat:  2,3 he triplets, 1,20 compton */
		if( pheat.lgPrintHeating )
		{
			/*fprintf( ioQQQ, "            %10.3e%10.3e    Lines:%10.2e%10.2e  Compton:%10.3e FF Heatig%10.3e\n", 
			  HeatingCom.heating[0][1]/heat.htot, HeatingCom.heating[1][1]/
			  heat.htot, HeatingCom.heating[22][0]/heat.htot, HeatingCom.heating[2][1]/
			  heat.htot, HeatingCom.heating[19][0]/heat.htot, HeatingCom.heating[11][0]/
			  heat.htot );*/
			fprintf( ioQQQ, "              ");
			fprintf(ioQQQ,PrintEfmt("%9.2e",HeatingCom.heating[0][1]/heat.htot));
			fprintf(ioQQQ,PrintEfmt("%9.2e",HeatingCom.heating[1][1]/heat.htot));
			fprintf( ioQQQ, "   Lines:");
			fprintf(ioQQQ,PrintEfmt("%9.2e",HeatingCom.heating[22][0]/heat.htot));
			fprintf(ioQQQ,PrintEfmt("%9.2e",HeatingCom.heating[2][1]/heat.htot));
			fprintf( ioQQQ, " Compton:");
			fprintf(ioQQQ,PrintEfmt("%9.2e",HeatingCom.heating[19][0]/heat.htot));
			fprintf( ioQQQ, " FFHeatig");
			fprintf(ioQQQ,PrintEfmt("%9.2e",HeatingCom.heating[11][0]/heat.htot));
			fprintf( ioQQQ, "\n");
		}

		/* helium singlets and triplets */
		fac = xIonFracs[2][ipHE]/xIonFracs[0][ipHE];
		fprintf( ioQQQ, " He singlet   " );
		for( i=0; i < 6; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", phe1lv.he1n[i]*fac ));
		}
		fprintf( ioQQQ, " He tripl" );
		for( i=0; i < NHE3LVL; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", he3nCom.he3n[i]*fac ));
		}
		fprintf( ioQQQ, "\n" );

		/* option to print departure coef */
		if( pbn.lgPrtBN )
		{
			fprintf( ioQQQ, " He singlet   " );
			for( i=0; i < (NHE1LVL + 1); i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", he1bnCOM.he1bn[i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " He triplet   " );
			for(i=0; i < NHE3LVL; i++)
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", he3bnCom.he3bn[i] ));
			}
			fprintf( ioQQQ, "\n" );
		}

		/* print helium ion level populations and departure coef */
		PrtHydro(1);
	}


	fprintf( ioQQQ , "               Texc(La)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  LyaExT.TexcLya ));
	fprintf( ioQQQ , " T(contn)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  coninc ));
	fprintf( ioQQQ , " T(diffs)");
	fprintf(ioQQQ,PrintEfmt("%9.2e",  cdif ));
	/* print the radiation to gas pressure */
	fprintf( ioQQQ , " PRad/Gas");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.pbeta ));
	fprintf( ioQQQ, "\n" );

	if( GrainVar.lgDustOn )
	{
		/*  Change things so the quantum heated dust species are marked with an
		 *  asterisk just after the name (K Volk) */
		for( nd=0; nd < NDUST; nd++ )
		{
			if( GrainVar.lgDustOn1[nd] )
			{
				if( GrainVar.lgQHeat[nd] )
				{
					chQHMark = '*';
				}
				else
				{
					chQHMark = ' ';
				}
				/*   (K Volk) added QHMARK here and in the write statement */
				fprintf( ioQQQ, "%10.10s%c    DustTemp",chGrainVar.chDstLab[nd], chQHMark);
				fprintf(ioQQQ,PrintEfmt("%9.2e",GrainVar.tedust[nd] ));
				fprintf( ioQQQ, " Pot Volt");
				fprintf(ioQQQ,PrintEfmt("%9.2e", GrainVar.dstpot[nd]*EVRYD));
				fprintf( ioQQQ, " Nelectrn");
				fprintf(ioQQQ,PrintEfmt("%9.2e", GrainVar.dstq[nd]));
				fprintf( ioQQQ, " drft cm ");
				fprintf(ioQQQ,PrintEfmt("%9.2e", GrainVar.DustDftVel[nd]));
				fprintf( ioQQQ, " Heating:");
				fprintf(ioQQQ,PrintEfmt("%9.2e", GrainVar.dheat[nd]));
				fprintf( ioQQQ, " Frac tot");
				fprintf(ioQQQ,PrintEfmt("%9.2e", GrainVar.dheat[nd]/ heat.htot));
				fprintf( ioQQQ, "\n" );
			}
		}
		/*  (K Volk) added QHEAT punch out of temperature values
		 * */
		if( GrainVar.iqp != NULL )
		{
			if( GrainVar.nqout1 == 0 )
			{

				for( nd=0; nd < NDUST; nd++ )
				{
					if( GrainVar.lgQHeat[nd] && GrainVar.lgDustOn1[nd] )
					{
						GrainVar.nqout1 += 1;
						GrainVar.nqp1[GrainVar.nqout1-1] = nd+1;
					}
				}

				if( GrainVar.nqout1 == 0 )
				{
					fprintf( GrainVar.iqp, "**** no grains were quantum heated.\n" );
					GrainVar.nqout1 = -1;
				}
			}

			if( GrainVar.nqout1 > 0 )
			{
				fprintf( GrainVar.iqp, "Dust Temperature Distribution: zone %4ld", 
				  ZoneCnt.nzone );
				fprintf( GrainVar.iqp, " --- dust types:" );
				for( nd=0; nd < GrainVar.nqout1; nd++ )
				{
					fprintf( GrainVar.iqp, "%3ld", GrainVar.nqp1[nd] );
				}
				fprintf( GrainVar.iqp, "\n" );

				fprintf( GrainVar.iqp, "Weighted mean temperatures: " );
				for( nd=0; nd < GrainVar.nqout1; nd++ )
				{
					fprintf( GrainVar.iqp, "%8.2f", GrainVar.tedust[GrainVar.nqp1[nd]-1] );
				}
				fprintf( GrainVar.iqp, "\n" );

				fprintf( GrainVar.iqp, "Equilibrium temperatures: " );
				for( nd=0; nd < GrainVar.nqout1; nd++ )
				{
					fprintf( GrainVar.iqp, "%8.2f", GrainVar.tequil1[GrainVar.nqp1[nd]-1] );
				}
				fprintf( GrainVar.iqp, "\n" );

				for( nq=0; nq < GrainVar.nqhg1; nq++ )
				{
					for( nd=0; nd < GrainVar.nqout1; nd++ )
					{
						qp1[nd] = GrainVar.qhprob[nq][GrainVar.nqp1[nd]-1];
						qt1[nd] = GrainVar.qtemp[nq][GrainVar.nqp1[nd]-1];
					}
					for( nd=0; nd < GrainVar.nqout1; nd++ )
					{
						fprintf( GrainVar.iqp, "%8.2f%10.3e", 
						  qt1[nd], qp1[nd] );
					}
					fprintf( GrainVar.iqp, "\n" );
				}
			}
		}
		/* End of changes here (K Volk) */
	}

	/* gas pressure, pressure due to incident radiation field, rad accel */
	fprintf( ioQQQ, " Pressure      NgasTgas");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.GasPres/BOLTZMANN));
	fprintf( ioQQQ, " P(total)");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.RadPres + pressure.GasPres));
	fprintf( ioQQQ, " P( gas )");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.GasPres));
	fprintf( ioQQQ, " P(Radtn)");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.RadPres));
	fprintf( ioQQQ, " Rad accl");
	fprintf(ioQQQ,PrintEfmt("%9.2e", wind.AccelTot));
	fprintf( ioQQQ, " ForceMul");
	fprintf(ioQQQ,PrintEfmt("%9.2e", wind.fmul));
	fprintf( ioQQQ, "\n" );

	/* heavy element molecules */
	if( hevmolec.hevmol[IPCO-1] > 0. )
	{
		fprintf( ioQQQ, " Molecules     CH/Ctot:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPCH-1]/xIonFracs[0][5]));
		fprintf( ioQQQ, " CH+/Ctot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPCHP-1]/xIonFracs[0][5]));
		fprintf( ioQQQ, " CO/Ctot:");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPCO-1]/xIonFracs[0][5]));
		fprintf( ioQQQ, " CO+/Ctot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPCOP-1]/xIonFracs[0][5]));
		fprintf( ioQQQ, " H2O/Otot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPH2O-1]/xIonFracs[0][7]));
		fprintf( ioQQQ, " OH/Ototl");
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPOH-1]/xIonFracs[0][7]));
		fprintf( ioQQQ, "\n");
	}

	/* Lithium, berylium */
	if( elmton.lgElmtOn[2] || elmton.lgElmtOn[3] )
	{
		fprintf( ioQQQ, " Lithium      " );
		for( i=1; i <= 4; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", xIonFracs[i][2]/MAX2(1e-35,xIonFracs[0][2]) ));
		}
		fprintf( ioQQQ, " Berylliu" );
		for( i=1; i <= 5; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e",  xIonFracs[i][3]/MAX2(1e-35,xIonFracs[0][3])) );
		}
		fprintf( ioQQQ, "\n" );

		/* option to print heating due to these stages*/
		if( pheat.lgPrintHeating )
		{
			fprintf( ioQQQ, "              " );
			for( i=0; i < 3; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e",  HeatingCom.heating[i][2]/ heat.htot) );
			}
			fprintf( ioQQQ, "                    " );

			for( i=0; i < 4; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[i][3]/heat.htot ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* borom*/
	if( elmton.lgElmtOn[4] )
	{
		fprintf( ioQQQ, " Boron        " );
		for( i=1; i <= 6; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", xIonFracs[i][4]/MAX2(1e-35,xIonFracs[0][4]) ));
		}
		fprintf( ioQQQ, "\n" );

		/* option to print heating*/
		if( pheat.lgPrintHeating )
		{ 
			fprintf( ioQQQ, "              " );
			for( i=0; i < 5; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[i][4]/heat.htot ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* Carbon */
	if( elmton.lgElmtOn[5] )
	{
		fprintf( ioQQQ, " Carbon       " );
		for( i=1; i <= 7; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", xIonFracs[i][5]/MAX2(1e-35,xIonFracs[0][5])) );
		}
		/* some molecules trail the line */
		fprintf( ioQQQ, " H2O+/O  " );
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPH2OP-1]/MAX2(1e-35,xIonFracs[0][7]) ));
		fprintf( ioQQQ, " OH+/Otot" );
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPOHP-1]/ MAX2(1e-35,xIonFracs[0][7]) ));
		/* print extra heating, normally zero */
		fprintf( ioQQQ, " Hex(tot)" );
		fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[20][0] ));
		fprintf( ioQQQ, "\n" );

		/* option to print heating*/
		if( pheat.lgPrintHeating )
		{
			fprintf( ioQQQ, "              " );
			for( i=0; i < 6; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[i][5]/ heat.htot) );
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* Nitrogen */
	if( elmton.lgElmtOn[6] )
	{
		fprintf( ioQQQ, " Nitrogen     " );
		for( i=1; i <= 8; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e",xIonFracs[i][6]/ MAX2(1e-35,xIonFracs[0][6]) ));
		}
		fprintf( ioQQQ, " O2/Ototl" );
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPOTWO-1]/MAX2(1e-35,xIonFracs[0][7])));
		fprintf( ioQQQ, " O2+/Otot" );
		fprintf(ioQQQ,PrintEfmt("%9.2e", hevmolec.hevmol[IPO2P-1]/ MAX2(1e-35,xIonFracs[0][7]) ));
		fprintf( ioQQQ, "\n" );

		/* option to print heating*/
		if( pheat.lgPrintHeating )
		{
			fprintf( ioQQQ, "              " );
			for( i=0; i < 7; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[i][6]/ heat.htot ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* now print rest of elements inside loops */

	/* Oxygen through Magnesium */
	for( nelem=7; nelem< 12; ++nelem )
	{
		if( elmton.lgElmtOn[nelem] )
		{
			/* print the element name and amount of shift */
			fprintf( ioQQQ, " %10.10s   ", elementnames.chElementName[nelem]);

			for( i=1; i <= nelem+2; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", xIonFracs[i][nelem]/xIonFracs[0][nelem] ));
			}
			if( nelem==7 )
			{
				/* print secondary ionization rate for oxygen */
				fprintf( ioQQQ, " sec ion:" );
					fprintf(ioQQQ,PrintEfmt("%9.2e", Secondaries.csupra) );
			}
			fprintf( ioQQQ, "\n" );

			/* print heating but only if needed */
			if( pheat.lgPrintHeating )
			{
				fprintf( ioQQQ, "              " );
				for( i=1; i <= nelem+2; i++ )
				{
					fprintf(ioQQQ,PrintEfmt("%9.2e", HeatingCom.heating[i-1][nelem]/heat.htot ));
				}
				fprintf( ioQQQ, "\n" );
			}
		}
	}

	/* Aluminium through Zinc */
	for( nelem=12; nelem< LIMELM; ++nelem )
	{
		if( elmton.lgElmtOn[nelem] )
		{
			/* number of ionization stages to print across the page */
			enum {LINE=13};
			ishift = MAX2(0,IonRange.IonHigh[nelem]-LINE);

			/* print the element name and amount of shift */
			fprintf( ioQQQ, " %10.10s%2ld ", elementnames.chElementName[nelem],ishift );

			for( i=1; i <= LINE; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", xIonFracs[i+ishift][nelem]/xIonFracs[0][nelem]) );
			}
			fprintf( ioQQQ, "\n" );

			/* print heating but only if needed */
			if( pheat.lgPrintHeating )
			{
				fprintf( ioQQQ, "              " );
				for( i=1; i <= LINE; i++ )
				{
					fprintf(ioQQQ,
						PrintEfmt("%9.2e", HeatingCom.heating[i+ishift-1][nelem]/heat.htot ));
				}
				fprintf( ioQQQ, "\n" );
			}
		}
	}

	/* call FeII print routine if large atom is turned on */
	if( FeII.lgFeIION )
	{
		FeIIPrint();
	}

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