/* 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 */
/*PunchDo produce punch output during calculation,
 * chTime is 'MIDL' during calculation, 'LAST' at the end */
/*PunchNewContinuum produce the 'punch new continuum' output */
/*PunchLineStuff punch optical depths or source functions for all transferred lines */
/*pun1Line called by PunchLineStuff to produce output for one line */
/*PunchNewContinuum produce the 'punch new continuum' output */
/*PunLineIntensity produce the 'punch lines intensity' output */
/* punch h emiss, for chapt 4, routine is below */
/*PunResults1Line do single line of output for the punch results and punch line intensity commands */
/* the number of emission lines across one line of printout */
/*PunchSpecial generate output for the punch special command */
/*punResults punch results from punch results command */
/*PunResults1Line do single line of output for the punch results and punch line intensity commands */
/*PunchGaunts called by punch gaunts command to output gaunt factors */
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "mean.h"
#include "taulines.h"
#include "struc.h"
#include "iso.h"
#include "mole.h"
#include "hyperfine.h"
#include "rt.h"
#include "lines_service.h"
#include "doppvel.h"
#include "dense.h"
#include "h2.h"
#include "magnetic.h"
#include "hydrogenic.h"
#include "gaunt.h"
#include "secondaries.h"
#include "grainvar.h"
#include "lines.h"
#include "dynamics.h"
#include "colden.h"
#include "ionbal.h"
#include "yield.h"
#include "prt.h"
#include "iterations.h"
#include "heavy.h"
#include "conv.h"
#include "geometry.h"
#include "called.h"
#include "helike.h"
#include "opacity.h"
#include "rfield.h"
#include "phycon.h"
#include "timesc.h"
#include "radius.h"
#include "atomfeii.h"
#include "assertresults.h"
#include "thermal.h"
#include "wind.h"
#include "hmi.h"
#include "pressure.h"
#include "elementnames.h"
#include "ipoint.h"
#include "gammas.h"
#include "atmdat.h"
#include "map.h"
#include "input.h"
#include "punch.h"
#include "grid.h"

/*PunResults1Line do single line of output for the punch results and punch line intensity commands */
/* the number of emission lines across one line of printout */
static void PunResults1Line(
  FILE * io, 
  /* 4 char + null string */
  const char *chLab, 
  float wl, 
  double xInten, 
  const char *chFunction);/* null terminated string saying what to do */

/*PunchGaunts called by punch gaunts command to output gaunt factors */
static void PunchGaunts(FILE* io);

/*punResults punch results from punch results command */
/*PunResults1Line do single line of output for the punch results and punch line intensity commands */
static void punResults(FILE* io);

static void PunchLineStuff(
  FILE * io,
  const char *chJob , 
  float xLimit);

/* punch h emiss, for chapt 4, routine is below */
static void AGN_Hemis(FILE *io );

/*PunchNewContinuum produce the 'punch new continuum' output */
static void PunchNewContinuum(FILE * io , long ipCon );

/*PunLineIntensity produce the 'punch lines intensity' output */
static void PunLineIntensity(FILE * io);

void PunchDo(
	/* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
	const char *chTime) 
{
	int lgOK, 
	  lgTlkSav;
	long int
	  ipPun, 
	  i,
	  i1, 
	  ion, 
	  ipConMax, 
	  ipHi,
	  ipLinMax, 
	  ipLo,
	  ips, 
	  j, 
	  jj, 
	  limit, 
	  nd, 
	  nelem;
	double ConMax, 
	  RateInter, 
	  av, 
	  conem, 
	  eps, 
	  flxatt, 
	  flxcor, 
	  flxin, 
	  flxref, 
	  flxtrn, 
	  fout, 
	  fref, 
	  fsum, 
	  opConSum, 
	  opLinSum, 
	  pop, 
	  stage, 
	  sum, 
	  texc, 
	  xLinMax;

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

	/* 
	 * the "last" option on punch command, to punch on last iteration,
	 * is parsed at the top of the loop in only one place.  
	 * no further action is needed at all for punch last to work
	 * ok throughout this routine 
	 */

	/* 
	 * each branch can have a test whether chTime is or is not "LAST"
	 *
	 * if( strcmp(chTime,"LAST") == 0 )  <== print after calculation is complete 
	 *
	 * if "LAST" then this is last call to routine after calc complete
	 * punch only if "LAST" when results at end of cal are needed
	 *
	 * if( strcmp(chTime,"LAST") != 0 )  <== print for every zone 
	 *
	 * test for .not."LAST" is for every zone result, where you do not
	 * want to punch last zone twice
	 */

	/* punch results on unit ipPnunit  */
	if( punch.npunch < 1 )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->PunchDo()\n", debug_fp );
#		endif
		return;
	}

	for( ipPun=0; ipPun < punch.npunch; ipPun++ )
	{
		/* this global variable to remember where in the punch stack we are */
		punch.ipConPun = ipPun;
		/* 
		 * lgPunLstIter set true if 'last' key occured on punch command
		 * normally is false.  This will skip punching if last set and
		 * this is not last iteration
		 */

		if( iterations.lgLastIt || (!punch.lgPunLstIter[ipPun]) )
		{

			if( strcmp(punch.chPunch[ipPun],"ABUN") == 0 )
			{
				/* punch abundances vs depth */
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], "%.2f", 
						log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) );
					for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
					{
						/* >>chng 05 feb 03, protect against non-positive abundances,
						 * bug caught by Marcelo Castellanos */
						fprintf( punch.ipPnunit[ipPun], "\t%.2f", 
						  log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) );
					}
					fprintf( punch.ipPnunit[ipPun], "\n" );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"21CM") == 0 )
			{
				/* punch information about 21 cm line */
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], 
					  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
					  /* depth, cm */
					  radius.depth_mid_zone,
					  hyperfine.Tspin21cm ,
					  phycon.te ,
					  /* temperature from Lya - 21 cm optical depth ratio */
					  3.84e-7* EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon /
					  SDIV( HFLines[0].TauCon ),
					  /*TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ),*/
					  HFLines[0].PopLo ,
					  HFLines[0].PopHi ,
					  OccupationNumberLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ),
					  HFLines[0].TauCon , 
					  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon,
					  HFLines[0].PopOpc,
					  /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */
					  (dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s])/SDIV( hyperfine.Tspin21cm),
					  /* why was above multiplied by this following term? */
					  /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/
					  HFLines[0].TauIn,
					  TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) ,
					  colden.H0_ov_Tspin
					  );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"AGES") == 0 )
			{
				/* punch timescales vs depth */
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], "%.5e%10.2e%10.2e%10.2e%10.2e\n", 
					  /* depth, cm */
					  radius.depth_mid_zone,
					  /* cooling timescale */
					  dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot, 
					  /* H2 formation timescale */
					  timesc.time_H2_Dest_here, 
					  /* CO formation timescale */
					  timesc.AgeCOMoleDest, 
					  /* H recombination timescale */
					  1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) );
				}
			}

			else if( strcmp(punch.chPunch[ipPun]," AGN") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					if( strcmp( punch.chPunchArgs[ipPun], "HECS" ) == 0 )
					{
						/* this routine is in helike.c */
						AGN_He1_CS(punch.ipPnunit[ipPun]);
					}
					if( strcmp( punch.chPunchArgs[ipPun], "HEMI" ) == 0 )
					{
						/* punch h emiss, for chapt 4, routine is below */
						AGN_Hemis(punch.ipPnunit[ipPun]);
					}
					else
					{
						fprintf( ioQQQ, " PunchDo does not recognize flag %4.4s set for AGN punch.  This is impossible.\n", 
						  punch.chPunch[ipPun] );
						ShowMe();
						puts( "[Stop in PunchDo]" );
						cdEXIT(EXIT_FAILURE);
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"ASSE") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch the assert output */
					/*lint -e534 ignore return value of following  - nothing to do with it */
					lgCheckAsserts(punch.ipPnunit[ipPun]);
					/*lint +e534 ignore return value of following  - nothing to do with it */
				}
			}

			else if( strncmp(punch.chPunch[ipPun],"CHA",3) == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* one of the charge transfer options, all in chargtran.c */
					ChargTranPun( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"COOL") == 0 )
			{
				/* punch cooling, routine in file of same name */
				if( strcmp(chTime,"LAST") != 0 )
					CoolPunch(punch.ipPnunit[ipPun]);
			}

			else if( strncmp(punch.chPunch[ipPun],"DYN" , 3) == 0 )
			{
				/* punch information about dynamics solutions */
				if( strcmp(chTime,"LAST") != 0 )
					DynaPunch(punch.ipPnunit[ipPun] ,punch.chPunch[ipPun][3] );
			}

			else if( strcmp(punch.chPunch[ipPun],"ENTH") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
					fprintf( punch.ipPnunit[ipPun],
						"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
						radius.depth_mid_zone,
						phycon.EnthalpyDensity,
						phycon.EnergyExcitation,
						phycon.EnergyIonization,
						phycon.EnergyBinding ,
						0.5*POW2(wind.windv)*dense.xMassDensity ,	/* KE */
						5./2.*pressure.PresGasCurr ,				/* thermal plus PdV work */
						magnetic.EnthalpyDensity);						/* magnetic terms */
			}

			else if( strcmp(punch.chPunch[ipPun],"COLU") == 0 )
			{
				/* punch column densities */
				if( strcmp(chTime,"LAST") == 0 )
					PrtColumns(punch.ipPnunit[ipPun]);
			}

			else if( strcmp(punch.chPunch[ipPun],"COMP") == 0 )
			{
				/* compton energy exchange coefficients */
				if( strcmp(chTime,"LAST") != 0 )
				{
					for( jj=0; jj<rfield.nflux; jj = jj + punch.ncPunchSkip)
					{
						fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e\n", 
						  rfield.anu[jj], rfield.comup[jj]/rfield.widflx[jj], 
						  rfield.comdn[jj]/rfield.widflx[jj] );
					}
				}
			}

			/* punch continuum commands */
			else if( strcmp(punch.chPunch[ipPun],"CON ") == 0 )
			{
				/* this branch is the usual "punch continuum" case */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */
				if( strcmp(chTime,"LAST") == 0 )
				{
					for( j=0; j<rfield.nflux ; j = j+punch.ncPunchSkip)
					{
						/* four continua predicted here;
						 * incident, attenuated incident, emitted,
						 * then attenuated incident + emitted, last reflected
						 * reflected continuum is stored relative to inner radius
						 * others are stored for this radius */

						/* NB this code also used in punch emitted,
						 * transmitted continuum commands */

						/* the incident continuum */
						flxin = rfield.FluxSave[j]*rfield.anu2[j]*
						  EN1RYD/rfield.widflx[j];

						/* the reflected continuum */
						flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/rfield.widflx[j] +
							rfield.anu[j]*punch.PunchLWidth*rfield.reflin[j])*EN1RYD;

						/* the attenuated incident continuum */
						flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/
						  rfield.widflx[j]*radius.r1r0sq;

						/* the outward emitted continuum */
						conem = ((rfield.ConEmitOut[j])/
						  rfield.widflx[j]*rfield.anu2[j] + punch.PunchLWidth*
						  rfield.outlin[j]*rfield.anu[j])*radius.r1r0sq*
						  EN1RYD*geometry.covgeo;

						/* sum of emitted and transmitted continua */
						flxtrn = conem + flxatt;

						/* photon energy in appropriate energy or wavelength units */
						fprintf( punch.ipPnunit[ipPun],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) );
						/* incident continuum */
						fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxin ); 
						/* trans cont */
						fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxatt ); 
						/* DiffOut cont */
						fprintf( punch.ipPnunit[ipPun],"%.3e\t", conem ); 
						/* net trans cont */
						fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn ); 
						/* reflc cont */
						fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref ); 
						/* total cont */
						fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxref + flxtrn ); 
						fprintf( punch.ipPnunit[ipPun], "%4.4s\t%4.4s\t", 
						/* line	label */
						  rfield.chLineLabel[j] ,
						/* cont label*/
						  rfield.chContLabel[j] );
						/* number of lines within that cell over cell width
						 * punch raw continuum has number by itself */
						fprintf( punch.ipPnunit[ipPun], "%.2f\n", rfield.line_count[j]/rfield.widflx[j]*rfield.anu[j] );
					}
				}
			}

			/* this is the punch spectrum command, 
			 * the new continuum command that will replace the previous one */
			else if( strcmp(punch.chPunch[ipPun],"CONN") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
					/* io unit and which new continuum this is (was set when punch read in */
					PunchNewContinuum( punch.ipPnunit[ipPun] , (long)punch.punarg[0][ipPun]);
			}

			else if( strcmp(punch.chPunch[ipPun],"CONC") == 0 )
			{
				/* punch incident continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* incident continuum */
					for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
					{
						flxin = rfield.FluxSave[j]*rfield.anu2[j]*
						  EN1RYD/rfield.widflx[j];
						/* >>chng 96 oct 22, format of anu to 11.5 to resolve grid near 1 ryd */
						fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\n", 
						  AnuUnit(rfield.AnuOrg[j]), flxin );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONG") == 0 )
			{
				/* punch emitted grain continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* tried to remember emission from the two grain types in rtdiffuse */
					for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
					{
						fsum = gv.GraphiteEmission[j]*rfield.anu2[j]*
						  EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;

						fout = gv.SilicateEmission[j]*rfield.anu2[j]*
						  EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;

						/* >>chng 96 oct 22, format of anu to 11.5 to resolve grid near 1 ryd */
						fprintf( punch.ipPnunit[ipPun], "%12.4e %10.3e %10.3e\n", 
						  AnuUnit(rfield.AnuOrg[j]) , fsum , fout );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONR") == 0 )
			{
				/* punch reflected continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */
				if( strcmp(chTime,"LAST") == 0 )
				{
					if( geometry.lgSphere )
					{
						fprintf( punch.ipPnunit[ipPun], " Reflected continuum not predicted when SPHERE is set.\n" );
						fprintf( ioQQQ , 
							"\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" );
						puts( "[Stop in PunchDo]" );
						cdEXIT(EXIT_FAILURE);
					}

					for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
					{
						/* reflected continuum */
						flxref = rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/
						  rfield.widflx[j]*EN1RYD;
						/* reflected lines */
						fref = rfield.anu[j]*punch.PunchLWidth*
						  rfield.reflin[j]*EN1RYD;
						/* ratio of reflected to incident continuum, the albedo */
						if( rfield.FluxSave[j] > 1e-25 )
						{
							av = rfield.ConRefIncid[j]/rfield.FluxSave[j];
						}
						else
						{
							av = 0.;
						}
						/* >>chng 96 oct 22, format of anu to 12.5 to resolve grid near 1 ryd */
						fprintf( punch.ipPnunit[ipPun], "%12.5e%12.4e%12.4e%12.4e%12.4e %4.4s\n", 
						  AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref, 
						  av, rfield.chContLabel[j] );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CNVE") == 0 )
			{
				/* the punch convergence error command */
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], 
						"%.5e\t%li\t%.4e\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n", 
						radius.depth_mid_zone, 
						conv.nPres2Ioniz,
						pressure.PresTotlCorrect, 
						pressure.PresTotlCurr, 
						(pressure.PresTotlCorrect - pressure.PresTotlCurr)*100./pressure.PresTotlCorrect, 
						dense.EdenTrue,
						dense.eden,
						(dense.EdenTrue - dense.eden)*100./dense.EdenTrue,
						thermal.htot,
						thermal.ctot,
						(thermal.htot - thermal.ctot)*100./thermal.htot );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONB") == 0 )
			{
				/* punch continuum bins binning */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */
				if( strcmp(chTime,"LAST") != 0 )
				{
					for( j=0; j<rfield.nupper; j = j + punch.ncPunchSkip)
					{
						fprintf( punch.ipPnunit[ipPun], "%14.5e %14.5e %14.5e\n", 
						  AnuUnit(rfield.AnuOrg[j]), rfield.anu[j], rfield.widflx[j] );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"COND") == 0 )
			{
				/* punch diffuse continuum the local thermal emission */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* this option to punch diffuse continuum */
					for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
					{
						/* >>chng 96 oct 22, format of anu to 12.5 to resolve grid near 1 ryd */
						fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\n", 
						  AnuUnit(rfield.AnuOrg[j]), 
						  rfield.ConEmitLocal[j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONE") == 0 )
			{
				/* punch emitted continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch emitted continuum */
					for( j=0; j<rfield.nflux;j = j +punch.ncPunchSkip)
					{
						/* this is the reflected component */
						flxref = (rfield.anu2[j]*(rfield.ConRefIncid[j]+rfield.ConEmitReflec[j])/
						  rfield.widflx[j] + rfield.anu[j]*punch.PunchLWidth*
						  rfield.reflin[j])*EN1RYD;

						/* this is the total emission in the outward direction */
						conem = (rfield.ConEmitOut[j])/rfield.widflx[j]*rfield.anu2[j] + 
							punch.PunchLWidth*rfield.outlin[j]*rfield.anu[j];

						conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;

						/* output: photon energy, reflected, outward, total emission
						 *  >>chng 96 oct 22, format of anu to 11.5 to resolve grid near 1 ryd */
						fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n", 
						  AnuUnit(rfield.AnuOrg[j]), 
						  flxref, 
						  conem, 
						  flxref + conem, 
						  rfield.chLineLabel[j], 
						  rfield.chContLabel[j]
						   );
					}
				}
			}

			/* punch fine continuum command */
			else if( strcmp(punch.chPunch[ipPun],"CONf") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					long nu_hi , nskip;
					if( punch.punarg[0][ipPun] > 0. )
					{
						/* possible lower bounds to energy range - will be zero if not set */
						j = ipFineCont( punch.punarg[0][ipPun] );
					}
					else
					{
						j = 0;
					}

					nskip = (long)punch.punarg[2][ipPun];

					/* upper limit */
					if( punch.punarg[1][ipPun]> 0. )
					{
						nu_hi = ipFineCont( punch.punarg[1][ipPun]);
					}
					else
					{
						nu_hi = rfield.nfine;
					}
					do
					{
						float sum1 = rfield.fine_opt_depth[j];
						float xnu = rfield.fine_anu[j+nskip/2];
						for( jj=1; jj<nskip; ++jj )
						{
							sum1 += rfield.fine_opt_depth[j+jj];
						}
						fprintf( punch.ipPnunit[ipPun], 
							"%.6e\t%.3e\n", 
							AnuUnit(xnu), 
							sexp(sum1/nskip) );
						j += nskip;
					} while( j < nu_hi );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONi") == 0 )
			{
				/* punch continuum interactions */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */

				/* continuum interactions */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* this is option to set lowest energy */
					if( punch.punarg[0][ipPun] <= 0. )
					{
						i1 = 1;
					}
					else if( punch.punarg[0][ipPun] < 100. )
					{
						i1 = ipoint(punch.punarg[0][ipPun]);
					}
					else
					{
						i1 = (long int)punch.punarg[0][ipPun];
					}

					fref = 0.;
					fout = 0.;
					fsum = 0.;
					sum = 0.;
					flxin = 0.;

					for( j=i1-1; j < rfield.nflux; j++ )
					{
						fref += rfield.flux[j]*opac.opacity_abs[j];
						fout += rfield.otslin[j]*opac.opacity_abs[j];
						fsum += rfield.otscon[j]*opac.opacity_abs[j];
						sum += rfield.ConInterOut[j]*opac.opacity_abs[j];
						flxin += (rfield.outlin[j] + rfield.outlin_noplot[j])*opac.opacity_abs[j];
					}
					fprintf( punch.ipPnunit[ipPun], "%10.2e%10.2e%10.2e%10.2e%10.2e\n", 
					  fref, fout, fsum, sum, flxin );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONI") == 0 )
			{
				/* punch ionizing continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */

				if( (punch.punarg[2][ipPun]>0) || (strcmp(chTime,"LAST") == 0) )
				{
					/* this flag will remember whether we have ever printed anything */
					int lgPrt=FALSE;
					if( punch.punarg[2][ipPun]>0 )
						fprintf(punch.ipPnunit[ipPun],"#punch every zone %li\n", nzone);

					/* punch ionizing continuum command
					 * this is option to set lowest energy,
					 * if no number was entered then this was zero */
					if( punch.punarg[0][ipPun] <= 0. )
					{
						i1 = 1;
					}
					else if( punch.punarg[0][ipPun] < 100. )
					{
						i1 = ipoint(punch.punarg[0][ipPun]);
					}
					else
					{
						i1 = (long int)punch.punarg[0][ipPun];
					}

					sum = 0.;
					for( j=i1-1; j < rfield.nflux; j++ )
					{
						flxcor = rfield.flux[j] + 
						  rfield.otslin[j] + 
						  rfield.otscon[j] + 
						  rfield.ConInterOut[j] +
						  rfield.outlin[j] + rfield.outlin_noplot[j];

						sum += flxcor*opac.opacity_abs[j];
					}

					if( sum > 0. )
					{
						sum = 1./sum;
					}
					else
					{
						sum = 1.;
					}

					fsum = 0.;

					for( j=i1-1; j<rfield.nflux ; ++j)
					{
						flxcor = rfield.flux[j] + 
						  rfield.otslin[j] + 
						  rfield.otscon[j] + 
						  rfield.ConInterOut[j]+
						  rfield.outlin[j] + rfield.outlin_noplot[j];

						fsum += flxcor*opac.opacity_abs[j];

						/* punched quantities are freq, flux, flux*cross sec,
						 * fraction of total, integral fraction of total */
						RateInter = flxcor*opac.opacity_abs[j]*sum;

						/* punage(ipPun,2) is lowest interaction rate to consider, def=0.01 (1 percnt) */
						/* >>chng 01 nov 22, format to c-friendly */
						if( (RateInter >= punch.punarg[1][ipPun]) && (flxcor > SMALLFLOAT) )
						{
							lgPrt = TRUE;
							/* >>chng 96 oct 22, format of anu to 11.5 to resolve grid near 1 ryd */
							fprintf( punch.ipPnunit[ipPun], 
								"%li\t%.5e\t%.2e\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2e\t%.2e\t%.4s\t%.4s\n", 
							  j,
							  AnuUnit(rfield.AnuOrg[j]), 
							  flxcor, 
							  flxcor*opac.opacity_abs[j], 
							  rfield.flux[j]/flxcor, 
							  rfield.otslin[j]/flxcor, 
							  rfield.otscon[j]/flxcor, 
							  (rfield.outlin[j] + rfield.outlin_noplot[j])/flxcor, 
							  rfield.ConInterOut[j]/flxcor, 
							  RateInter, 
							  fsum*sum, 
							  rfield.chLineLabel[j], 
							  rfield.chContLabel[j] );
						}
					}
					if( !lgPrt )
					{
						/* entered logical block but did not print anything */
						fprintf(punch.ipPnunit[ipPun],
							" punchdo, the PUNCH IONIZING CONTINUUM command did not find a strong point, sum fsum were %.2e %.2e\n",
							sum,fsum);
						fprintf(punch.ipPnunit[ipPun],
							" punchdo, the low-frequency energy was %.5e Ryd\n",rfield.anu[i1-1]);
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CORA") == 0 )
			{
				/* punch raw continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */

				if( strcmp(chTime,"LAST") == 0 )
				{
					/* this option to punch all raw ionizing continuum */
					for( j=0;j<rfield.nflux;j = j + punch.ncPunchSkip)
					{
						fprintf( punch.ipPnunit[ipPun], 
							"%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t", 
						  AnuUnit(rfield.AnuOrg[j]), 
						  rfield.flux[j], 
						  rfield.otslin[j], 
						  rfield.otscon[j], 
						  rfield.ConRefIncid[j],
						  rfield.ConEmitReflec[j], 
						  rfield.ConInterOut[j],
						  rfield.outlin[j]+rfield.outlin_noplot[j], 
						  rfield.ConEmitOut[j] ,
						  rfield.chLineLabel[j], 
						  rfield.chContLabel[j]
						  );
						/* number of lines within that cell */
						fprintf( punch.ipPnunit[ipPun], "%li\n", rfield.line_count[j] );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONo") == 0 )
			{
				/* punch outward local continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */

				if( strcmp(chTime,"LAST") == 0 )
				{
					/* option to punch out outward continuum here */
					for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip)
					{
						fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e\n", 
						  AnuUnit(rfield.AnuOrg[j]), 
						  rfield.ConEmitOut[j]+ rfield.outlin[j] + rfield.outlin_noplot[j], 
						  (rfield.flux[j] + rfield.otscon[j] + rfield.otslin[j] + 
						  rfield.ConInterOut[j])*opac.opacity_abs[j]*
						  rfield.anu[j] );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONO") == 0 )
			{
				/* punch outward continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */

				if( strcmp(chTime,"LAST") == 0 )
				{
					/* option to punch out outward continuum */
					for( j=0; j<rfield.nflux; j = j + punch.ncPunchSkip)
					{
						fprintf( punch.ipPnunit[ipPun], "%11.5e%10.2e%10.2e%10.2e%10.2e\n", 
						  AnuUnit(rfield.AnuOrg[j]), 
						  rfield.flux[j]*rfield.anu2[j]* EN1RYD/rfield.widflx[j]*radius.r1r0sq, 
						  rfield.ConInterOut[j]/rfield.widflx[j]*rfield.anu2[j]*radius.r1r0sq*EN1RYD, 
						  punch.PunchLWidth* (rfield.outlin[j]+rfield.outlin_noplot[j])*rfield.anu[j]*radius.r1r0sq*EN1RYD, 
						  (rfield.ConInterOut[j]/rfield.widflx[j]*
						  rfield.anu2[j] + punch.PunchLWidth*(rfield.outlin[j]+rfield.outlin_noplot[j])*
						  rfield.anu[j])*radius.r1r0sq*EN1RYD );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"CONT") == 0 )
			{
				/* punch transmitted continuum */
				/* set pointer for possible change in units of energy in continuum
				 * AnuUnit will give anu in whatever units were set with punch units */

				if( strcmp(chTime,"LAST") == 0 )
				{
					/* this option to punch transmitted continuum */
					for( j=0;j<rfield.nflux; j = j + punch.ncPunchSkip)
					{
						/* attenuated incident continuum
						 * >>chng 97 jul 10, remove PunchLWidth from this one only since
						 * we must conserve energy even in lines */
						flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/
						  rfield.widflx[j]*radius.r1r0sq;

						/*conem = (rfield.ConOutNoInter[j] + rfield.ConInterOut[j]+rfield.outlin[j])*
						  rfield.anu2[j] ;
						conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;*/
						/* >>chng 00 jan 03, above did not include all contributors.  
						 * Pasted in below from usual
						 * punch continuum command */
						/* >>chng 04 jul 15, removed factor of punch.PunchLWidth -
						 * this should not be there to conserve energy, as explained in hazy
						 * where command was documented, and in comment above.  caught by PvH */
						/* >>chng 04 jul 23, incorrect use of outlin - before mult by an2,
						 * quantity should be photons per Ryg, since init quantity is
						 * photons per cell.  Must div by widflx.  caught by PvH  */
						/*conem = (rfield.ConEmitOut[j]/rfield.widflx[j]*rfield.anu2[j] +
						  rfield.outlin[j]*rfield.anu[j])*radius.r1r0sq*
						  EN1RYD*geometry.covgeo;*/
						conem = (rfield.ConEmitOut[j] + rfield.outlin[j]) / rfield.widflx[j]*
							rfield.anu2[j]*radius.r1r0sq*EN1RYD*geometry.covgeo;

						flxtrn = conem + flxatt;

						/* use AnuOrg here instead of anu since probably
						 * going to be used by table read
						 * and we want good anu array for sanity check*/
						fprintf( punch.ipPnunit[ipPun],"%10.3e\t", AnuUnit(rfield.AnuOrg[j]) );
						fprintf( punch.ipPnunit[ipPun],"%10.3e\t", flxtrn );
						fprintf( punch.ipPnunit[ipPun],"%10.3e\n", rfield.trans_coef[j] );
					}
				}
			}


			else if( strcmp(punch.chPunch[ipPun],"CON2") == 0 )
			{
				/* punch total two-pohoton continuum  */
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* this option to punch diffuse continuum */
					for( j=0; j<rfield.nflux; j = j+punch.ncPunchSkip)
					{
						fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\n", 
						  AnuUnit(rfield.AnuOrg[j]), 
						  rfield.TotDiff2Pht[j]/rfield.widflx[j] , 
						  rfield.TotDiff2Pht[j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUSE") == 0 )
			{
				/* punch grain extinction - includes only grain opacity, not total */
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], " %.5e\t", 
						radius.depth_mid_zone );

					/* visual extinction of an extended source (like a PDR)*/
					fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);

					/* visual extinction of point source (star)*/
					fprintf( punch.ipPnunit[ipPun], "%.2e\n" , rfield.extin_mag_V_point);
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUSO") == 0 )
			{
				/* punch grain opacity */
				if( strcmp(chTime,"LAST") == 0 )
				{
					for( j=0; j < rfield.nflux; j++ )
					{
						double scat;
						fprintf( punch.ipPnunit[ipPun], 
						  "%.5e\t%.2e\t%.2e\t%.2e\t", 
						  /* photon energy or wavelength */
						  AnuUnit(rfield.AnuOrg[j]), 
						  /* total opacity, discount forward scattering */
						  gv.dstab[j] + gv.dstsc[j], 
						  /* absorption opacity */
						  gv.dstab[j], 
						  /* scatter, with forward discounted */
						  gv.dstsc[j] );
						/* add together total scattering, discounting 1-g */
						scat = 0.;
						/* sum over all grain species */
						for( nd=0; nd < gv.nBin; nd++ )
						{
							scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund;
						}
						/* finally, scattering including effects of forward scattering */
						fprintf( punch.ipPnunit[ipPun], 
						  "%.2e", scat );
						fprintf( punch.ipPnunit[ipPun], 
						  "%.2e\n", gv.dstsc[j]/(gv.dstab[j] + gv.dstsc[j]) );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUSP") == 0 )
			{
				/* grain potential */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* used to print header exactly one time */
					static int lgMustPrtHeader = TRUE;
					/* punch grain velocity, but do labels first if this is first zone */
					if( lgMustPrtHeader )
					{
						/* first print string giving grain id */
						fprintf( punch.ipPnunit[ipPun], "#Depth" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
						fprintf( punch.ipPnunit[ipPun], "\n" );

						/* now print grain radius converting from cm to microns */
						fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
						fprintf( punch.ipPnunit[ipPun], "\n" );
						
						/* don't need to do this, ever again */
						lgMustPrtHeader = FALSE;
					}
					fprintf( punch.ipPnunit[ipPun], " %.5e", 
						radius.depth_mid_zone );
					/* grain potential in eV */
					for( nd=0; nd < gv.nBin; ++nd ) 
						fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD );
					fprintf( punch.ipPnunit[ipPun], "\n" );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUSR") == 0 )
			{
				/* grain H2 formation rates */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* used to print header exactly one time */
					static int lgMustPrtHeader = TRUE;
					/* punch grain velocity, but do labels first if this is first zone */
					if( lgMustPrtHeader )
					{
						/* first print string giving grain id */
						fprintf( punch.ipPnunit[ipPun], "#Depth" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
						fprintf( punch.ipPnunit[ipPun], "\n" );

						/* now print grain radius converting from cm to microns */
						fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
						fprintf( punch.ipPnunit[ipPun], "\n" );
						
						/* don't need to do this, ever again */
						lgMustPrtHeader = FALSE;
					}
					fprintf( punch.ipPnunit[ipPun], " %.5e", 
						radius.depth_mid_zone );
					/* grain formation rate for H2 */
					for( nd=0; nd < gv.nBin; ++nd ) 
						fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used );
					fprintf( punch.ipPnunit[ipPun], "\n" );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUST") == 0 )
			{
				/* grain temperatures - K*/
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* used to print header exactly one time */
					static int lgMustPrtHeader = TRUE;
					/* punch grain temperature, but do labels first if this is first zone */
					if( lgMustPrtHeader )
					{
						/* first print string giving grain id */
						fprintf( punch.ipPnunit[ipPun], "#Depth" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
						fprintf( punch.ipPnunit[ipPun], "\n" );

						/* now print grain radius converting from cm to microns */
						fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
						fprintf( punch.ipPnunit[ipPun], "\n" );

						/* don't need to do this, ever again */
						lgMustPrtHeader = FALSE;
					}
					fprintf( punch.ipPnunit[ipPun], " %.5e", 
						radius.depth_mid_zone );
					for( nd=0; nd < gv.nBin; ++nd ) 
						fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust );
					fprintf( punch.ipPnunit[ipPun], "\n" );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUSC") == 0 )
			{
				/* punch charge per grain in electrons / grain */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* used to print header exactly one time */
					static int lgMustPrtHeader = TRUE;
					/* punch grain charge, but do labels first if this is first zone */
					if( lgMustPrtHeader )
					{
						/* first print string giving grain id */
						fprintf( punch.ipPnunit[ipPun], "#Depth\tne(grn)" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
						fprintf( punch.ipPnunit[ipPun], "\n" );

						/* now print grain radius converting from cm to microns */
						fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)\tne\t" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
						fprintf( punch.ipPnunit[ipPun], "\n" );
						
						/* don't need to do this, ever again */
						lgMustPrtHeader = FALSE;
					}
					fprintf( punch.ipPnunit[ipPun], " %.5e\t%.4e", 
						radius.depth_mid_zone ,
						gv.TotalEden );
					/* average charge per grain in electrons */
					for( nd=0; nd < gv.nBin; ++nd ) 
						fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ );
					fprintf( punch.ipPnunit[ipPun], "\n" );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUSH") == 0 )
			{
				/* grain heating */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* used to print header exactly one time */
					static int lgMustPrtHeader = TRUE;
					/* punch grain charge, but do labels first if this is first zone */
					if( lgMustPrtHeader )
					{
						/* first print string giving grain id */
						fprintf( punch.ipPnunit[ipPun], "#Depth" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
						fprintf( punch.ipPnunit[ipPun], "\n" );

						/* now print grain radius converting from cm to microns */
						fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
						fprintf( punch.ipPnunit[ipPun], "\n" );
						
						/* don't need to do this, ever again */
						lgMustPrtHeader = FALSE;
					}
					fprintf( punch.ipPnunit[ipPun], " %.5e", 
						radius.depth_mid_zone );
					/* grain heating */
					for( nd=0; nd < gv.nBin; ++nd ) 
						fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl );
					fprintf( punch.ipPnunit[ipPun], "\n" );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"DUSV") == 0 )
			{
				/* grain drift velocities */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* used to print header exactly one time */
					static int lgMustPrtHeader = TRUE;
					/* punch grain velocity, but do labels first if this is first zone */
					if( lgMustPrtHeader )
					{
						/* first print string giving grain id */
						fprintf( punch.ipPnunit[ipPun], "#Depth" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
						fprintf( punch.ipPnunit[ipPun], "\n" );

						/* now print grain radius converting from cm to microns */
						fprintf( punch.ipPnunit[ipPun], "#grn rad (mic)" );
						for( nd=0; nd < gv.nBin; ++nd ) 
							fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
						fprintf( punch.ipPnunit[ipPun], "\n" );
						
						/* don't need to do this, ever again */
						lgMustPrtHeader = FALSE;
					}
					fprintf( punch.ipPnunit[ipPun], " %.5e", 
						radius.depth_mid_zone );
					/* grain drift velocity in km/s */
					for( nd=0; nd < gv.nBin; ++nd ) 
						fprintf( punch.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 );
					fprintf( punch.ipPnunit[ipPun], "\n" );
				}
			}

			/* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
			else if( strcmp(punch.chPunch[ipPun],"DUSQ") == 0 )
			{
				/* punch grain Qs */
				if( strcmp(chTime,"LAST") == 0 )
				{
					for( j=0; j < rfield.nflux; j++ )
					{
						fprintf( punch.ipPnunit[ipPun], "%11.4e", 
						  rfield.anu[j] );
						for( nd=0; nd < gv.nBin; nd++ )
						{
							fprintf( punch.ipPnunit[ipPun], "%9.1e%9.1e", 
							   gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea,
							   gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea );
						}
						fprintf( punch.ipPnunit[ipPun], "\n" );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"ELEM") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					float renorm = 1.f;

					/* this is the index for the atomic number on the physical scale */
					/* >>chng 04 nov 23, use c scale throughout */
					nelem = (long int)punch.punarg[0][ipPun];
					ASSERT( nelem >= ipHYDROGEN );

					/* don't do this if element is not turned on */
					if( dense.lgElmtOn[nelem] )
					{
						/* >>chng 04 nov 23, add density option, leave as cm-3 
						* default is still norm to total of that element */
						if( punch.punarg[1][ipPun] == 0 )
							renorm = dense.gas_phase[nelem];

						fprintf( punch.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone );

						for( j=0; j <= (nelem + 1); ++j)
						{
							fprintf( punch.ipPnunit[ipPun], "\t%.2e", 
							dense.xIonDense[nelem][j]/renorm );
						}
						if( nelem==ipHYDROGEN )
						{
							/* H2 */
							fprintf( punch.ipPnunit[ipPun], "\t%.2e", 
								hmi.H2_total/renorm );
						}
						/* >>chng 04 nov 23 add C and O fine structure pops */
						else if( nelem==ipCARBON )
						{
							fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e", 
								colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm);
							fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e", 
								colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm);
							fprintf( punch.ipPnunit[ipPun], "\t%.2e", 
								co.hevmol[ipCO]/renorm );
						}
						else if( nelem==ipOXYGEN )
						{
							fprintf( punch.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e", 
								colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm);
						}
						fprintf( punch.ipPnunit[ipPun], "\n" );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"RECA") == 0 )
			{
				/* this will create table for agn then exit, 
				 * routine is in makerecom.c */
				ion_recombAGN( punch.ipPnunit[ipPun] );
				cdEXIT(EXIT_FAILURE);
			}

			else if( strcmp(punch.chPunch[ipPun],"RECE") == 0 )
			{
				/* punch recombination efficiencies, 
				 * option turned on with the  "punch recombination efficiencies" command
				 * output for the punch recombination coefficients command is actually
				 * produced by a series of routines, as they generate the recombination
				 * coefficients.  these include 
				 * dielsupres, helium, hydrorecom, iibod, and makerecomb*/
				fprintf( punch.ipPnunit[ipPun], 
					"%12.4e %12.4e %12.4e %12.4e\n", 
				  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][0][ipRecRad], 
				  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][0][ipRecNetEsc] ,
				  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecRad],
				  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][2][ipRecNetEsc]);

			}

			else if( strcmp(punch.chPunch[ipPun],"FEco") == 0 )
			{
				/* FeII continuuum, check that FeII is turned on */
				if( strcmp(chTime,"LAST") == 0 )
				{
					for( j=0; j < nFeIIConBins; j++ )
					{
						/* [j][0] is flux, [j][1] and [j][2] are upper and
						 * lower bounds in Angstroms.  
						 * these are set in FeIIZero */
						/* emission from large FeII atom, integrated over continuum intervals 
						 * units are always in units, erg cm-2 s-1 for intensity case */
						fprintf( punch.ipPnunit[ipPun], "%.2f\t%e \n", 
							(FeII_Cont[j][1]+FeII_Cont[j][2])/2.,
							FeII_Cont[j][0] );
					}
				}
			}

			/* punch column densities */
			else if( strcmp(punch.chPunch[ipPun],"FEcl") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch FeII level energies and stat weights, followed by column density */
					FeIIPunchColden( punch.ipPnunit[ipPun] );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"FE2l") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch FeII level energies and stat weights */
					FeIIPunchLevels( punch.ipPnunit[ipPun] );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"FEli") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch line intensities */
					FeIIPunchLines( punch.ipPnunit[ipPun] );
				}
			}
			else if( strcmp(punch.chPunch[ipPun],"FE2f") == 0 )
			{
				/* set with punch feii fred command, this punches some stuff from
				* fred hamann's feii atom */
				if( strcmp(chTime,"LAST") != 0 )
					pfeii(punch.ipPnunit[ipPun]);

			}
			else if( strcmp(punch.chPunch[ipPun],"FE2d") == 0 )
			{
				/* punch some departure cefficients for large feii atom */
				if( strcmp(chTime,"LAST") != 0 )
					FeIIPunDepart(punch.ipPnunit[ipPun],FALSE);

			}
			else if( strcmp(punch.chPunch[ipPun],"FE2D") == 0 )
			{
				/* punch all departure cefficients for large feii atom */
				if( strcmp(chTime,"LAST") != 0 )
					FeIIPunDepart(punch.ipPnunit[ipPun],TRUE);

			}
			else if( strcmp(punch.chPunch[ipPun],"FE2p") == 0 )
			{
				/* punch level populations for large feii atom */
				if( strcmp(chTime,"LAST") != 0 )
					FeIIPunPop(punch.ipPnunit[ipPun],FALSE);

			}
			else if( strcmp(punch.chPunch[ipPun],"FE2P") == 0 )
			{
				/* punch all level populations for large feii atom */
				if( strcmp(chTime,"LAST") != 0 )
					FeIIPunPop(punch.ipPnunit[ipPun],TRUE);

			}
			/* punch spectra in fits format */
			else if( strcmp(punch.chPunch[ipPun],"FITS") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					punchFITSfile( punch.ipPnunit[ipPun], "FITS" );
				}
			}
			/* punch gammas (but without element) */
			else if( strcmp(punch.chPunch[ipPun],"GAMt") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					long ns;
					/* punch photoionization rates, with the PUNCH GAMMAS command */
					for( nelem=0; nelem < LIMELM; nelem++ )
					{
						for( ion=0; ion <= nelem; ion++ )
						{
							for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
							{
								fprintf( punch.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e", 
									nelem+1, ion+1, ns+1, 
									ionbal.PhotoRate_Shell[nelem][ion][ns][0], 
									ionbal.PhotoRate_Shell[nelem][ion][ns][1] ,
									ionbal.PhotoRate_Shell[nelem][ion][ns][2] );

								for( j=0; j < yield.n_elec_eject[nelem][ion][ns]; j++ )
								{
									fprintf( punch.ipPnunit[ipPun], "%5.2f", yield.frac_elec_eject[nelem][ion][ns][j] );
								}
								fprintf( punch.ipPnunit[ipPun], "\n" );
							}
						}
					}
				}
			}
			/* punch gammas element, ion */
			else if( strcmp(punch.chPunch[ipPun],"GAMe") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					int ns;
					nelem = (long)punch.punarg[0][ipPun];
					ion = (long)punch.punarg[1][ipPun];
					/* valence shell */
					ns = Heavy.nsShells[nelem][ion]-1;
					/* show what some of the ionization sources are */
					GammaPrt( 
						opac.ipElement[nelem][ion][ns][0] , 
						opac.ipElement[nelem][ion][ns][1] , 
						opac.ipElement[nelem][ion][ns][2] , 
						punch.ipPnunit[ipPun], 
						ionbal.PhotoRate_Shell[nelem][ion][ns][0] , 
						ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"GAUN") == 0 )
			{
				/* punch gaunt factors */
				if( strcmp(chTime,"LAST") != 0 )
					PunchGaunts(punch.ipPnunit[ipPun]);
			}

			else if( strcmp(punch.chPunch[ipPun],"HISp") == 0 )
			{
				/* punch pressure history of current zone */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* note if pressure convergence faliure occured in history that follows */
					if( !conv.lgConvPres )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"#PROBLEM  Pressure not converged iter %li zone %li density-pressure follows:\n",
							iteration , nzone );
					}
					/* note if temperature convergence faliure occured in history that follows */
					if( !conv.lgConvTemp )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"#PROBLEM  Temperature not converged iter %li zone %li density-pressure follows:\n",
							iteration , nzone );
					}
					for(i=0; i<conv.hist_pres_npres; ++i )
					{
						/* punch history of density - pressure, with correct pressure */
						fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
							iteration,
							nzone,
							conv.hist_pres_density[i],
							conv.hist_pres_current[i],
							conv.hist_pres_correct[i]);
					}
				}

			}

			else if( strcmp(punch.chPunch[ipPun],"HISt") == 0 )
			{
				/* punch temperature history of current zone */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* note if pressure convergence faliure occured in history that follows */
					if( !conv.lgConvPres )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"#PROBLEM  Pressure not converged iter %li zone %li temp heat cool follows:\n",
							iteration , nzone );
					}
					/* note if temperature convergence faliure occured in history that follows */
					if( !conv.lgConvTemp )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"#PROBLEM  Temperature not converged iter %li zone %li temp heat cool follows:\n",
							iteration , nzone );
					}
					for(i=0; i<conv.hist_temp_ntemp; ++i )
					{
						/* punch history of density - pressure, with correct pressure */
						fprintf( punch.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
							iteration,
							nzone,
							conv.hist_temp_temp[i],
							conv.hist_temp_heat[i],
							conv.hist_temp_cool[i]);
					}
				}
			}

			else if( strncmp(punch.chPunch[ipPun],"H2",2) == 0 )
			{
				/* all punch info on large H2 molecule include H2 PDR pdr */
				H2_PunchDo( punch.ipPnunit[ipPun] , punch.chPunch[ipPun] , chTime, ipPun );

			}
			else if( strcmp(punch.chPunch[ipPun],"HEAT") == 0 )
			{
				/* punch heating, routine in file of same name */
				if( strcmp(chTime,"LAST") != 0 )
					HeatPunch(punch.ipPnunit[ipPun]);
			}

			/* punch hummer, results needed for Lya transport, to feed into David's routine */
			else if( strcmp(punch.chPunch[ipPun],"HUMM") == 0 )
			{
				eps = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul/
				  (EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL *dense.eden);
				fprintf( punch.ipPnunit[ipPun], 
					" %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n", 
				  radius.depth_mid_zone, 
				  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn, 
				  iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][1], 
				  iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]*dense.xIonDense[ipHYDROGEN][1], 
				  phycon.te, EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].damp, eps );
			}

			else if( strcmp(punch.chPunch[ipPun],"HYDc") == 0 )
			{
				/* punch hydrogen physical conditions */
				fprintf( punch.ipPnunit[ipPun], 
				  " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
				  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden, 
				  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN], 
				  dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN], 
				  hmi.H2_total/dense.gas_phase[ipHYDROGEN], 
				  hmi.Hmolec[ipMH2p]/dense.gas_phase[ipHYDROGEN], 
				  hmi.Hmolec[ipMH3p]/dense.gas_phase[ipHYDROGEN], 
				  hmi.Hmolec[ipMHm]/dense.gas_phase[ipHYDROGEN] );
			}

			else if( strcmp(punch.chPunch[ipPun],"HYDi") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* punch hydrogen ionization
					 * this will be total decays to ground */
					RateInter = 0.;
					stage = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][0]*dense.eden*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
					fref = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
					fout = iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
					for( ion=ipH2s; ion < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; ion++ )
					{
						/* this is total decays to ground */
						RateInter += 
							EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Aul*
						  (EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Pesc + 
						  EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Pelec_esc + 
						  EmisLines[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Pdest);
						/* total photo from all levels */
						fref += iso.gamnc[ipH_LIKE][ipHYDROGEN][ion]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ion];
						/* total col ion from all levels */
						stage += iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ion]*dense.eden*
						  iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ion];
						fout += iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ion];
					}
					fprintf( punch.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e", 
					  nzone, 
					  /* photo and collis ion rates have units s-1 */
					  iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s], 
					  iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s]* dense.EdenHCorr,
					  ionbal.RateRecomTot[ipHYDROGEN][0] );

					fprintf( punch.ipPnunit[ipPun], "\t%.2e", 
						iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN] );

					fprintf( punch.ipPnunit[ipPun], 
						"\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
					  dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0], 
					  iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/(ionbal.RateRecomTot[ipHYDROGEN][0]), 
					  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][1][ipRecEsc], 
					  RateInter, 
					  fref/MAX2(1e-37,fout), 
					  stage/MAX2(1e-37,fout), 
					  /* simple H+ */
					  iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][0]/(dense.eden* dense.xIonDense[ipHYDROGEN][1]) ,
					  secondaries.csupra[ipHYDROGEN][0]);

					GammaPrt(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0],rfield.nflux,iso.ipOpac[ipH_LIKE][ipHYDROGEN][ipH1s],
					  punch.ipPnunit[ipPun],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s],iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*
					  0.05);
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"HYDp") == 0 )
			{
				/* punch hydrogen populations */
				fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e", 
				  radius.depth_mid_zone, 
				  dense.xIonDense[ipHYDROGEN][0], 
				  dense.xIonDense[ipHYDROGEN][1] );
				for( j=ipH1s; j < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]-1; j++ )
				{
					fprintf( punch.ipPnunit[ipPun], "\t%.2e", 
						iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][j]*dense.xIonDense[ipHYDROGEN][1] );
				}
				fprintf( punch.ipPnunit[ipPun], "\n" );
			}

			else if( strcmp(punch.chPunch[ipPun],"HYD2") == 0 )
			{
				/* punch hydrogen 21 cm optical depth and spin temperature */
				fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\n", 
				radius.depth_mid_zone, HFLines[0].TauIn , phycon.te );
			}

			else if( strcmp(punch.chPunch[ipPun],"HYDl") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch hydrogen line intensities and optical depths */
					for( ipHi=4; ipHi<iso.numLevels_max[ipH_LIKE][ipHYDROGEN]-1 ; ++ipHi )
					{
						for( ipLo=ipHi-5; ipLo<ipHi; ++ipLo )
						{
							if( ipLo<0 )
								continue;
							fprintf(punch.ipPnunit[ipPun], "%li\t%li\t%.4e\t%.2e\n",
								ipHi, ipLo,
								/* wl in cm */
								1./EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].EnergyWN,
								EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].TauIn );

						}
					}
				}
			}

			/* punch lyman alpha - some details about Lya */
			else if( strcmp(punch.chPunch[ipPun],"HYDL") == 0 )
			{
				/* the population ratio for Lya */
				pop = iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]/SDIV(iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]);
				/* the excitation temperature of Lya */
				texc = TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );
				fprintf( punch.ipPnunit[ipPun], 
				  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
				  radius.depth_mid_zone,
				  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn, 
				  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot, 
				  pop, 
				  texc, 
				  phycon.te, 
				  texc/phycon.te ,
				  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc, 
				  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest, 
				  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].pump, 
				  opac.opacity_abs[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1],
				  opac.albedo[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] );
			}

			else if( strcmp(punch.chPunch[ipPun],"HYDr") == 0 )
			{
				/* punch hydrogen recomb cooling for agn */
				phycon.te = 2500.;
				tfidle(FALSE);
				while( phycon.te <= 20000. )
				{
					double r1;
					double ThinCoolingCaseB; 

					tfidle(FALSE);

					r1 = HydroRecCool(1,0);
					ThinCoolingCaseB = pow(10.,((-25.859117 + 
					  0.16229407*phycon.telogn[0] + 
					  0.34912863*phycon.telogn[1] - 
					  0.10615964*phycon.telogn[2])/(1. + 
					  0.050866793*phycon.telogn[0] - 
					  0.014118924*phycon.telogn[1] + 
					  0.0044980897*phycon.telogn[2] + 
					  6.0969594e-5*phycon.telogn[3])))/phycon.te;

					fprintf( punch.ipPnunit[ipPun], " %10.2e\t", 
						phycon.te);
					fprintf( punch.ipPnunit[ipPun], " %10.2e\t", 
						(r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) );

					fprintf( punch.ipPnunit[ipPun], " %10.2e\t", 
						r1/(BOLTZMANN*phycon.te));

					fprintf( punch.ipPnunit[ipPun], " %10.2e\n", 
						ThinCoolingCaseB/(BOLTZMANN*phycon.te));

					phycon.te = phycon.te *2.f;
					tfidle(FALSE);
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"IONI") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch mean ionization distribution */
					PrtMeanIon( 'i', FALSE , punch.ipPnunit[ipPun] );
				}
			}

			/* punch ionization rates */
			else if( strcmp(punch.chPunch[ipPun],"IONR") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* this is element number */
					nelem = (long)punch.punarg[0][ipPun];
					fprintf( punch.ipPnunit[ipPun], 
						"%.5e\t%.4e\t%.4e", 
						radius.depth_mid_zone,
						dense.eden ,
						dynamics.Rate);
					/* >>chng 04 oct 15, from nelem+2 to nelem+1 - array overread -
					 * caught by PnH */
					for( ion=0; ion<nelem+1; ++ion )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"\t%.4e\t%.4e\t%.4e\t%.4e", 
							dense.xIonDense[nelem][ion] ,
							ionbal.RateIonizTot[nelem][ion] ,
							ionbal.RateRecomTot[nelem][ion] ,
							dynamics.Source[nelem][ion] );
					}
					fprintf( punch.ipPnunit[ipPun], "\n");
				}
			}

			else if( strcmp(punch.chPunch[ipPun]," IP ") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch valence shell ip's */
					for( nelem=0; nelem < LIMELM; nelem++ )
					{
						int ion_big;
						double energy;

						/* this is the largest number of ion stages per line */
#						define NELEM_LINE	10
						/* this loop in case all ions do not fit across page */
						for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE )
						{
							int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem);

							/* new line then element name */
							fprintf( punch.ipPnunit[ipPun], 
								"\n%2.2s", elementnames.chElementSym[nelem]);

							/* print ion stages across line */
							for( ion=ion_big; ion <= ion_limit; ++ion )
							{
								fprintf( punch.ipPnunit[ipPun], "\t%4ld", ion+1 );
							}
							fprintf( punch.ipPnunit[ipPun], "\n" );

							/* this loop is over all shells */
							ASSERT( ion_limit < LIMELM );
							/* upper limit is number of shells in atom */
							for( ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ )
							{

								/* print shell label */
								fprintf( punch.ipPnunit[ipPun], "%2.2s", Heavy.chShell[ips]);

								/* loop over possible ions */
								for( ion=ion_big; ion<=ion_limit; ++ion )
								{

									/* does this subshell exist for this ion? break if it does not*/
									/*if( Heavy.nsShells[nelem][ion]<Heavy.nsShells[nelem][0] )*/
									if( ips >= Heavy.nsShells[nelem][ion] )
										break;

									/* array elements are shell, numb of electrons, element, 0 */
									energy = PH1COM.PH1[ips][nelem-ion][nelem][0];

									/* now print threshold with correct format */
									if( energy < 10. )
									{
										fprintf( punch.ipPnunit[ipPun], "\t%6.3f", energy );
									}
									else if( energy < 100. )
									{
										fprintf( punch.ipPnunit[ipPun], "\t%6.2f", energy );
									}
									else if( energy < 1000. )
									{
										fprintf( punch.ipPnunit[ipPun], "\t%6.1f", energy );
									}
									else
									{
										fprintf( punch.ipPnunit[ipPun], "\t%6ld",  (long)(energy) );
									}
								}

								/* put cs at end of long line */
								fprintf( punch.ipPnunit[ipPun], "\n" );
							}
						}
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"LINC") == 0 )
			{
				/* punch line cumulative */
				/* lgOK not used, placdholder */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* not used for anything here, but should not pass baloney*/
					lgOK = TRUE;
					punch_line(punch.ipPnunit[ipPun],"PUNC",lgOK); 
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"LIND") == 0 )
			{
				/* punch line data, then stop */
				PunchLineData(punch.ipPnunit[ipPun]);
				/* above will exit so can't get down to here, but some compilers don't know that 
				puts( "[Normal stop in PunchDo]" );
				fprintf( ioQQQ, " Cloudy ends:\n" );
				cdEXIT(EXIT_FAILURE);*/
			}

			else if( strcmp(punch.chPunch[ipPun],"LINL") == 0 )
			{
				/* punch line labels */
				prt_LineLabels(punch.ipPnunit[ipPun]);
			}

			else if( strcmp(punch.chPunch[ipPun],"LINO") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch line optical depths, routine is below, file static */
					PunchLineStuff(punch.ipPnunit[ipPun],"optical" , punch.punarg[0][ipPun]);
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"LINP") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					static int lgFirst=TRUE;
					/* punch line populations, need to do this twice if very first
					 * call since first call to PunchLineStuff generates atomic parameters
					 * rather than level pops, routine is below, file static */
					PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[0][ipPun]);
					if( lgFirst )
					{
						lgFirst = FALSE;
						PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[0][ipPun]);
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"LINS") == 0 )
			{
				/* punch line structure */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* not acutally used, but should not pass baloney*/
					lgOK = TRUE;
					punch_line(punch.ipPnunit[ipPun],"PUNS",lgOK);
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"LINR") == 0 )
			{
				/* punch line RT */
				if( strcmp(chTime,"LAST") != 0 )
				{
					Punch_Line_RT( punch.ipPnunit[ipPun] , "PUNC" );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"LINA") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					/* punch out all lines with energies */
					for( j=0; j < LineSave.nsum; j++ )
					{
						/* xLineEnergy is sentinel saying to add to continuum only if positive */
						if( LineSv[j].xLineEnergy > 0. && 
							LineSv[j].sumlin > 0. )
						{
							fprintf( punch.ipPnunit[ipPun], "%12.5e\t%8.3f\t%4.4s ", 
							  AnuUnit(LineSv[j].xLineEnergy), 
							  log10(SDIV(LineSv[j].sumlin) ) + radius.Conv2PrtInten, 
							  LineSv[j].chALab );

							prt_wl( punch.ipPnunit[ipPun], LineSv[j].wavelength );

							fprintf( punch.ipPnunit[ipPun], " \t%c\n", 
							  LineSv[j].chSumTyp);
						}
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"LINI") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 && 
					(nzone/punch.LinEvery)*punch.LinEvery != nzone )
				{
					/* this is the last zone
					 * punch line intensities - but do not do last zone twice */
					PunLineIntensity(punch.ipPnunit[ipPun]);
				}
				else if( strcmp(chTime,"LAST") != 0 )
				{
					/* following so we only punch first zone if LinEvery reset */
					if( (punch.lgLinEvery && nzone == 1) || 
					  (nzone/punch.LinEvery)*punch.LinEvery == 
					  nzone )
					{
						/* this is middle of calculation
						 * punch line intensities */
						PunLineIntensity(punch.ipPnunit[ipPun]);
					}
				}
			}

			else if( strcmp( punch.chPunch[ipPun],"LEIL") == 0)
			{
				/* some line intensities for the Leiden PDR,
				 * but only do this when calculation is complete */
				if( strcmp(chTime,"LAST") == 0 )
				{
					double absval , rel ;
					long int n;
					/* the lines we will find,
					 * for a sample list of PDR lines look at LineList_PDR_H2.dat
					 * in the cloudy data dir */
					/* the number of H2 lines */
#					define NLINE_H2 31 
					/* the number of lines which are not H2 */
#					define NLINE_NOTH_H2 5 
					/* the labels and wavelengths for the lines that are not H2 */
					char chLabel[NLINE_NOTH_H2][5]=
					{"C  2", "O  1", "O  1","C  1", "C  1" };
					double Wl[NLINE_NOTH_H2]=
					{157.6 , 63.17 , 145.5 ,609.2 , 369.7 ,  };
					/* these are wavelengths in microns, conv to Ang before call */
					/* >>chng 05 sep 06, many of following wavelengths updated to agree
					 * with output - apparently not updated when energies changed */
					double Wl_H2[NLINE_H2]=
					{2.121 ,
					 28.21 , 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 ,	
					 4.408 , 4.180 , 3.996 , 3.845 , 3.724 , 3.626 , 3.547 , 3.485 , 3.437 , 3.403 ,
					 3.381 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603};
					/* print a header for the lines */
					for( n=0; n<NLINE_NOTH_H2; ++n )
					{
						fprintf(punch.ipPnunit[ipPun], "%s\t%.2f",chLabel[n] , Wl[n]);
						/* get the line, non pos return says didn't find it */
						/* args are 4-char label, wavelength, return log total inten, linear rel inten */
						if( cdLine( chLabel[n] , (float)(Wl[n]*1e4) , &absval , &rel ) <= 0 )
						{
							fprintf(punch.ipPnunit[ipPun], " did not find\n");
						}
						else
						{
							fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
						}
					}
					fprintf(punch.ipPnunit[ipPun], "\n\n\n");

					/* only print the H2 lines if the big molecule is turned on */
					if( h2.lgH2ON )
					{
						fprintf(punch.ipPnunit[ipPun], 
							"Here are some of the H2 Intensities, The first one is the\n"
							"1-0 S(0) line and the following ones are the 0-0 S(X)\n"
							"lines where X goes from 0 to 29\n\n");
						for( n=0; n<NLINE_H2; ++n )
						{
							fprintf(punch.ipPnunit[ipPun], "%s\t%.3f","H2  " , Wl_H2[n]);
							/* get the line, non pos return says didn't find it */
							if( cdLine( "H2  " , (float)(Wl_H2[n]*1e4) , &absval , &rel ) <= 0 )
							{
								fprintf(punch.ipPnunit[ipPun], " did not find\n");
							}
							else
							{
								fprintf(punch.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
							}
						}
					}
#					undef NLINE_H2 
#					undef NLINE_NOTH_H2 
				}

			}
			else if( strcmp( punch.chPunch[ipPun],"LEIS") == 0)
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* get some column densities we shall need */
					double col_ci , col_oi , col_cii, col_heii;
					if( cdColm("carb" , 1 , &col_ci ) )
						TotalInsanity();
					if( cdColm("carb" , 2 , &col_cii ) )
						TotalInsanity();
					if( cdColm("oxyg" , 1 , &col_oi ) )
						TotalInsanity();
					if( cdColm("heli" , 2 , &col_heii ) )
						TotalInsanity();
					/* punch leiden structure - some numbers for the Leiden PDR model comparisons */
					fprintf( punch.ipPnunit[ipPun], 
					"%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
					"%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
					"%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t" 
					"%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
					"%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
					/* depth of this point */
					radius.depth_mid_zone,
					/* A_V for an extended source */
					0.00,	
					/* A_V for a point source */
					rfield.extin_mag_V_point,
					/* temperature */
					phycon.te ,
					dense.xIonDense[ipHYDROGEN][0],
					hmi.H2_total,
					dense.xIonDense[ipCARBON][0],
					dense.xIonDense[ipCARBON][1],
					dense.xIonDense[ipOXYGEN][0],
					co.hevmol[ipCO],
					co.hevmol[ipO2],
					co.hevmol[ipCH],
					co.hevmol[ipOH],
					dense.eden,
					dense.xIonDense[ipHELIUM][1],
					dense.xIonDense[ipHYDROGEN][1],
					hmi.Hmolec[ipMH3p],
					colden.colden[ipCOL_H0],
					colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
					col_ci,
					col_cii,
					col_oi,
					co.hevcol[ipCO],
					co.hevcol[ipO2],
					co.hevcol[ipCH],
					co.hevcol[ipOH],
					colden.colden[ipCOL_elec],
					col_heii,
					colden.colden[ipCOL_Hp],
					colden.colden[ipCOL_H3p],
					hmi.H2_Solomon_dissoc_rate_used_H2g ,
					gv.rate_h2_form_grains_used_total,
					hmi.UV_Cont_rel2_Draine_DB96_depth,
					/* CO and C dissociation rate */
				       	co.CO_photo_dissoc_rate,
					/* total CI ionization rate */
					ionbal.PhotoRate_Shell[ipCARBON][0][2][0],
					/* total heating, erg cm-3 s-1 */
					thermal.htot,
					/* total cooling, erg cm-3 s-1 */
					thermal.ctot,
					/* GrnP grain photo heating */
					thermal.heating[0][13],
					/* grain collisional cooling */
					MAX2(0,gv.GasCoolColl),	
					/* grain collisional heating */
					-1*MIN2(0,gv.GasCoolColl),	
					/* COds - CO dissociation heating */
					thermal.heating[0][9],
					/* H2dH-Heating due to H2 dissociation */
					hmi.HeatH2Dish_used,
					/* H2vH-Heating due to collisions with H2 */
					hmi.HeatH2Dexc_used ,
					/* ChaT - charge transfer heating */
					thermal.heating[0][24] ,
					/* cosmic ray heating */
					thermal.heating[1][6] ,
					/* heating due to atoms of various heavy elements */
					thermal.heating[ipMAGNESIUM][0],
					thermal.heating[ipSULPHUR][0],
					thermal.heating[ipSILICON][0],
					thermal.heating[ipIRON][0],
					thermal.heating[ipSODIUM][0],
					thermal.heating[ipALUMINIUM][0],
					thermal.heating[ipCARBON][0],
					TauLines[ipT610].cool,
					TauLines[ipT370].cool,
					TauLines[ipT157].cool,
					TauLines[ipT63].cool,
					TauLines[ipT146].cool );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"MAP ") == 0 )
			{
				/* do the map now if we are at the zone, or if this
				 * is the LAST call to this routine and map not done yet */
				if(  (map.lgMapDone==FALSE ) &&
					(nzone == map.MapZone  ||  strcmp(chTime,"LAST") == 0 ) )
				{
					lgTlkSav = called.lgTalk;
					called.lgTalk = TRUE;
					map.lgMapBeingDone = TRUE;
					map_do(punch.ipPnunit[ipPun] , " map");
					called.lgTalk = lgTlkSav;
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"MOLE") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* molecules, especially for pdr, first give radius */
					fprintf( punch.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone );

					/* visual extinction of point source (star)*/
					fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);

					/* visual extinction of an extended source (like a PDR)*/
					fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);

					/* carbon monoxide photodissociation rate */
					fprintf( punch.ipPnunit[ipPun], "%.5e\t" , co.CO_photo_dissoc_rate );

					/* carbon recombination rate */
					fprintf( punch.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] );

					/* now do all the H molecules */
					for(j=0; j<N_H_MOLEC; ++j )
					{
						fprintf(punch.ipPnunit[ipPun],"\t%.2e",hmi.Hmolec[j] );
					}

					/* now all the C molecules */
					for(j=0; j<NUM_COMOLE_CALC; ++j )
					{
						fprintf(punch.ipPnunit[ipPun],"\t%.2e",co.hevmol[j] );
					}
					fprintf(punch.ipPnunit[ipPun],"\n");
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"OPAC") == 0 )
			{
				/* punch opacity - routine will parse which type of opacity punch to do */
				if( strcmp(chTime,"LAST") == 0 )
					punch_opacity(punch.ipPnunit[ipPun],ipPun);
			}

			/* punch coarse optical depths command */
			else if( strcmp(punch.chPunch[ipPun],"OPTc") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					for( j=0; j < rfield.nflux; j++ )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"%12.4e\t%.3e\t%12.4e\t%.3e\n", 
						  AnuUnit(rfield.AnuOrg[j]), 
						  opac.TauAbsFace[j]+opac.TauScatFace[j], 
						  opac.TauAbsFace[j], 
						  opac.TauScatFace[j] );
					}
				}
			}

			/* punch fine optical depths command */
			else if( strcmp(punch.chPunch[ipPun],"OPTf") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					long nu_hi , nskip;
					if( punch.punarg[0][ipPun] > 0. )
					{
						/* possible lower bounds to energy range - will be zero if not set */
						j = ipFineCont( punch.punarg[0][ipPun] );
					}
					else
					{
						j = 0;
					}

					/* upper limit */
					if( punch.punarg[1][ipPun]> 0. )
					{
						nu_hi = ipFineCont( punch.punarg[1][ipPun]);
					}
					else
					{
						nu_hi = rfield.nfine;
					}
					nskip = (long)punch.punarg[2][ipPun];
					do
					{
						float sum1=rfield.fine_opt_depth[j];
						float sum2=rfield.fine_opac[j];
						float xnu = rfield.fine_anu[j+nskip/2];
						for( jj=1; jj<nskip; ++jj )
						{
							sum1 += rfield.fine_opt_depth[j+jj];
							sum2 += rfield.fine_opac[j+jj];
						}
						if( sum2 > 0. )
							fprintf( punch.ipPnunit[ipPun], 
							  "%12.6e\t%.3e\t%.3e\n", 
							  AnuUnit(xnu), 
							  sum1/nskip , 
							  sum2/nskip);
						j += nskip;
					}while( j < nu_hi );
				}
			}

			else if( strcmp(punch.chPunch[ipPun]," OTS") == 0 )
			{
				ConMax = 0.;
				xLinMax = 0.;
				opConSum = 0.;
				opLinSum = 0.;
				ipLinMax = 1;
				ipConMax = 1;

				for( j=0; j < rfield.nflux; j++ )
				{
					opConSum += rfield.otscon[j]*opac.opacity_abs[j];
					opLinSum += rfield.otslin[j]*opac.opacity_abs[j];
					if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax )
					{
						xLinMax = rfield.otslin[j]*opac.opacity_abs[j];
						ipLinMax = j+1;
					}
					if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax )
					{
						ConMax = rfield.otscon[j]*opac.opacity_abs[j];
						ipConMax = j+1;
					}
				}
				fprintf( punch.ipPnunit[ipPun], 
				  "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n", 
				  opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1]
				  , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1]
				  , rfield.anu[ipConMax-1], ConMax );
			}

			else if( strcmp(punch.chPunch[ipPun],"OVER") == 0 )
			{
				/* punch overview
				 * this is the floor for the smallest ionization fractions printed */
				double toosmall = SMALLFLOAT ,
					hold;

				/* overview of model results,
				 * depth, te, hden, eden, ion fracs H, He, c, O */
				if( strcmp(chTime,"LAST") != 0 )
				{
					
					/* print the depth */
					fprintf( punch.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone );

					/* temperature, heating */
					if (dynamics.Cool > dynamics.Heat) 
					{
						fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f", 
						  log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat)) ); 
					}
					else
					{
						double diff = fabs(thermal.htot-dynamics.Cool);
						fprintf( punch.ipPnunit[ipPun], "%.4f\t%.3f", 
						  log10(phycon.te), log10(SDIV(diff)) ); 
					}

					/* hydrogen and electron densities */
					fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f", 
					  log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) );
					  
					/* molecular fraction of hydrogen */
					fprintf( punch.ipPnunit[ipPun], "\t%.4f", 
					  /*log10(MAX2(toosmall,2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN])) );*/
					  log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) );
					  
					/* ionization fractions of hydrogen */
					fprintf( punch.ipPnunit[ipPun], "\t%.4f\t%.4f", 
					  log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])), 
					  log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );

					/* ionization fractions of helium */
					for( j=1; j <= 3; j++ )
					{
						double arg = SDIV(dense.gas_phase[ipHELIUM]);
						arg = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg );
						fprintf( punch.ipPnunit[ipPun], "\t%.4f", 
						  log10(arg) );
					}

					/* carbon monoxide molecular fraction of CO */
					hold = SDIV(dense.gas_phase[ipCARBON]);
					hold = MAX2(toosmall, co.hevmol[ipCO]/hold);
					fprintf( punch.ipPnunit[ipPun], "\t%.4f", 
					  log10(hold) );

					/* ionization fractions of carbon */
					for( j=1; j <= 4; j++ )
					{
						hold = SDIV(dense.gas_phase[ipCARBON]);
						hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold);
						fprintf( punch.ipPnunit[ipPun], "\t%.4f", 
						  log10(hold) );
					}

					/* ionization fractions of oxygen */
					for( j=1; j <= 6; j++ )
					{
						hold = SDIV(dense.gas_phase[ipOXYGEN]);
						hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold);
						fprintf( punch.ipPnunit[ipPun], "\t%.4f", 
						  log10(hold) );
					}

					/* visual extinction of point source (star)*/
					fprintf( punch.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point);

					/* visual extinction of an extended source (like a PDR)*/
					fprintf( punch.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended);
				}
			}

			else if( strcmp(punch.chPunch[ipPun]," PDR") == 0 )
			{
				/* this is the punch pdr command */
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* convert optical depth at wavelength of V filter
					 * into magnitudes of extinction */
					/* >>chyng 03 feb 25, report extinction to illuminated face,
					 * rather than total extinction which included far side when
					 * sphere was set */
					/*av = opac.TauTotalGeo[0][rfield.ipV_filter-1]*1.08574;*/

					fprintf( punch.ipPnunit[ipPun], 
						"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t", 
					  radius.depth_mid_zone, 
					  /* total hydrogen column density, all forms */
					  colden.colden[ipCOL_HTOT], 
					  phycon.te, 
					  /* fraction of H that is atomic */
					  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN], 
					  /* ratio of n(H2) to total H, == 0.5 when fully molecular */
					  2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN], 
					  2.*hmi.Hmolec[ipMH2s]/dense.gas_phase[ipHYDROGEN], 
					  /* atomic to total carbon */
					  dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]), 
					  co.hevmol[ipCO]/SDIV(dense.gas_phase[ipCARBON]), 
					  co.hevmol[ipH2O]/SDIV(dense.gas_phase[ipOXYGEN]),
					  /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */
					  hmi.UV_Cont_rel2_Habing_TH85_depth);

					/* visual extinction due to dust alone, of point source (star)*/
					fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);

					/* visual extinction due to dust alone,  of an extended source (like a PDR)*/
					fprintf( punch.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);

					/* visual extinction (all sources) of a point source (like a PDR)*/
					fprintf( punch.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"PHYS") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					/* punch physical conditions */
					fprintf( punch.ipPnunit[ipPun], "%.5e%12.4e%11.3e%11.3e%11.3e%11.3e\n", 
					  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], 
					  dense.eden, thermal.htot, wind.AccelTot );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"PRES") == 0 )
			{
				/* the punch pressure command */
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], 
						"%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%c\n", 
					  radius.depth_mid_zone, 
					  /* Pcorrect */
					  pressure.PresTotlCorrect, 
					  /* Pcurrent */
					  pressure.PresTotlCurr, 
					  /* Pln + pintg */
					  pressure.PresTotlInit + pressure.PresInteg, 
					  /* pgas (0) */
					  pressure.PresTotlInit, 
					  /* Pgas */
					  pressure.PresGasCurr, 
					  /* Pram */
					  pressure.PresRamCurr, 
					  /* Prad */
					  pressure.PresRadCurr, 
					  /* subtract continuum rad pres which has already been added on */
					  pressure.PresInteg - pressure.pinzon, 
					  /* wind speed in km/s */
					  wind.windv/1e5,
					  /* sound speed in km/s */
					  timesc.sound_speed_adiabatic/1e5,
					  /* the magnetic pressure */
					  magnetic.pressure ,
					  /* the local turbulent velocity */
					  DoppVel.TurbVel ,
					  /* turbulent pressure */
					  pressure.PresTurbCurr ,
					  TorF(conv.lgConvPres) );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"RADI") == 0 )
			{
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], " %4ld%.5e%12.4e%12.4e\n", 
					  nzone, radius.Radius_mid_zone, radius.depth_mid_zone, 
					  radius.drad );
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"RESU") == 0 )
			{
				/*  punch results of the calculation */
				if( strcmp(chTime,"LAST") == 0 )
					punResults(punch.ipPnunit[ipPun]);
			}

			else if( strcmp(punch.chPunch[ipPun],"SECO") == 0 )
			{
				/*  punch secondary ionization */
				if( strcmp(chTime,"LAST") != 0 )
					fprintf(punch.ipPnunit[ipPun],
					"%.5e\t%.3e\t%.3e\t%.3e\n",
					radius.depth ,
					secondaries.csupra[ipHYDROGEN][0],
					secondaries.csupra[ipHYDROGEN][0]*2.02,
					secondaries.x12tot );
			}

			else if( strcmp(punch.chPunch[ipPun],"SOUS") == 0 )
			{
				/* full spectrum of continuum source function at 1 depth
				 *  command was "punch source spectrum" */
				if( strcmp(chTime,"LAST") != 0 )
				{
					limit = MIN2(rfield.ipMaxBolt,rfield.nflux);
					for( j=0; j < limit; j++ )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"%12.4e%12.4e%12.4e%12.4e%12.4e\n", 
						  rfield.anu[j], 
						  rfield.ConEmitLocal[j]/rfield.widflx[j], 
						  opac.opacity_abs[j], 
						  rfield.ConEmitLocal[j]/rfield.widflx[j]/MAX2(1e-35,opac.opacity_abs[j]), 
						  rfield.ConEmitLocal[j]/rfield.widflx[j]/MAX2(1e-35,opac.opacity_abs[j])/plankf(j) );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"SOUD") == 0 )
			{
				/* parts of continuum source function vs depth
				 * command was source depth */
				j = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] + 2;
				fprintf( punch.ipPnunit[ipPun], 
					"%12.4e%12.4e%12.4e%12.4e\n", 
				  opac.TauAbsFace[j-1], 
				  rfield.ConEmitLocal[j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]), 
				  rfield.otscon[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1], 
				  rfield.otscon[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/opac.opacity_abs[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] );
			}

			/* this is punch special option */
			else if( strcmp(punch.chPunch[ipPun],"SPEC") == 0 )
			{
				PunchSpecial(punch.ipPnunit[ipPun],chTime);
			}

			else if( strcmp(punch.chPunch[ipPun],"TEGR") == 0 )
			{
				/* punch history of heating and cooling */
				if( strcmp(chTime,"LAST") != 0 )
				{
					fprintf( punch.ipPnunit[ipPun], " " );
					for( j=0; j < NGRID; j++ )
					{
						fprintf( punch.ipPnunit[ipPun], 
							"%4ld%11.3e%11.3e%11.3e\n", 
						  thermal.nZonGrid[j], 
						  thermal.TeGrid[j], 
						  thermal.HtGrid[j], 
						  thermal.ClGrid[j] );
					}
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"TEMP") == 0 )
			{
				static double deriv_old=-1;
				double deriv=-1. , deriv_sec;
				/* temperature and its derivatives,
				 * if second zone then have one deriv */
				if( nzone >1 )
				{
					deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
					fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e", 
					radius.depth_mid_zone, 
					phycon.te, 
					deriv );
					/* if third zone then have second deriv */
					if( nzone > 2 )
					{
						deriv_sec = (deriv-deriv_old)/ radius.drad;
						fprintf( punch.ipPnunit[ipPun], "\t%.2e", 
						deriv_sec );
					}
					deriv_old = deriv;
					fprintf( punch.ipPnunit[ipPun], "\n");
				}
			}

			else if( strcmp(punch.chPunch[ipPun],"TIME") == 0 )
			{
				static double ElapsedTime , ZoneTime;
				if( nzone<=1 )
				{
					ElapsedTime = cdExecTime();
					ZoneTime = 0.;
				}
				else
				{
					double t = cdExecTime();
					ZoneTime = t - ElapsedTime;
					ElapsedTime = t;
				}

				/* zone, time for this zone, elapsed time */
				fprintf( punch.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n", 
				  nzone,  ZoneTime , ElapsedTime );
			}

			else if( strcmp(punch.chPunch[ipPun],"TPRE") == 0 )
			{
				/* temperature and its predictors, turned on with punch tprid */
				fprintf( punch.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n", 
				  nzone, phycon.TeInit, phycon.TeProp, phycon.te, 
				  (phycon.TeProp- phycon.te)/phycon.te );
			}

			else if( strcmp(punch.chPunch[ipPun],"WIND") == 0 )
			{
				/* wind velocity, radiative accel, and ratio total to e- accel */
				fprintf( punch.ipPnunit[ipPun], 
					"%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 
				  radius.Radius_mid_zone, 
				  radius.depth_mid_zone, 
				  wind.windv, 
				  wind.AccelTot, 
				  wind.AccelLine,
				  wind.AccelCont ,
				  wind.fmul );
			}
			else if( strcmp(punch.chPunch[ipPun],"XSPA") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					if( grid.lgGrid )
                        punchFITSfile( punch.ipPnunit[ipPun], "XSPA" );
					else
					{
						fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
						puts( "[Stop in PunchDo]" );
                        cdEXIT(EXIT_FAILURE);
					}
				}
			}
			else if( strcmp(punch.chPunch[ipPun],"XSPM") == 0 )
			{
				if( strcmp(chTime,"LAST") == 0 )
				{
					if( grid.lgGrid )
                        punchFITSfile( punch.ipPnunit[ipPun], "XSPM" );
					else
					{
						fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
						puts( "[Stop in PunchDo]" );
                        cdEXIT(EXIT_FAILURE);
					}
				}
			}
			else
			{
				/* this is insanity, internal flag set in getpunch not seen here */
				fprintf( ioQQQ, " PunchDo does not recognize flag %4.4s set by ParsePunch.  This is impossible.\n", 
				  punch.chPunch[ipPun] );
				ShowMe();
				puts( "[Stop in PunchDo]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* at start of loop, broke to here if not last iteration,
			 * and punch last was set
			 * throw special line if punching every iteration */
			if( strcmp(chTime,"LAST") == 0 && !iterations.lgLastIt && punch.lgHashEndIter )
			{
				/*fprintf( punch.ipPnunit[ipPun], " ###########################\n" );*/
				fprintf( punch.ipPnunit[ipPun], "%s\n",punch.chHashString );
			}
			if( punch.lgFLUSH )
				fflush( punch.ipPnunit[ipPun] );
		}
	}

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

/*PunLineIntensity produce the 'punch lines intensity' output */
static void PunLineIntensity(FILE * io)
{
	char chCard[INPUT_LINE_LENGTH];
	int lgEOF;
	long int i;

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

	/* used to punch out all the emission line intensities
	 * first initialize the line image reader */

	/* start the file with the input commands */
	input_init();
	fprintf( io, "**********************************************************************************************************************************\n" );
	
	lgEOF = FALSE;
	while( !lgEOF )
	{
		input_readarray(chCard,&lgEOF);
		if( !lgEOF )
		{
			fprintf( io, "%s\n", chCard );
		}
	}

	/* now print any cautions or warnings */
	cdWarnings( io);
	cdCautions( io);
	fprintf( io, "zone=%5ld\n", nzone );
	fprintf( io, "**********************************************************************************************************************************\n" );
	fprintf( io, "begin emission lines\n" );

	/* >>chng 96 jul 03, only punch non-zero intensities */
	PunResults1Line(io,"    ",0,0.,"Start");
	for( i=0; i < LineSave.nsum; i++ )
	{
		if( LineSv[i].sumlin > 0. )
		{
			PunResults1Line(io,(char*)LineSv[i].chALab,LineSv[i].wavelength,
			  LineSv[i].sumlin, "Line ");
		}
	}

	PunResults1Line(io,"    ",0,0.,"Flush");

	fprintf( io, "     \n" );
	fprintf( io, "**********************************************************************************************************************************\n" );


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

/* lgPunchOpticalDepths true says punch optical depths */
static int lgPopsFirstCall , lgPunchOpticalDepths ;

/*PunchLineStuff punch optical depths or source functions for all transferred lines */
static void PunchLineStuff(
  FILE * io,
  const char *chJob , 
  float xLimit )
{

	long int nelem , ipLo , ipHi , i , ipISO;
	long index = 0;
	static int lgFirst=TRUE;

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

	/*find out which job this is and set a flag to use later */
	if( strcmp( &*chJob , "optical" ) == 0 )
	{
		/* punch line optical depths */
		lgPunchOpticalDepths = TRUE;
		lgPopsFirstCall = FALSE;
	}
	else if( strcmp( &*chJob , "populat" ) == 0 )
	{
		lgPunchOpticalDepths = FALSE;
		/* level population information */
		if( lgFirst )
		{
			lgPopsFirstCall = TRUE;
			fprintf(io,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
			lgFirst = FALSE;
		}
		else
		{
			lgPopsFirstCall = FALSE;
		}
	}
	else
	{
		fprintf( ioQQQ, " insane job in PunchLineStuff =%s\n", 
		  &*chJob );
		puts( "[Stop in PunchLineStuff]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* loop over all lines, calling put1Line to create info (routine located below) */
	/* hydrogen like lines */
	/* >>chng 02 may 16, had been explicit H and He-like loops */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem]  )
			{
				for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
				{
					for( ipLo=0; ipLo <ipHi; ipLo++ )
					{
						++index;
						pun1Line( &EmisLines[ipISO][nelem][ipHi][ipLo] , io , xLimit  , index , dense.xIonDense[nelem][nelem+1-ipISO] );
					}
				}
				/* also do extra Lyman lines if optical depths are to be done,
				 * these are line that are included only for absorption, not in the
				 * model atoms */
				if( lgPunchOpticalDepths )
				{
					/* >>chng 02 aug 23, for he-like, had starting on much too high a level since
					 * index was number of levels - caught by Adrian Turner */
					/* now output extra line lines, starting one above those already done above */
					/*for( ipHi=iso.numLevels_max[ipISO][nelem]; ipHi < iso.nLyman[ipISO]; ipHi++ )*/
					for( ipHi=iso.quant_desig[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n+1; ipHi < iso.nLyman[ipISO]; ipHi++ )
					{
						++index;
						pun1Line( &iso.ExtraLymanLines[ipISO][nelem][ipHi] , io , xLimit  , index , dense.xIonDense[nelem][nelem+1-ipISO]);
					}
				}
			}
		}
	}

	/* index from 1 to skip over dummy line */
	for( i=1; i < nLevel1; i++ )
	{
		++index;
		pun1Line( &TauLines[i] , io , xLimit  , index , 1. );
	}

	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			++index;
			pun1Line( &TauLine2[i] , io , xLimit  , index , 1. );
		}
	}

	for( i=0; i < nUTA; i++ )
	{
		++index;
		pun1Line( &UTALines[i] , io , xLimit  , index , 1. );
	}


	/* do optical depths of FeII lines, if large atom is turned on */
	FeIIPunchLineStuff( io , xLimit  , index);

	/* punch optical depths of H2 lines */
	H2_PunchLineStuff( io , xLimit  , index);

	/* C12O16 */
	for( i=0; i < nCORotate; i++ )
	{
		++index;
		pun1Line( &C12O16Rotate[i] , io , xLimit  , index , 1. );
	}

	/* C13O16 */
	for( i=0; i < nCORotate; i++ )
	{
		++index;
		pun1Line( &C13O16Rotate[i] , io , xLimit  , index , 1. );
	}
	/*fprintf(io, "##################################\n"); */
	fprintf( io , "%s\n",punch.chHashString );

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

/*pun1Line called by PunchLineStuff to produce output for one line */
void pun1Line( EmLine * t , FILE * io , float xLimit  , long index , double abundance)
{

	if( lgPunchOpticalDepths )
	{
		/* optical depths, no special first time, only print them */
		if( t->TauIn >= xLimit )
		{
			/* label like "C  4" or "H  1" */
			fprintf( io , "%2.2s%2.2s\t", 
			  elementnames.chElementSym[t->nelem-1] ,
			  elementnames.chIonStage[t->IonStg-1]  );

			/* print wavelengths, either line in main printout labels, 
			 * or in various units in exponential notation - prt_wl is in prt.c */
			if( strcmp( punch.chConPunEnr[punch.ipConPun], "labl" )== 0 )
			{
				prt_wl( io , t->WLAng );
			}
			else
			{
				/* this converts energy in rydbergs into any of the other units */
				fprintf( io , "%.7e", AnuUnit((float)(t->EnergyWN * WAVNRYD)) );
			}
			/* print the optical depth */
			fprintf( io , "\t%.3f", t->TauIn  );
			/* damping constant */
			fprintf(io, "\t%.3e", 
				t->dampXvel / DoppVel.doppler[ t->nelem -1 ] );
			fprintf(io, "\n");
		}
	}
	else if( lgPopsFirstCall )
	{
		char chLbl[11];
		/* first call to line populations, print atomic parameters and indices */
		strcpy( chLbl, chLineLbl(t) );
		fprintf(io, "%li\t%s" , index , chLbl );
		/* stat weights */
		fprintf(io, "\t%.0f\t%.0f", 
			t->gLo ,t->gHi);
		/* energy difference, gf */
		fprintf(io, "\t%.2f\t%.3e", 
			t->EnergyWN ,t->gf);
		fprintf(io, "\n");
	}
	else
	{
		/* not first call, so do level populations and indices defined above */
		if( t->PopHi > xLimit )
		{
			fprintf(io,"%li\t%.2e\t%.2e\n", index,
				/* >>chng 05 may 08, add abundances, which for iso-seq species is
				 * the density of the parent ion, for other lines, is unity.
				 * had not been included so pops for iso seq were rel to parent ion.
				 * caught by John Evrett */
				t->PopLo*abundance , 
				t->PopHi*abundance );
		}
	}
}

/*PunchNewContinuum produce the 'punch new continuum' output */
static void PunchNewContinuum(FILE * io, long ipCon )
{
	long int ipLo, ipHi,
		j ,
		ncells;

	double 
		wllo=3500. , 
		wlhi=7000. , 
		resolution = 2.;

	double *energy, 
		*cont_incid,
		*cont_atten,
		*diffuse_in,
		*diffuse_out,
		*emis_lines_out,
		*emis_lines_in;

	/* get the low limit, cp_range_min is zero if not set */
	wllo = MAX2( rfield.anu[0] , punch.cp_range_min[ipCon] );
	if( punch.cp_range_max[ipCon] > 0. )
	{
		/* set to smaller of these two  */
		wlhi = MIN2( rfield.anu[rfield.nflux-1] , punch.cp_range_max[ipCon]  );
	}
	else
	{
		/* use high energy limit since not set */
		wlhi = rfield.anu[rfield.nflux-1] ;
	}

	if( punch.cp_resolving_power[ipCon] != 0. )
	{
		/* next two bogus to fool lint - they cannot happen */
		ipLo = LONG_MAX;
		ipHi = LONG_MAX;
		/* wwe will set a new continuum mesh */
		ncells = (long)((wlhi-wllo)/resolution + 1);
	}
	else
	{
		/* use native continuum mesh */
		ipLo = ipoint(wllo)-1;
		ipHi = ipoint(wlhi)-1;
		ncells = ipHi - ipLo + 1;
	}

	/* now allocate the space */
	if( (energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	if( (cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	if( (cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	if( (diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	if( (diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	if( (emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	if( (emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	/*if( (emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();
	if( (emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) ) )==NULL )
		BadMalloc();*/

	/* was the resolution changed from the default native resolution ? */
	if( punch.cp_resolving_power[ipCon] != 0. )
	{
		/* now create energy grid */
		energy[0] = wlhi;
		j = 0;
		while( energy[j]-resolution > wllo )
		{
			++j;
			ASSERT( j< ncells+1 );
			energy[j] = energy[j-1] - resolution;
		}

		ncells = j;
		/* now convert to rydbergs */
		for( j=0; j<ncells; ++j )
		{
			energy[j] = RYDLAM / energy[j];
		}
	}
	else 
	{
		/* now convert to rydbergs */
		for( j=0; j<ncells; ++j )
		{
			energy[j] = rfield.AnuOrg[j+ipLo] - rfield.widflx[j+ipLo]/2.;
		}
	}

	/* for cdSPEC the energy vector is the lower edge of the energy cell */
	/* get incident continuum */
	cdSPEC( 1 , energy , ncells , cont_incid );
	/* get attenuated incident continuum */
	cdSPEC( 2 , energy , ncells , cont_atten );
	/* get attenuated incident continuum */
	cdSPEC( 5 , energy , ncells , diffuse_in );
	/* get attenuated incident continuum */
	cdSPEC( 4 , energy , ncells , diffuse_out );
	/*TODO	2	- NB - if continuum resolution changed the lines WILL NOT WORK */
	/* different parts of lines */
	cdSPEC( 6 , energy , ncells , emis_lines_out );
	cdSPEC( 7 , energy , ncells , emis_lines_in );
	/*cdSPEC( 8 , energy , ncells , emis_lines_pump_out );
	cdSPEC( 9 , energy , ncells , emis_lines_pump_in );*/

	/* for this example we will do a wavelength range */
	for( j=0; j<ncells-1 ; ++j )
	{
		/* photon energy in appropriate energy or wavelength units */
		fprintf( io,"%.3e\t", AnuUnit((float)(energy[j]+rfield.widflx[j+ipLo]/2.) ) );
		fprintf( io,"%.3e\t", cont_incid[j] );
		fprintf( io,"%.3e\t", cont_atten[j] );
		fprintf( io,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
		fprintf( io,"%.3e", 
			emis_lines_out[j]+emis_lines_in[j]/*+emis_lines_pump_out[j]+emis_lines_pump_in[j]*/ );
		fprintf( io,"\n" );
	}

	free(energy);
	free(cont_incid);
	free(diffuse_in);
	free(diffuse_out);
	free(cont_atten);
	free(emis_lines_out);
	free(emis_lines_in);
	/*free(emis_lines_pump_out);
	free(emis_lines_pump_in );*/
}

/* punch agn hemiss, for chapt 4, routine is below */
static void AGN_Hemis(FILE *io )
{
#	define NTE 4
	float te[NTE] = {5000., 10000., 15000., 20000.};
	float *continuum[NTE];
	double TempSave = phycon.te;
	long i , j;

	/* make table of continuous emission at various temperatuers */
	/* first allocate space */
	for( i=0;i<NTE; ++i)
	{
		if( (continuum[i] = (float *)MALLOC((unsigned)rfield.nflux*sizeof(float) ))==NULL)
			BadMalloc();

		/* set the next temperature */
		phycon.te = te[i];
		/* recompute everything at this new temp */
		tfidle(TRUE);
		/* converge the pressure-temperature-ionization solution for this zone */
		/*lint -e534 ignore return value of function */
		ConvPresTempEdenIoniz();
		/*lint +e534 ignore return value of function */

		/* now get the thermal emission */
		RT_diffuse();
		for(j=0;j<rfield.nflux; ++j )
		{
			continuum[i][j] = rfield.ConEmitLocal[j]/(float)dense.eden/
				(dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] );
		}
	}

	/* print title for line */
	fprintf(io,"wl");
	for( i=0;i<NTE; ++i)
	{
		fprintf(io,"\t%.0f",te[i]);
	}
	fprintf( io , "\tcont\n"); 

	/* not print all n temperatures across a line */
	for(j=0;j<rfield.nflux; ++j )
	{
		fprintf( io , "%12.5e", 
		  AnuUnit(rfield.AnuOrg[j]) );
		/* loop over the temperatures, and for each, calculate a continuum */
		for( i=0;i<NTE; ++i)
		{
			fprintf(io,"\t%.3e",continuum[i][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
		}
		/* cont label and end of line*/
		fprintf( io , "\t%s\n" , rfield.chContLabel[j]); 
	}

	/* now free the continua */
	for( i=0;i<NTE; ++i)
	{
		free( continuum[i] );
	}

	/* Restore temperature stored before this routine was called	*/
	phycon.te = TempSave;
	/* and force update */
	tfidle(TRUE);

	fprintf( ioQQQ, " I have left the code in a disturbed state, and will now exit.\n");
	puts( "[Stop in PunchDo]" );
	cdEXIT(EXIT_FAILURE);
#	undef NTE
}

/*punResults punch results from punch results command */
/*PunResults1Line do single line of output for the punch results and punch line intensity commands */
static void punResults(FILE* io)
{
	char chCard[INPUT_LINE_LENGTH];
	int lgEOF;
	long int i , nelem , ion;

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

	/* used to punch out line intensities, optical depths,
	 * and column densities */
	lgEOF = FALSE;
	/* first initialize the line image reader */
	input_init();

	fprintf( io, "**********************************************************************************************************************************\n" );
	lgEOF = FALSE;
	input_readarray(chCard,&lgEOF);
	while( !lgEOF )
	{
		fprintf( io, "%s\n", chCard );
		input_readarray(chCard,&lgEOF);
	}

	/* first print any cautions or warnings */
	cdWarnings(io);
	cdCautions(io);
	fprintf( io, "**********************************************************************************************************************************\n" );

	fprintf( io, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );

	fprintf( io, "BEGIN EMISSION LINES\n" );
	PunResults1Line(io,"    ",0,0.,"Start");

	for( i=0; i < LineSave.nsum; i++ )
	{
		if( LineSv[i].sumlin > 0. )
		{
			PunResults1Line(io,(char*)LineSv[i].chALab,LineSv[i].wavelength,
			  LineSv[i].sumlin,"Line ");
		}
	}

	PunResults1Line(io,"    ",0,0.,"Flush");

	fprintf( io, "     \n" );

	fprintf( io, "BEGIN COLUMN DENSITIES\n" );

	/* this dumps out the whole array,*/
	/* following loop relies on LIMELM being 30, assert it here in case
	 * this is ever changed */
	ASSERT( LIMELM == 30 );
	/* this order of indices is to keep 30 as the fastest variable,
	 * and the 32 (LIMELM+1) as the slower one */
	for( nelem=0; nelem<LIMELM; nelem++ )
	{
		for(ion=0; ion < nelem+1; ion++)
		{
			fprintf( io, " %10.3e", mean.xIonMeans[0][nelem][ion] );
			/* throw line feed every 10 numbers */
			if( nelem==9|| nelem==19 || nelem==29 )
			{
				fprintf( io, "\n" );
			}
		}
	}

	fprintf( io, "END OF RESULTS\n" );
	fprintf( io, "**********************************************************************************************************************************\n" );

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

#define	LINEWIDTH	6

/*PunResults1Line do single line of output for the punch results and punch line intensity commands */
/* the number of emission lines across one line of printout */
static void PunResults1Line(
  FILE * io, 
  /* 4 char + null string */
  const char *chLab, 
  float wl, 
  double xInten, 
  /* null terminated string saying what to do */
  const char *chFunction)
{

	long int i;
	static float wavelength[LINEWIDTH];
	static long ipLine;
	static double xIntenSave[LINEWIDTH];
	static char chLabSave[LINEWIDTH][5];

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

	/* if LineWidth is changed then change format in write too */

	if( strcmp(chFunction,"Start") == 0 )
	{
		ipLine = 0;
	}
	else if( strcmp(chFunction,"Line ") == 0 )
	{
		/* save results in array so that they can be printed when done */
		wavelength[ipLine] = wl;
		strcpy( chLabSave[ipLine], chLab );
		xIntenSave[ipLine] = xInten;

		/* now increment the counter and then check if we have filled the line, 
		 * and so should print it */
		++ipLine;
		/* do print now if we are in column mode (one line per line) or if we have filled up
		 * the line */
		if( ( strcmp(punch.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
		{
			/* "array " is usual array 6 wide, "column" is one line per line */
			for( i=0; i < ipLine; i++ )
			{
				fprintf( io, " %4.4s ", chLabSave[i] );
				prt_wl( io, wavelength[i] );
				fprintf( io,"\t%.3e", xIntenSave[i]);
				/* >>chng 02 apr 24, do not print type */
				/* single column for input into data base */
				if( strcmp(punch.chPunRltType,"column") == 0 )				
					fprintf( io, "\n" );
			}
			/* only put cr if we did not just put one */
			if( strcmp(punch.chPunRltType,"array ") == 0 )				
				fprintf( io, " \n" );
			ipLine = 0;
		}
	}
	else if( strcmp(chFunction,"Flush") == 0 )
	{
		if( ipLine > 0 )
		{
			/* this is an option to print many emission lines across an output line,
			 * the array option, or a single column of numbers, the column option
			 * that appears on the "punch results" and "punch intensity" commands
			 */
			/* usual array 6 wide */
			for( i=0; i < ipLine; i++ )
			{
				fprintf( io, " %4.4s", chLabSave[i] );
				prt_wl( io, wavelength[i] );
				fprintf( io,"\t%.3e", xIntenSave[i]);
				/* >>chng 02 apr 24, do not print type */
				/* single column for input into data base */
				if( strcmp(punch.chPunRltType,"column") == 0 )				
					fprintf( io, "\n" );
			}
			if( strcmp(punch.chPunRltType,"array ") == 0 )
				fprintf( io, " \n" );
		}
	}
	else
	{
		fprintf( ioQQQ, " PunResults1Line called with insane request=%5.5s\n", 
		  chFunction );
		puts( "[Stop in result]" );
		cdEXIT(EXIT_FAILURE);
	}

#	ifdef DEBUG_FUN
	fputs( " <->PunResults1Line()\n", debug_fp );
#	endif
	return;
}
#define	NENR_GAUNT	37
#define	NTE_GAUNT	21
/*PunchGaunts called by punch gaunts command to output gaunt factors */
static void PunchGaunts(FILE* io)
{
	long int i, 
	  charge,
	  ite, 
	  j;

	float ener[NENR_GAUNT], 
	  ste[NTE_GAUNT],
	  z;
	double tempsave;
    double g[NENR_GAUNT][NTE_GAUNT], gfac;
	  

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

	/* this routine is called from the PUNCH GAUNTS command
	 * it drives the gaunt factor routine to punch gaunts over full range */
	tempsave = phycon.te;

	for( i=0; i < NTE_GAUNT; i++ )
	{
		ste[i] = 0.5f*i;
	}

	for( i=0; i < NENR_GAUNT; i++ )
	{
		ener[i] = 0.5f*i - 8.f;
	}

	for( charge = 1; charge<LIMELM; charge++ )
	{
		/* nuclear charge */
		z = (float)log10((double)charge);

		/* ener is log of energy */
		for( ite=0; ite < NTE_GAUNT; ite++ )
		{
			phycon.alogte = ste[ite];
			phycon.te = pow(10.,phycon.alogte);

			for( j=0; j < NENR_GAUNT; j++ )
			{
				gfac = CalcThermAveGaunt( phycon.te, (double)charge, pow(10.,(double)ener[j]) );
				/*fprintf(ioQQQ, "z %.2e ener %.2e temp %.2e gfac %.3e \n",
					pow(10.,z), pow(10.,ener[j]), (double)phycon.te, gfac );*/
				g[j][ite] = gfac;
			}
		}

		/* now punch out the results */
		fprintf( io, "      " );
		for( i=1; i <= NTE_GAUNT; i++ )
		{
			fprintf( io, "\t%6.1f", ste[i-1] );
		}
		fprintf( io, "\n" );

		for( j=0; j < NENR_GAUNT; j++ )
		{
			fprintf( io, "\t%6.2f", ener[j] );
			for( ite=0; ite < NTE_GAUNT; ite++ )
			{
				fprintf( io, "\t%6.2f", log10(g[j][ite]) );
			}
			fprintf( io, "\n" );
		}

		fprintf( io, "      " );
		for( i=0; i < NTE_GAUNT; i++ )
		{
			fprintf( io, "\t%6.1f", ste[i] );
		}
		fprintf( io, "\n" );

		/* now do the same thing as triplets */
		fprintf( io, "      " );
		for( i=0; i < NTE_GAUNT; i++ )
		{
			fprintf( io, "\t%6.1f", ste[i] );
		}
		fprintf( io, "\n" );

		for( i=0; i < NTE_GAUNT; i++ )
		{
			for( j=0; j < NENR_GAUNT; j++ )
			{
				fprintf( io, "\t%6.2f\t%6.2f\t%6.2e\n", ste[i], ener[j], 
				  log10(g[j][i]) );
			}
		}

		fprintf( io, "Below is log(u), log(gamma**2), gff\n" );
		/* do the same thing as above, but this time print log(u) and log(gamma2) instead of temp and energy.	*/
		for( i=0; i < NTE_GAUNT; i++ )
		{
			for( j=0; j < NENR_GAUNT; j++ )
			{
				fprintf( io, "\t%6.2f\t%6.2f\t%6.2e\n", 2.*z + log10(TE1RYD) - ste[i] , log10(TE1RYD)+ener[j]-ste[i], 
				log10(g[j][i]) );
			}
		}
		fprintf( io, "end of charge = %li\n", charge);
		fprintf( io, "****************************\n");
	}

	phycon.te = tempsave;
    phycon.alogte = log10( phycon.te );

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

#undef	NENR_GAUNT
#undef	NTE_GAUNT


