/home66/gary/public_html/cloudy/c08_branch/source/punch_do.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  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 emission, for AGN3 chapter 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 "secondaries.h"
#include "grainvar.h"
#include "lines.h"
#include "dynamics.h"
#include "colden.h"
#include "continuum.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 "hcmap.h"
#include "input.h"
#include "punch.h"
#include "optimize.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 * ioPUN, 
  /* 4 char + null string */
  const char *chLab, 
  realnum 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* ioPUN);

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

STATIC void PunchLineStuff(
  FILE * ioPUN,
  const char *chJob , 
  realnum xLimit);

/* punch h emission, for chapter 4, routine is below */
STATIC void AGN_Hemis(FILE *ioPUN );

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

/*PunLineIntensity produce the 'punch lines intensity' output */
STATIC void PunLineIntensity(FILE * ioPUN);

char *chDummy;

void PunchDo(
        /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
        const char *chTime) 
{
        bool 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, 
          stage, 
          sum, 
          texc, 
          xLinMax;

        DEBUG_ENTRY( "PunchDo()" );

        /* 
         * 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 iteration is complete 
         *
         * if "LAST" then this is last call to routine after iteration complete
         * punch only if "LAST" when results at end of iteration 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
         */

        /* return if no punch to do */
        if( punch.npunch < 1 )
        { 
                return;
        }

        /* during a grid calculation this routine saves grid points after
         * cloudy is called.  we may output it below */
        if( grid.lgGrid )
        {
                if( chTime[0]=='L' )
                        GridGatherInCloudy();

                /* save grid information */
                GridGatherAfterCloudy( chTime );
        }

        for( ipPun=0; ipPun < punch.npunch; ipPun++ )
        {
                /* this global variable to remember where in the punch stack we are */
                punch.ipConPun = ipPun;

                /* iterations.lgLastIt is true if this is last iteration
                 * lgPunLstIter set true if 'last' key occurred 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\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* Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon /
                                          SDIV( HFLines[0].Emis->TauCon ),
                                          /*TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ),*/
                                          HFLines[0].Lo->Pop ,
                                          HFLines[0].Hi->Pop ,
                                          OccupationNumberLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ),
                                          HFLines[0].Emis->TauCon , 
                                          Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon,
                                          HFLines[0].Emis->PopOpc,
                                          /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */
                                          (dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop)/SDIV( hyperfine.Tspin21cm),
                                          /* why was above multiplied by this following term? */
                                          /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/
                                          HFLines[0].Emis->TauIn,
                                          TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) ,
                                          colden.H0_ov_Tspin,
                                          /*>>chng 27 mar, GS, integrated 21cm spin temperature*/
                                          colden.H0_21cm_lower,
                                          colden.H0_21cm_upper,
                                          -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower)
                                          );
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"AGES") == 0 )
                        {
                                /* punch timescales vs depth */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
                                          /* depth, cm */
                                          radius.depth_mid_zone,
                                          /* cooling timescale */
                                          dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot, 
                                          /* H2 destruction timescale */
                                          timesc.time_H2_Dest_here, 
                                          /* CO destruction timescale */
                                          timesc.AgeCOMoleDest[findspecies("CO")->index], 
                                          /* OH destruction timescale */
                                          timesc.AgeCOMoleDest[findspecies("OH")->index], 
                                          /* 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 chapter 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();
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"ASSE") == 0 )
                        {
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        /* punch the assert output */
                                        lgCheckAsserts(punch.ipPnunit[ipPun]);
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"AVER") == 0 )
                        {
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        /* punch the averages output */
                                        punch_average( ipPun , "PUNCH", chDummy );
                                }
                        }

                        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 dynamics xxx, information about dynamical 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],"COLS") == 0 )
                        {
                                /* punch some column densities */
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        char chHeader[2];
                                        punch_colden(chHeader , punch.ipPnunit[ipPun] , "PUNS" );
                                }
                        }

                        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 is the usual "punch continuum" case */
                                /* >>chng 06 apr 03, add every option to do every zone */
                                /* if punch every is set then nPunchEveryZone must be positive 
                                 * was init to -1 */
                                bool lgPrintThis =false;
                                if( punch.lgPunchEveryZone[ipPun] )
                                {
                                        /* this branch, every option is on line so want to print every n zone */
                                        if( strcmp(chTime,"LAST") != 0 )
                                        {
                                                /* not last zone - print first and intermediate cases */
                                                if( nzone==1 )
                                                {
                                                        lgPrintThis = true;
                                                }
                                                else if( nzone%punch.nPunchEveryZone[ipPun]==0 )
                                                {
                                                        lgPrintThis = true;
                                                }
                                        }
                                        else
                                        {
                                                /* this is last zone, print only if did not trip on above */
                                                if( nzone!=1 && nzone%punch.nPunchEveryZone[ipPun]!=0 )
                                                {
                                                        lgPrintThis = true;
                                                }
                                        }
                                }
                                else
                                {
                                        /* this branch, not "every", so only print the last zone */
                                        if( strcmp(chTime,"LAST") == 0 )
                                                lgPrintThis = true;
                                }
                                ASSERT( !punch.lgPunchEveryZone[ipPun] || punch.nPunchEveryZone[ipPun]>0 );
                                if( lgPrintThis )
                                {
                                        if( punch.lgPunchEveryZone[ipPun] && nzone!=1)
                                                fprintf( punch.ipPnunit[ipPun], "%s\n",
                                                        punch.chHashString );

                                        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.flux_total_incident[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 * rfield.trans_coef_total[j];

                                                /* 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 ); 
                                                /* reflected 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[ipPun][0]);
                        }

                        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.flux_total_incident[j]*rfield.anu2[j]*
                                                  EN1RYD/rfield.widflx[j];
                                                /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh 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 in optically thin limit */
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        /* GrainMakeDiffuse broke out emission into types 
                                         * according to matType */
                                        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;

                                                /* anu is .5e format to resolve energy mesh near 1 Ryd 
                                                 * AnuUnit gives anu in whatever units were set with units option */
                                                fprintf( punch.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n", 
                                                  AnuUnit(rfield.AnuOrg[j]) , fsum , fout , 
                                                  /* total emission per unit volume defined in GrainMakeDiffuse
                                                   * used in RT_diffuse to form total grain emission */
                                                  gv.GrainEmission[j]*rfield.anu2[j]*
                                                  EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo );
                                        }
                                }
                        }

                        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" );
                                                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.flux_total_incident[j] > 1e-25 )
                                                {
                                                        av = rfield.ConRefIncid[j]/rfield.flux_total_incident[j];
                                                }
                                                else
                                                {
                                                        av = 0.;
                                                }
                                                /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
                                                fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.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 .5e to resolve energy mesh near 1 Ryd */
                                                fprintf( punch.ipPnunit[ipPun], "%.5e\t%.5e\n", 
                                                  AnuUnit(rfield.AnuOrg[j]), 
                                                  rfield.ConEmitLocal[nzone][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 .5e to resolve energy mesh 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[ipPun][0] > 0. )
                                                /* possible lower bounds to energy range - 
                                                 * 0 if not set with range option*/
                                                j = ipFineCont( punch.punarg[ipPun][0] );
                                        else
                                                j = 0;

                                        /* upper limit set with range option */
                                        if( punch.punarg[ipPun][1]> 0. )
                                                nu_hi = ipFineCont( punch.punarg[ipPun][1]);
                                        else
                                                nu_hi = rfield.nfine;

                                        /* number of cells to bring together, default is 10 */
                                        nskip = (long)punch.punarg[ipPun][2];

                                        do
                                        {
                                                realnum sum1 = rfield.fine_opt_depth[j];
                                                realnum xnu = rfield.fine_anu[j];
                                                for( jj=1; jj<nskip; ++jj )
                                                {
                                                        xnu += rfield.fine_anu[j+jj];
                                                        sum1 += rfield.fine_opt_depth[j+jj];
                                                }
                                                fprintf( punch.ipPnunit[ipPun], 
                                                        "%.6e\t%.3e\n", 
                                                        AnuUnit(xnu/nskip), 
                                                        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[ipPun][0] <= 0. )
                                        {
                                                i1 = 1;
                                        }
                                        else if( punch.punarg[ipPun][0] < 100. )
                                        {
                                                i1 = ipoint(punch.punarg[ipPun][0]);
                                        }
                                        else
                                        {
                                                i1 = (long int)punch.punarg[ipPun][0];
                                        }

                                        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[ipPun][2]>0) || (strcmp(chTime,"LAST") == 0) )
                                {
                                        /* this flag will remember whether we have ever printed anything */
                                        bool lgPrt=false;
                                        if( punch.punarg[ipPun][2]>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[ipPun][0] <= 0. )
                                        {
                                                i1 = 1;
                                        }
                                        else if( punch.punarg[ipPun][0] < 100. )
                                        {
                                                i1 = ipoint(punch.punarg[ipPun][0]);
                                        }
                                        else
                                        {
                                                i1 = (long int)punch.punarg[ipPun][0];
                                        }

                                        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 percent) */
                                                /* >>chng 01 nov 22, format to c-friendly */
                                                if( (RateInter >= punch.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) )
                                                {
                                                        lgPrt = true;
                                                        /* >>chng 96 oct 22, format of anu to 11.5 to resolve energy mesh 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 and fsum were %.2e %.2e\n",
                                                        sum,fsum);
                                                fprintf(punch.ipPnunit[ipPun],
                                                        " punchdo, the low-frequency energy was %.5e Ryd\n",
                                                        rfield.anu[i1-1]);
                                                fprintf(punch.ipPnunit[ipPun],
                                                        " You can reset the threshold for the lowest fractional "
                                                        "interaction to print with the second number of the punch command\n"
                                                        " The fraction was %.3f and this was too large.\n",
                                                        punch.punarg[ipPun][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 - this is not the main "punch continuum"
                                 * command - search on "CON " above 
                                 * 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 
                                                 * >>chng 07 apr 26 include transmission coefficient */
                                                flxatt = rfield.flux[j]*rfield.anu2[j]*EN1RYD/
                                                  rfield.widflx[j]*radius.r1r0sq * rfield.trans_coef_total[j];

                                                /*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 multiplied by an2,
                                                 * quantity should be photons per Ryd, 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],"%.3e\t", AnuUnit(rfield.AnuOrg[j]) );
                                                fprintf( punch.ipPnunit[ipPun],"%.3e\t", flxtrn);
                                                fprintf( punch.ipPnunit[ipPun],"%.3e\n", rfield.trans_coef_total[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]) );
                                        }
                                }
                        }

                        /* punch grain abundance and punch grain D/G ratio commands */
                        else if( strcmp(punch.chPunch[ipPun],"DUSA") == 0 ||
                                 strcmp(punch.chPunch[ipPun],"DUSD") == 0 )
                        {
                                bool lgDGRatio = ( strcmp(punch.chPunch[ipPun],"DUSD") == 0 );

                                /* grain abundance */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        /* used to print header exactly one time */
                                        static bool lgMustPrtHeaderDRRatio = true,
                                                lgMustPrtHeaderAbundance=true;
                                        /* print grain headder first if this has not yet been done */
                                        if( ( lgMustPrtHeaderDRRatio && lgDGRatio ) || 
                                            ( lgMustPrtHeaderAbundance && !lgDGRatio ) )
                                        {
                                                /* only print one header for each case, but must
                                                 * track separately if both used in same sim */
                                                if( lgMustPrtHeaderDRRatio && lgDGRatio )
                                                        lgMustPrtHeaderDRRatio = false;
                                                else if( lgMustPrtHeaderAbundance &&!lgDGRatio )
                                                        lgMustPrtHeaderAbundance = false;

                                                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], "\ttotal\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], "\txxxx\n" );
                                        }
                                        fprintf( punch.ipPnunit[ipPun], " %.5e", 
                                                radius.depth_mid_zone );
                                        /* grain abundance per bin in g/cm^3 */
                                        double total = 0.;
                                        for( nd=0; nd < gv.nBin; ++nd ) 
                                        {
                                                double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]*
                                                        gv.bin[nd]->cnv_H_pCM3;
                                                if( lgDGRatio )
                                                        abund /= dense.xMassDensity;
                                                fprintf( punch.ipPnunit[ipPun], "\t%.3e", abund );
                                                total += abund;
                                        }
                                        fprintf( punch.ipPnunit[ipPun], "\t%.3e\n", total );
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"DUSP") == 0 )
                        {
                                /* grain potential */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        /* used to print header exactly one time */
                                        static bool lgMustPrtHeader = true;
                                        /* 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 bool lgMustPrtHeader = true;
                                        /* 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 bool lgMustPrtHeader = true;
                                        /* 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 grain charge - eden from grains and 
                                 * charge per grain in electrons / grain */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        /* used to print header exactly one time */
                                        static bool lgMustPrtHeader = true;
                                        /* 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 ,
                                                /* electron density contributed by grains, in e/cm^3, 
                                                 * positive number means grain supplied free electrons */
                                                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 bool 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 bool 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 )
                                {
                                        realnum 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[ipPun][0];
                                        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[ipPun][1] == 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", 
                                                                findspecies("CO")->hevmol/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 AGN3 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 continuum, check that FeII is turned on */
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        for( j=0; j < nFeIIConBins; j++ )
                                        {
                                                /* emission from large FeII atom, integrated over continuum intervals 
                                                 * units are always in units, erg cm-2 s-1 as in intensity case,
                                                 * even if luminosity case is considered 
                                                 * [j][0] is flux, [j][1] and [j][2] are upper and
                                                 * lower bounds in Angstroms.  
                                                 * these are set in FeIIZero */
                                                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, routine is in atom_feii.c */
                                        FeIIPunchLines( punch.ipPnunit[ipPun] );
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"FE2o") == 0 )
                        {
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        /* punch line optical depths, routine is in atom_feii.c */
                                        FeIIPunchOpticalDepth( punch.ipPnunit[ipPun] );
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"FRED") == 0 )
                        {
                                /* set with punch Fred command, this punches some stuff from
                                 * Fred Hamann's dynamics project */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        /* Fred's list */
                                        fprintf( punch.ipPnunit[ipPun], "%.5e\t%.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\n", 
                                                radius.Radius, radius.depth ,wind.windv/1e5,
                                                dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te,
                                                wind.AccelLine , wind.AccelCont ,
                                                wind.fmul , 
                                                mean.xIonMeans[0][ipHYDROGEN][0] , mean.xIonMeans[0][ipHYDROGEN][1] ,
                                                mean.xIonMeans[0][ipHELIUM][0] , mean.xIonMeans[0][ipHELIUM][1] ,
                                                mean.xIonMeans[0][ipHELIUM][2] ,
                                                mean.xIonMeans[0][ipCARBON][1] , mean.xIonMeans[0][ipCARBON][2] ,
                                                mean.xIonMeans[0][ipCARBON][3] ,
                                                mean.xIonMeans[0][ipOXYGEN][0] , mean.xIonMeans[0][ipOXYGEN][1] ,
                                                mean.xIonMeans[0][ipOXYGEN][2] , mean.xIonMeans[0][ipOXYGEN][3] ,
                                                mean.xIonMeans[0][ipOXYGEN][4] , mean.xIonMeans[0][ipOXYGEN][5] ,
                                                mean.xIonMeans[0][ipOXYGEN][6] , mean.xIonMeans[0][ipOXYGEN][7] ,
                                                dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] ,
                                                dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] ,
                                                dense.xIonDense[ipHELIUM][2] ,
                                                dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] ,
                                                dense.xIonDense[ipCARBON][3] ,
                                                dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] ,
                                                dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] ,
                                                dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] ,
                                                dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] ,
                                                mean.xIonMeans[0][ipMAGNESIUM][1] , dense.xIonDense[ipMAGNESIUM][1]);
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"FE2d") == 0 )
                        {
                                /* punch some departure coefficients 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 coefficients for large FeII atom */
                                if( strcmp(chTime,"LAST") != 0 )
                                        FeIIPunDepart(punch.ipPnunit[ipPun],true);
                        }

                        else if( strcmp(punch.chPunch[ipPun],"FE2p") == 0 )
                        {
                                bool lgFlag = false;
                                if( punch.punarg[ipPun][2] )
                                        lgFlag = true;
                                /* punch small subset of level populations for large FeII atom */
                                if( strcmp(chTime,"LAST") != 0 )
                                        FeIIPunPop(punch.ipPnunit[ipPun],false,0,0,
                                        lgFlag);
                        }

                        else if( strcmp(punch.chPunch[ipPun],"FE2P") == 0 )
                        {
                                bool lgFlag = false;
                                if( punch.punarg[ipPun][2] )
                                        lgFlag = true;
                                /* punch range of level populations for large FeII atom */
                                if( strcmp(chTime,"LAST") != 0 )
                                        FeIIPunPop(punch.ipPnunit[ipPun],
                                        true,
                                        (long int)punch.punarg[ipPun][0],
                                        (long int)punch.punarg[ipPun][1],
                                        lgFlag);
                        }

                        /* punch spectra in fits format */
                        else if( strcmp(punch.chPunch[ipPun],"FITS") == 0 )
                        {
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        punchFITSfile( punch.ipPnunit[ipPun], NUM_OUTPUT_TYPES );
                                }
                        }
                        /* 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 < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ )
                                                                {
                                                                        fprintf( punch.ipPnunit[ipPun], "%5.2f",
                                                                                 t_yield::Inst().elec_eject_frac(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[ipPun][0];
                                        ion = (long)punch.punarg[ipPun][1];
                                        /* 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]);
                        }

                        /* punch parameters of grid calculation */
                        else if( strcmp(punch.chPunch[ipPun],"GRID") == 0 )
                        {
                                /* one time punch of all grid parameters after grid is complete */
                                if( (strcmp(chTime,"LAST") == 0 ) && grid.lgGridDone )
                                {
                                        /* start of line gives abort and warning summary */     
                                        fprintf(punch.ipPnunit[ipPun], "#Abort?\tWarnings?" );
                                        /* print start of each variable command line */
                                        for( i=0; i< grid.nintparm; i++ )
                                        {
                                                char chStr[10];
                                                strncpy( chStr , optimize.chVarFmt[i] , 9 );
                                                /* make sure this small bit of string is terminated */
                                                chStr[9] = '\0';
                                                fprintf(punch.ipPnunit[ipPun], "\t%s", chStr );
                                        }
                                        fprintf(punch.ipPnunit[ipPun], "\n" );
                                        for( i=0; i<grid.totNumModels; i++ )
                                        {
                                                /* abort / warning summary for this sim */
                                                fprintf(punch.ipPnunit[ipPun], "%c\t%c",
                                                        TorF(grid.lgAbort[i]) , 
                                                        TorF(grid.lgWarn[i]) );
                                                /* the grid parameters */
                                                for( j=0; j< grid.nintparm; j++ )
                                                {
                                                        fprintf(punch.ipPnunit[ipPun], "\t%f", grid.interpParameters[i][j] );
                                                }
                                                fprintf(punch.ipPnunit[ipPun], "\n" );
                                        }
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"HISp") == 0 )
                        {
                                /* punch pressure history of current zone */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        /* note if pressure convergence failure occurred 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 failure occurred 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 failure occurred 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 failure occurred 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]);
                        }

                        else if( strncmp(punch.chPunch[ipPun],"HE",2) == 0 )
                        {
                                /* various punch helium commands */
                                /* punch helium line wavelengths */
                                if( strcmp(punch.chPunch[ipPun] , "HELW") == 0 )
                                {
                                        if( strcmp(chTime,"LAST") == 0 )
                                        {
                                                /* punch helium & he-like wavelengths, first header */
                                                fprintf( punch.ipPnunit[ipPun], 
                                                        "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S"
                                                        "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" );
                                                fprintf( punch.ipPnunit[ipPun], "\n" );
                                                for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
                                                {
                                                        /* print element name, nuclear charge */
                                                        fprintf( punch.ipPnunit[ipPun], "%li\t%s", 
                                                                nelem+1 , elementnames.chElementSym[nelem] );
                                                        /*prt_wl print floating wavelength in Angstroms, in output format */
                                                        fprintf( punch.ipPnunit[ipPun], "\t" );
                                                        prt_wl( punch.ipPnunit[ipPun] , 
                                                                Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].WLAng );
                                                        fprintf( punch.ipPnunit[ipPun], "\t" );
                                                        prt_wl( punch.ipPnunit[ipPun] , 
                                                                Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe1s1S].WLAng );
                                                        fprintf( punch.ipPnunit[ipPun], "\t" );
                                                        prt_wl( punch.ipPnunit[ipPun] , 
                                                                Transitions[ipHE_LIKE][nelem][ipHe2p3P2][ipHe1s1S].WLAng );
                                                        fprintf( punch.ipPnunit[ipPun], "\t" );
                                                        prt_wl( punch.ipPnunit[ipPun] , 
                                                                Transitions[ipHE_LIKE][nelem][ipHe2s3S][ipHe1s1S].WLAng );
                                                        fprintf( punch.ipPnunit[ipPun], "\t" );
                                                        prt_wl( punch.ipPnunit[ipPun] , 
                                                                Transitions[ipHE_LIKE][nelem][ipHe2p3P2][ipHe2s3S].WLAng );
                                                        fprintf( punch.ipPnunit[ipPun], "\t" );
                                                        prt_wl( punch.ipPnunit[ipPun] , 
                                                                Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe2s3S].WLAng );
                                                        fprintf( punch.ipPnunit[ipPun], "\t" );
                                                        prt_wl( punch.ipPnunit[ipPun] , 
                                                                Transitions[ipHE_LIKE][nelem][ipHe2p3P0][ipHe2s3S].WLAng );
                                                        fprintf( punch.ipPnunit[ipPun], "\n"); 
                                                }
                                        }
                                }
                                else
                                        TotalInsanity();
                        }

                        /* punch hummer, results needed for Lya transport, to feed into David's routine */
                        else if( strcmp(punch.chPunch[ipPun],"HUMM") == 0 )
                        {
                                eps = Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul/
                                  (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.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, 
                                  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn, 
                                  StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1], 
                                  StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*dense.xIonDense[ipHYDROGEN][1], 
                                  phycon.te, Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->damp, eps );
                        }

                        else if( strncmp( punch.chPunch[ipPun] , "HYD", 3 ) == 0 )
                        {
                                /* various punch hydrogen commands */
                                if( strcmp(punch.chPunch[ipPun],"HYDc") == 0 )
                                {
                                        if( strcmp(chTime,"LAST") != 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*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop;
                                                fref = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop;
                                                fout = StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop;
                                                /* 06 aug 28, from numLevels_max to _local. */
                                                for( ion=ipH2s; ion < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ion++ )
                                                {
                                                        /* this is total decays to ground */
                                                        RateInter += 
                                                                Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Aul*
                                                          (Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pesc + 
                                                          Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pelec_esc + 
                                                          Transitions[ipH_LIKE][ipHYDROGEN][ion][ipH1s].Emis->Pdest);
                                                        /* total photo from all levels */
                                                        fref += iso.gamnc[ipH_LIKE][ipHYDROGEN][ion]*StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop;
                                                        /* total col ion from all levels */
                                                        stage += iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ion]*dense.eden*
                                                          StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop;
                                                        fout += StatesElem[ipH_LIKE][ipHYDROGEN][ion].Pop;
                                                }
                                                fprintf( punch.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e", 
                                                  nzone, 
                                                  /* photo and collision 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+ */
                                                  safe_div( 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 )
                                {
                                        if( strcmp(chTime,"LAST") != 0 )
                                        {
                                                /* punch hydrogen populations 
                                                 * first give total atom and ion density [cm-3]*/
                                                fprintf( punch.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e", 
                                                  radius.depth_mid_zone, 
                                                  dense.xIonDense[ipHYDROGEN][0], 
                                                  dense.xIonDense[ipHYDROGEN][1] );

                                                /* next give state-specific densities [cm-3] */
                                                for( j=ipH1s; j < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]-1; j++ )
                                                {
                                                        fprintf( punch.ipPnunit[ipPun], "\t%.2e", 
                                                                StatesElem[ipH_LIKE][ipHYDROGEN][j].Pop*dense.xIonDense[ipHYDROGEN][1] );
                                                }
                                                fprintf( punch.ipPnunit[ipPun], "\n" );
                                        }
                                }

                                else if( strcmp(punch.chPunch[ipPun],"HYDl") == 0 )
                                {
                                        if( strcmp(chTime,"LAST") == 0 )
                                        {
                                                /* punch hydrogen line 
                                                 * gives intensities and optical depths */
                                                for( ipHi=4; ipHi<iso.numLevels_local[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./Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].EnergyWN,
                                                                        Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Emis->TauIn );
                                                        }
                                                }
                                        }
                                }

                                /* punch hydrogen Lya - some details about Lya */
                                else if( strcmp(punch.chPunch[ipPun],"HYDL") == 0 )
                                {
                                        if( strcmp(chTime,"LAST") != 0 )
                                        {
                                                /* the population ratio for Lya */
                                                double popul = StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop/SDIV(StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop);
                                                /* the excitation temperature of Lya */
                                                texc = TexcLine( &Transitions[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,
                                                  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn, 
                                                  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauTot, 
                                                  popul, 
                                                  texc, 
                                                  phycon.te, 
                                                  texc/phycon.te ,
                                                  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pesc, 
                                                  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest, 
                                                  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->pump, 
                                                  opac.opacity_abs[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1],
                                                  opac.albedo[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] );
                                        }
                                }

                                else if( strcmp(punch.chPunch[ipPun],"HYDr") == 0 )
                                {
                                        /* punch hydrogen recc - recombination cooling for AGN3 */
                                        TempChange(2500.f, false);
                                        while( phycon.te <= 20000. )
                                        {
                                                double r1;
                                                double ThinCoolingCaseB; 

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

                                                TempChange(phycon.te *2.f , false);
                                        }
                                        /* must exit since we have disturbed the solution */
                                        fprintf(ioQQQ , "punch agn now exits since solution is disturbed.\n");
                                        cdEXIT( EXIT_SUCCESS );
                                }
                                else
                                        TotalInsanity();
                        }

                        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[ipPun][0];
                                        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 over read -
                                         * 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 */
                                                const int 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 = t_ADfA::Inst().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,chDummy); 
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"LIND") == 0 )
                        {
                                /* punch line data, then stop */
                                PunchLineData(punch.ipPnunit[ipPun]);
                        }

                        else if( strcmp(punch.chPunch[ipPun],"LINL") == 0 )
                        {
                                /* punch line labels */
                                bool lgPrintAll=false;
                                /* LONG keyword on punch line labels command sets this to 1 */
                                if( punch.punarg[ipPun][0]>0. )
                                        lgPrintAll = true;
                                prt_LineLabels(punch.ipPnunit[ipPun] , lgPrintAll );
                        }

                        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[ipPun][0]);
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"LINP") == 0 )
                        {
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        static bool 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[ipPun][0]);
                                        if( lgFirst )
                                        {
                                                lgFirst = false;
                                                PunchLineStuff(punch.ipPnunit[ipPun],"populat" , punch.punarg[ipPun][0]);
                                        }
                                }
                        }

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

                        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 )
                        {
                                /* punch line array */
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        /* punch out all lines with energies */
                                        for( j=0; j < LineSave.nsum; j++ )
                                        {
                                                if( LineSv[j].wavelength > 0. && 
                                                        LineSv[j].sumlin[LineSave.lgLineEmergent] > 0. )
                                                {
                                                        /* line energy, in units set with units option */
                                                        fprintf( punch.ipPnunit[ipPun], "%12.5e", 
                                                          AnuUnit((realnum)RYDLAM/LineSv[j].wavelength) );
                                                        /* line label */
                                                        fprintf( punch.ipPnunit[ipPun], "\t%4.4s ",
                                                                LineSv[j].chALab );
                                                        /* wavelength */
                                                        prt_wl( punch.ipPnunit[ipPun], LineSv[j].wavelength );
                                                        /* intrinsic intensity */
                                                        fprintf( punch.ipPnunit[ipPun], "\t%8.3f", 
                                                                log10(SDIV(LineSv[j].sumlin[0]) ) + radius.Conv2PrtInten );
                                                        /* emergent line intensity, r recombination  */
                                                        fprintf( punch.ipPnunit[ipPun], "\t%8.3f", 
                                                                log10(SDIV(LineSv[j].sumlin[1]) ) + radius.Conv2PrtInten );
                                                        /* type of line, i for info, etc */
                                                        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 */
                                        const int NLINE_H2 = 31; 
                                        /* the number of lines which are not H2 */
                                        const int 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 Angstroms 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 positive return says didn't find it */
                                                /* arguments are 4-char label, wavelength, return log total intensity, linear rel inten */
                                                if( cdLine( chLabel[n] , (realnum)(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 positive return says didn't find it */
                                                        if( cdLine( "H2  " , (realnum)(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);
                                                        }
                                                }
                                        }
                                }
                        }

                        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],
                                        findspecies("CO")->hevmol,
                                        findspecies("O2")->hevmol,
                                        findspecies("CH")->hevmol,
                                        findspecies("OH")->hevmol,
                                        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,
                                        findspecies("CO")->hevcol,
                                        findspecies("O2")->hevcol,
                                        findspecies("CH")->hevcol,
                                        findspecies("OH")->hevcol,
                                        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_findrk("PHOTON,CO=>C,O"),
                                        /* 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].Coll.cool,
                                        TauLines[ipT370].Coll.cool,
                                        TauLines[ipT157].Coll.cool,
                                        TauLines[ipT63].Coll.cool,
                                        TauLines[ipT146].Coll.cool );
                                }
                        }

                        else if( strcmp( punch.chPunch[ipPun],"LLST") == 0)
                        {
                                /* punch linelist command - do on last iteration */
                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        fprintf( punch.ipPnunit[ipPun], 
                                                "%li" , iteration );

                                        /* -1 is flag saying that this punch command was not set */
                                        if( punch.nLineList[ipPun] < 0 )
                                                TotalInsanity();

                                        /* loop over all lines in the file we read */
                                        for( j=0; j<punch.nLineList[ipPun]; ++j )
                                        {
                                                double relative , absolute, PrtQuantity;
                                                if( (cdLine( punch.chLineListLabel[ipPun][j] , 
                                                        punch.wlLineList[ipPun][j]  , 
                                                        &relative , &absolute ))<=0 )
                                                {
                                                        if( !h2.lgH2ON && strncmp( punch.chLineListLabel[ipPun][j] , "H2  " , 4 )==0 )
                                                        {
                                                                static bool lgMustPrintFirstTime = true;
                                                                if( lgMustPrintFirstTime )
                                                                {
                                                                        /* it's an H2 line and H2 is not being done - ignore it */
                                                                        fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
                                                                                "included, so I will ignore it.  Log intensity set to -30.\n" );
                                                                        fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
                                                                        lgMustPrintFirstTime = false;
                                                                }
                                                                relative = -30.f;
                                                                absolute = -30.f;
                                                        }
                                                        else if( lgAbort )
                                                        {
                                                                /* we are in abort mode */
                                                                relative = -30.f;
                                                                absolute = -30.f;
                                                        }
                                                        else
                                                        {
                                                                fprintf(ioQQQ,"DISASTER - did not find a line in the Line List table\n");
                                                                cdEXIT(EXIT_FAILURE);
                                                        }
                                                }

                                                /* options to do either relative or absolute intensity
                                                 * default is relative, is absolute keyword on line then
                                                 * punarg set to 1 */
                                                /* straight line intensities */
                                                if( punch.punarg[ipPun][0] > 0 )
                                                        PrtQuantity = absolute;
                                                else
                                                        PrtQuantity = relative;
                                                /* if taking ratio print every other line as ratio
                                                 * with previous line */
                                                if( punch.lgLineListRatio[ipPun] )
                                                {
                                                        /* do line pair ratios */
                                                        static double SaveQuantity = 0;
                                                        if( is_odd(j) )
                                                                fprintf( punch.ipPnunit[ipPun], "\t%.3e" , 
                                                                        SaveQuantity / SDIV( PrtQuantity ) );
                                                        else
                                                                SaveQuantity = PrtQuantity;
                                                }
                                                else
                                                {
                                                        fprintf( punch.ipPnunit[ipPun], "\t%.3e" , PrtQuantity );
                                                }
                                        }
                                        fprintf( punch.ipPnunit[ipPun], "\n" );
                                }
                        }

                        else if( strcmp( punch.chPunch[ipPun],"RCO ") == 0)
                        {
                                /* punch chemistry CO rates command */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        CO_punch_mol(punch.ipPnunit[ipPun],"CO",NULL,radius.depth_mid_zone);
                                }
                        }

                        /*>>chng 06 jul 13, Nick Abel add OH rates */
                        else if( strcmp( punch.chPunch[ipPun],"ROH ") == 0)
                        {
                                /* punch chemistry OH rates command */
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        CO_punch_mol(punch.ipPnunit[ipPun],"OH",NULL,radius.depth_mid_zone);
                                }
                        }
                        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(  !hcmap.lgMapDone &&
                                        (nzone == hcmap.MapZone  ||  strcmp(chTime,"LAST") == 0 ) )
                                {
                                        lgTlkSav = called.lgTalk;
                                        called.lgTalk = true;
                                        hcmap.lgMapBeingDone = true;
                                        map_do(punch.ipPnunit[ipPun] , " map");
                                        called.lgTalk = lgTlkSav;
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"MOLE") == 0 )
                        {
                                if( strcmp(chTime,"LAST") != 0 )
                                {
                                        static bool lgMustPrintHeader = true;
                                        if( lgMustPrintHeader )
                                        {
                                                lgMustPrintHeader = false;
                                                fprintf( punch.ipPnunit[ipPun], 
                                                        "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate");

                                                for(i=0; i<N_H_MOLEC; ++i )
                                                {
                                                        fprintf( punch.ipPnunit[ipPun], "\t%s", hmi.chLab[i] );
                                                }
                                                for(i=0; i<mole.num_comole_calc; ++i )
                                                {
                                                        fprintf( punch.ipPnunit[ipPun], "\t%s",COmole[i]->label );
                                                }
                                                fprintf( punch.ipPnunit[ipPun], "\n");
                                        }
                                        /* 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_findrk("PHOTON,CO=>C,O") );

                                        /* 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<mole.num_comole_calc; ++j )
                                        {
                                                fprintf(punch.ipPnunit[ipPun],"\t%.2e",COmole[j]->hevmol );
                                        }
                                        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[ipPun][0] > 0. )
                                                /* possible lower bounds to energy range - will be zero if not set */
                                                j = ipFineCont( punch.punarg[ipPun][0] );
                                        else
                                                j = 0;

                                        /* upper limit */
                                        if( punch.punarg[ipPun][1]> 0. )
                                                nu_hi = ipFineCont( punch.punarg[ipPun][1]);
                                        else
                                                nu_hi = rfield.nfine;

                                        /* we will bring nskip cells together into one printed
                                         * number to make output smaller - default is 10 */
                                        nskip = (long)punch.punarg[ipPun][2];
                                        do
                                        {
                                                realnum sum1 = rfield.fine_opt_depth[j];
                                                realnum sum2 = rfield.fine_opac_zone[j];
                                                /* want to report the central wavelength of the cell */
                                                realnum xnu = rfield.fine_anu[j];
                                                for( jj=1; jj<nskip; ++jj )
                                                {
                                                        sum1 += rfield.fine_opt_depth[j+jj];
                                                        sum2 += rfield.fine_opac_zone[j+jj];
                                                        xnu += rfield.fine_anu[j+jj];
                                                }
                                                if( sum2 > 0. )
                                                        fprintf( punch.ipPnunit[ipPun], 
                                                          "%12.6e\t%.3e\t%.3e\n", 
                                                          AnuUnit(xnu/nskip), 
                                                          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 fractions 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 arg1 = SDIV(dense.gas_phase[ipHELIUM]);
                                                arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 );
                                                fprintf( punch.ipPnunit[ipPun], "\t%.4f", 
                                                  log10(arg1) );
                                        }

                                        /* carbon monoxide molecular fraction of CO */
                                        hold = SDIV(dense.gas_phase[ipCARBON]);
                                        hold = findspecies("CO")->hevmol/hold;
                                        hold = MAX2(toosmall, 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]), 
                                          findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 
                                          findspecies("H2O")->hevmol/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\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
                                          radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], 
                                          dense.eden, thermal.htot, wind.AccelTot, geometry.FillFac );
                                }
                        }

                        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", 
                                          /*A 1 #P depth */
                                          radius.depth_mid_zone, 
                                          /*B 2 Pcorrect */
                                          pressure.PresTotlCorrect, 
                                          /*C 3 Pcurrent */
                                          pressure.PresTotlCurr, 
                                          /*D 4 Pln + pintg 
                                           * >>chng 06 apr 19, subtract pinzon the acceleration added in this zone
                                           * since is not total at outer edge of zone, above is at inner edge */
                                          pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon, 
                                          /*E 5 pgas (0) */
                                          pressure.PresTotlInit, 
                                          /*F 6 Pgas */
                                          pressure.PresGasCurr, 
                                          /*G 7 Pram */
                                          pressure.PresRamCurr, 
                                          /*H 8 P rad in lines */
                                          pressure.pres_radiation_lines_curr, 
                                          /*I 9 Pinteg subtract continuum rad pres which has already been added on */
                                          pressure.PresInteg - pressure.pinzon, 
                                          /*J 10 V(wind km/s) wind speed in km/s */
                                          wind.windv/1e5,
                                          /*K cad(km/s) sound speed in km/s */
                                          timesc.sound_speed_adiabatic/1e5,
                                          /* the magnetic pressure */
                                          magnetic.pressure ,
                                          /* the local turbulent velocity in km/s */
                                          DoppVel.TurbVel/1e5 ,
                                          /* turbulent pressure */
                                          pressure.PresTurbCurr*DoppVel.lgTurb_pressure ,
                                          TorF(conv.lgConvPres) );
                                }
                        }

                        else if( punch.chPunch[ipPun][0]=='R' )
                        {
                                /* work around internal limits to Microsoft vs compiler */
                                if( strcmp(punch.chPunch[ipPun],"RADI") == 0 )
                                {
                                        /* punch radius information for all zones */
                                        if( strcmp(chTime,"LAST") != 0 )
                                        {
                                                fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n", 
                                                  nzone, radius.Radius_mid_zone, radius.depth_mid_zone, 
                                                  radius.drad );
                                        }
                                }

                                else if( strcmp(punch.chPunch[ipPun],"RADO") == 0 )
                                {
                                        /* punch radius information for only the last zone */
                                        if( strcmp(chTime,"LAST") == 0 )
                                        {
                                                fprintf( punch.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.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
                                {
                                        /* this can't happen */
                                        TotalInsanity();
                                }
                        }

                        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], 
                                                        "%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 
                                                  rfield.anu[j], 
                                                  rfield.ConEmitLocal[nzone][j]/rfield.widflx[j], 
                                                  opac.opacity_abs[j], 
                                                  rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j]), 
                                                  rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j])/plankf(j) );
                                        }
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"SOUD") == 0 )
                        {
                                /* parts of continuum source function vs depth
                                 * command was punch source function depth */
                                j = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] + 2;
                                fprintf( punch.ipPnunit[ipPun], 
                                        "%.4e\t%.4e\t%.4e\t%.4e\n", 
                                  opac.TauAbsFace[j-1], 
                                  rfield.ConEmitLocal[nzone][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 */
                                fprintf( punch.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e", 
                                        radius.depth_mid_zone, 
                                        phycon.te, 
                                        thermal.dCooldT );
                                /* if second zone then have one deriv */
                                if( nzone >1 )
                                {
                                        deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
                                        fprintf( punch.ipPnunit[ipPun], "\t%.2e", 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");
                        }

                        /* time dependent model */
                        else if( strcmp(punch.chPunch[ipPun],"TIMD") == 0 )
                        {
                                if( strcmp(chTime,"LAST") == 0 )
                                        DynaPunchTimeDep( punch.ipPnunit[ipPun] , "END" );
                        }

                        /* execution time per zone */
                        else if( strcmp(punch.chPunch[ipPun],"XTIM") == 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 acceleration, and ratio total
                                 * to electron scattering acceleration */
                                /* first test only punch last zone */
                                if( (punch.punarg[ipPun][0] == 0 && strcmp(chTime,"LAST") == 0)
                                        ||
                                        /* this test punch all zones */
                                        (punch.punarg[ipPun][0] == 1  && strcmp(chTime,"LAST") != 0 ) )
                                {
                                        fprintf( punch.ipPnunit[ipPun], 
                                                "%.5e\t%.5e\t%.4e\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 ,
                                                wind.AccelGravity );
                                }
                        }

                        else if( strcmp(punch.chPunch[ipPun],"XATT") == 0 )
                        {
                                /* attenuated incident continuum */
                                ASSERT( grid.lgOutputTypeOn[2] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 2 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XRFI") == 0 )
                        {
                                /* reflected incident continuum */
                                ASSERT( grid.lgOutputTypeOn[3] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 3 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XINC") == 0 )
                        {
                                /* incident continuum */
                                ASSERT( grid.lgOutputTypeOn[1] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 1 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XDFR") == 0 )
                        {
                                /* reflected diffuse continuous emission */
                                ASSERT( grid.lgOutputTypeOn[5] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 5 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XDFO") == 0 )
                        {
                                /* diffuse continuous emission outward */
                                ASSERT( grid.lgOutputTypeOn[4] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 4 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XLNR") == 0 )
                        {
                                /* reflected lines */
                                ASSERT( grid.lgOutputTypeOn[7] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 7 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XLNO") == 0 )
                        {
                                /* outward lines */
                                ASSERT( grid.lgOutputTypeOn[6] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 6 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XREF") == 0 )
                        {
                                /* total reflected, lines and continuum */
                                ASSERT( grid.lgOutputTypeOn[9] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 9 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XTRN") == 0 )
                        {
                                /* total outward, lines and continuum */
                                ASSERT( grid.lgOutputTypeOn[8] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 8 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        else if( strcmp(punch.chPunch[ipPun],"XSPM") == 0 )
                        {
                                /* exp(-tau) to the illuminated face */
                                ASSERT( grid.lgOutputTypeOn[10] );

                                if( strcmp(chTime,"LAST") == 0 )
                                {
                                        if( grid.lgGrid )
                                                punchFITSfile( punch.ipPnunit[ipPun], 10 );
                                        else
                                        {
                                                fprintf( ioQQQ," Cannot punch xspec files unless doing a grid.\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                        /* there are a few "punch" commands that are handled elsewhere
                         * punch dr is an example.  These will have lgReadPunch set
                         * false */
                        else if( punch.lgRealPunch[ipPun] )
                        {
                                /* this is insanity, internal flag set in ParsePunch not seen here */
                                fprintf( ioQQQ, " PROBLEM DISASTER PunchDo does not recognize flag %4.4s set by ParsePunch.  This is impossible.\n", 
                                  punch.chPunch[ipPun] );
                                TotalInsanity();
                        }

                        /* print special hash string to separate out various iterations 
                         * chTime is LAST on last iteration
                         * punch.lgHashEndIter flag is true by default, set false
                         * with "no hash" keyword on punch command
                         * punch.lg_separate_iterations is true by default, set false
                         * when punch time dependent calcs since do not want special
                         * character between time steps 
                         * grid.lgGrid is only true when doing a grid of calculations */
                        if( strcmp(chTime,"LAST") == 0 && 
                                !(iterations.lgLastIt && !grid.lgGrid ) && 
                                punch.lgHashEndIter[ipPun] &&
                                punch.lg_separate_iterations[ipPun] &&
                                !punch.lgFITS[ipPun] )
                        {
                                if( dynamics.lgStatic && strcmp( punch.chHashString , "TIME_DEP" )==0 )
                                {
                                        fprintf( punch.ipPnunit[ipPun], "\"time=%f\n",
                                                dynamics.time_elapsed );
                                }
                                else
                                {
                                        fprintf( punch.ipPnunit[ipPun], "%s\n",
                                                punch.chHashString );
                                }
                        }
                        if( punch.lgFLUSH )
                                fflush( punch.ipPnunit[ipPun] );
                }
        }
        return;
}

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

        DEBUG_ENTRY( "PunLineIntensity()" );

        /* 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( ioPUN, "**********************************************************************************************************************************\n" );

        lgEOF = false;
        while( !lgEOF )
        {
                input_readarray(chCard,&lgEOF);
                if( !lgEOF )
                {
                        fprintf( ioPUN, "%s\n", chCard );
                }
        }

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

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

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

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

        return;
}

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

/*PunchLineStuff punch optical depths or source functions for all transferred lines */
STATIC void PunchLineStuff(
  FILE * ioPUN,
  const char *chJob , 
  realnum xLimit )
{

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

        DEBUG_ENTRY( "PunchLineStuff()" );

        /*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(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
                        lgFirst = false;
                }
                else
                {
                        lgPopsFirstCall = false;
                }
        }
        else
        {
                fprintf( ioQQQ, " insane job in PunchLineStuff =%s\n", 
                  &*chJob );
                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]  )
                        {
                                /* 06 aug 28, from numLevels_max to _local. */
                                for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
                                {
                                        for( ipLo=0; ipLo <ipHi; ipLo++ )
                                        {
                                                if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
                                                        continue;

                                                ++index;
                                                pun1Line( &Transitions[ipISO][nelem][ipHi][ipLo] , ioPUN , 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++ )*/
                                        /* 06 aug 28, from numLevels_max to _local. */
                                        for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1].n+1; ipHi < iso.nLyman[ipISO]; ipHi++ )
                                        {
                                                ++index;
                                                pun1Line( &ExtraLymanLines[ipISO][nelem][ipHi] , ioPUN , 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] , ioPUN , xLimit  , index , 1. );
        }

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

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


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

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

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

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

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

        if( lgPunchOpticalDepths )
        {
                /* optical depths, no special first time, only print them */
                if( t->Emis->TauIn >= xLimit )
                {
                        /* label like "C  4" or "H  1" */
                        fprintf( ioPUN , "%2.2s%2.2s\t", 
                          elementnames.chElementSym[t->Hi->nelem-1] ,
                          elementnames.chIonStage[t->Hi->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( ioPUN , t->WLAng );
                        }
                        else
                        {
                                /* this converts energy in Rydbergs into any of the other units */
                                fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t->EnergyWN * WAVNRYD)) );
                        }
                        /* print the optical depth */
                        fprintf( ioPUN , "\t%.3f", t->Emis->TauIn  );
                        /* damping constant */
                        fprintf(ioPUN, "\t%.3e", 
                                t->Emis->dampXvel / DoppVel.doppler[ t->Hi->nelem -1 ] );
                        fprintf(ioPUN, "\n");
                }
        }
        else if( lgPopsFirstCall )
        {
                char chLbl[11];
                /* first call to line populations, print atomic parameters and indices */
                strcpy( chLbl, chLineLbl(t) );
                fprintf(ioPUN, "%li\t%s" , index , chLbl );
                /* stat weights */
                fprintf(ioPUN, "\t%.0f\t%.0f", 
                        t->Lo->g ,t->Hi->g);
                /* energy difference, gf */
                fprintf(ioPUN, "\t%.2f\t%.3e", 
                        t->EnergyWN ,t->Emis->gf);
                fprintf(ioPUN, "\n");
        }
        else
        {
                /* not first call, so do level populations and indices defined above */
                if( t->Hi->Pop > xLimit )
                {
                        fprintf(ioPUN,"%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->Lo->Pop*abundance , 
                                t->Hi->Pop*abundance );
                }
        }
}

/*PunchNewContinuum produce the 'punch new continuum' output */
STATIC void PunchNewContinuum(FILE * ioPUN, 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;
                /* we 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 */
        energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
        cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
        cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
        diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
        diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
        emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
        emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
        /*emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));
          emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));*/

        /* 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 );
        /* 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( ioPUN,"%.3e\t", AnuUnit((realnum)(energy[j]+rfield.widflx[j+ipLo]/2.) ) );
                fprintf( ioPUN,"%.3e\t", cont_incid[j] );
                fprintf( ioPUN,"%.3e\t", cont_atten[j] );
                fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
                fprintf( ioPUN,"%.3e", 
                        emis_lines_out[j]+emis_lines_in[j]/*+emis_lines_pump_out[j]+emis_lines_pump_in[j]*/ );
                fprintf( ioPUN,"\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 AGN3 hemiss, for Chapter 4, routine is below */
STATIC void AGN_Hemis(FILE *ioPUN )
{
        const int NTE = 4;
        realnum te[NTE] = {5000., 10000., 15000., 20000.};
        realnum *agn_continuum[NTE];
        double TempSave = phycon.te;
        long i , j;

        DEBUG_ENTRY( "AGN_Hemis()" );

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

                /* set the next temperature */
                /* recompute everything at this new temp */
                TempChange(te[i] , true);
                /* converge the pressure-temperature-ionization solution for this zone */
                ConvPresTempEdenIoniz();

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

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

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

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

        /* Restore temperature stored before this routine was called    */
        /* and force update */
        TempChange(TempSave , true);

        fprintf( ioQQQ, "AGN_Hemis - result of punch AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
        cdEXIT(EXIT_FAILURE);
}

/*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* ioPUN)
{
        char chCard[INPUT_LINE_LENGTH];
        bool lgEOF;
        long int i , nelem , ion;

        DEBUG_ENTRY( "punResults()" );

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

        fprintf( ioPUN, "**********************************************************************************************************************************\n" );
        lgEOF = false;
        input_readarray(chCard,&lgEOF);
        while( !lgEOF )
        {
                /* will get copy of input line image as all capital letters here */
                char chCAPS[INPUT_LINE_LENGTH];
                strcpy( chCAPS , chCard );
                caps( chCAPS );
                /* keyword HIDE means to hide the command - do not print it */
                if( !nMatch( "HIDE" , chCAPS ) )
                        fprintf( ioPUN, "%s\n", chCard );
                input_readarray(chCard,&lgEOF);
        }

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

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

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

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

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

        fprintf( ioPUN, "     \n" );

        fprintf( ioPUN, "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( ioPUN, " %10.3e", mean.xIonMeans[0][nelem][ion] );
                        /* throw line feed every 10 numbers */
                        if( nelem==9|| nelem==19 || nelem==29 )
                        {
                                fprintf( ioPUN, "\n" );
                        }
                }
        }

        fprintf( ioPUN, "END OF RESULTS\n" );
        fprintf( ioPUN, "**********************************************************************************************************************************\n" );
        return;
}

static const int 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 * ioPUN, 
  /* 4 char + null string */
  const char *chLab, 
  realnum wl, 
  double xInten, 
  /* null terminated string saying what to do */
  const char *chFunction)
{

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

        DEBUG_ENTRY( "PunResults1Line()" );

        /* 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( ioPUN, " %4.4s ", chLabSave[i] );
                                prt_wl( ioPUN, wavelength[i] );
                                fprintf( ioPUN,"\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( ioPUN, "\n" );
                        }
                        /* only put cr if we did not just put one */
                        if( strcmp(punch.chPunRltType,"array ") == 0 )                          
                                fprintf( ioPUN, " \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( ioPUN, " %4.4s", chLabSave[i] );
                                prt_wl( ioPUN, wavelength[i] );
                                fprintf( ioPUN,"\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( ioPUN, "\n" );
                        }
                        if( strcmp(punch.chPunRltType,"array ") == 0 )
                                fprintf( ioPUN, " \n" );
                }
        }
        else
        {
                fprintf( ioQQQ, " PunResults1Line called with insane request=%5.5s\n", 
                  chFunction );
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

static const int NENR_GAUNT = 37;
static const int NTE_GAUNT = 21;

/*PunchGaunts called by punch gaunts command to output gaunt factors */
STATIC void PunchGaunts(FILE* ioPUN)
{
        long int i, 
          charge,
          ite, 
          j;

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


        DEBUG_ENTRY( "PunchGaunts()" );

        /* 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 = (realnum)log10((double)charge);

                /* energy 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 = cont_gaunt_calc( 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( ioPUN, "      " );
                for( i=1; i <= NTE_GAUNT; i++ )
                {
                        fprintf( ioPUN, "\t%6.1f", ste[i-1] );
                }
                fprintf( ioPUN, "\n" );

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

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

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

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

                fprintf( ioPUN, "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( ioPUN, "\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( ioPUN, "end of charge = %li\n", charge);
                fprintf( ioPUN, "****************************\n");
        }

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

---