/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*H2_ParsePunch parse the punch h2 command */
/*H2_PunchDo punch some properties of the large H2 molecule */
/*chMolBranch returns a char with the spectroscopic branch of a transition */
/*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
/*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from PunchLineStuff */
/*H2_Punch_line_data punch line data for H2 molecule */
/*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */
/*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
/*H2_ReadEnergies read energies for all electronic levels */
/*H2_ReadTransprob read transition probabilities */
/*H2_ReadCollRates read collision rates */
/*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
/*H2_ParsePunch parse the punch h2 command */
/*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */
/*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
 /*cdH2_Line returns 1 if we found the line, 
  * or FALSE==0 if we did not find the line because ohoto-para transition
  * or upper level has lower energy than lower level */
#include "cddefines.h" 
#include "physconst.h" 
#include "punch.h" 
#include "hmi.h"
#include "path.h"
#include "prt.h"
#include "secondaries.h"
#include "grainvar.h"
#include "phycon.h"
#include "rfield.h"
#include "hyperfine.h"
#include "thermal.h"
#include "lines.h"
#include "lines_service.h"
#include "dense.h"
#include "radius.h"
#include "colden.h"
#include "taulines.h"
#include "h2_priv.h"
#include "h2.h"
#include "cddrive.h"
#include "mole.h"
/*lint -e662 possible access of out-of-bounds pointer */

/* this will say whether ortho or para,
 * H2_lgOrtho is 0 or 1 depending on whether or not ortho, 
 * so chlgPara[H2_lgOrtho] gives P or O for printing */
static char chlgPara[2]={'P','O'};

/* intensity, relative to normalization line, for faintest line to punch */
static float thresh_punline_h2;

/*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
void H2_LinesAdd(void)
{
	/* these are the quantum designations of the lines we will output */
	int iRotHi, iVibHi, iElecHi ,iRotLo, iVibLo, iElecLo;

	/* H2 not on, so space not allocated */
	if( !h2.lgH2ON )
		return;

	/* >>chng 05 nov 04, make info copies of these lines up here
	 * these are among the strongest lines in the 2 mic window and some have nearly the same
	 * wavelength as far weaker lines that may come before them in the line stack.  in that case
	 * cdLine would find the much weaker line with the same wavelength.
	 * put strong H2 lines first so that line search will find these, and not far weaker
	 * lines with nearly the same wavelength - these will be duplicated in the output but 
	 * these are here for into (the 'i) so does no harm
	 *
	 * the array indices in the following strucutures give upper and lower elec, jiv, rot quantum
	 * indices. 
	 * H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]
	  * >>chng 05 dec 22, had hand entered wavelength in A as second parameter.  This gave
	  * rounded off result when set line precision 5 was used.  now uses same logic that
	  * PutLine will eventually use - simply enter same wl in Ang 
	 * 1-0 S(4) - 18910 */
	linadd( H2Lines[0][1][6][0][0][4].xIntensity , H2Lines[0][1][6][0][0][4].WLAng , "H2  " ,'i');
	/* 1-0 S(3) - 19570 */
	linadd( H2Lines[0][1][5][0][0][3].xIntensity , H2Lines[0][1][5][0][0][3].WLAng , "H2  " ,'i');
	/* 1-0 S(2) - 20330 */
	linadd( H2Lines[0][1][4][0][0][2].xIntensity , H2Lines[0][1][4][0][0][2].WLAng , "H2  " ,'i');
	/* 1-0 S(1) - 21210 */
	linadd( H2Lines[0][1][3][0][0][1].xIntensity , H2Lines[0][1][3][0][0][1].WLAng , "H2  " ,'i');
	/* 1-0 S(0) - 22230 */
	linadd( H2Lines[0][1][2][0][0][0].xIntensity , H2Lines[0][1][2][0][0][0].WLAng , "H2  " ,'i');
	/* 1-0 Q(2) - 24130 */
	linadd( H2Lines[0][1][2][0][0][2].xIntensity , H2Lines[0][1][2][0][0][2].WLAng , "H2  " ,'i');
	/* 1-0 Q(1) - 24060 */
	linadd( H2Lines[0][1][1][0][0][1].xIntensity , H2Lines[0][1][1][0][0][1].WLAng , "H2  " ,'i');
	/* 1-0 Q(0) - 24020 */
	linadd( H2Lines[0][1][0][0][0][0].xIntensity , H2Lines[0][1][0][0][0][0].WLAng , "H2  " ,'i');

	/* print all lines from lowest n levels within X */
	/* loop over all possible lines and set H2_populations, 
	 * and quantities that depend on them */
	for( iElecHi=0; iElecHi<h2.nElecLevelOutput; ++iElecHi )
	{
		for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
		{
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{
				long int lim_elec_lo = 0;
				/* now the lower levels */
				/* NB - X is the only lower level considered here, since we are only 
				* concerned with excited electronic levels as a photodissociation process
				* code exists to relax this assumption - simply change following to iElecHi */
				for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
				{
					/* want to include all vib states in lower level if different elec level,
					* but only lower vib levels if same elec level */
					long int nv = nVib_hi[iElecLo];
					if( iElecLo==iElecHi )
						nv = iVibHi;
					for( iVibLo=0; iVibLo<=nv; ++iVibLo )
					{
						long nr = nRot_hi[iElecLo][iVibLo];
						if( iElecLo==iElecHi && iVibHi==iVibLo )
							nr = iRotHi-1;

						for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
						{
							if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
							{
								/* all ground vib state rotation lines - first is J to J-2 */
								PutLine(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]);
								if( LineSave.ipass == 0 )
								{
									H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] =  0.;
								}
								else if( LineSave.ipass == 1 )
								{
									H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] += (float)(
										radius.dVeff*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity);
								}
							}
						}
					}
				}
			}
		}
	}
	return;
}

/*H2_ParsePunch parse the punch h2 command */
void H2_ParsePunch( char *chCard )
{
	long int i , nelem ;
	int lgEOL;

	i = 5;
	nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	if( nelem != 2 )
	{
		fprintf( ioQQQ, " The first number on this line must be the 2 in H2\n Sorry.\n" );
		puts( "[Stop in ParsePunch]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* this provides info on the large H2 molecule */
	if( lgMatch("COLU",chCard) )
	{
		/* punch column density */
		strcpy( punch.chPunch[punch.npunch], "H2cl" );

		/* this is an option to scan off highest vib and rot states 
		 * to punch pops - first is limit to vibration, then rotation 
		 * if no number is entered then 0 is set and all levels punched */
		/* now get vib lim */
		punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* highest rotation */
		punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		/* this says whether to punch triplets or a matrix for output -
		 * default is triplets, so only check for matrix */
		if( lgMatch( "MATR" , chCard ) )
		{
			/* matrix */
			punch.punarg[2][punch.npunch] = 1;
			fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\tcolumn density\n" );
		}
		else
		{
			/* triplets */
			punch.punarg[2][punch.npunch] = -1;
			fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\tEner(K)\tcolden\tcolden/stat wght\tLTE colden\tLTE colden/stat wght\n" );
		}
	}
	else if( lgMatch("COOL",chCard) )
	{
		/* heating and cooling rates */
		strcpy( punch.chPunch[punch.npunch], "H2co" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#H2 depth\ttot cool\tTH Sol\tBig Sol\tTH pht dis\tpht dis\tTH Xcool\tXcool \n" );
	}

	else if( lgMatch("CREA",chCard) )
	{
		/* H2 creation rates */
		strcpy( punch.chPunch[punch.npunch], "H2cr" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#H2 depth\tH2_rate_create\tH2_rate_destroy\trate_h2_form_grains_used_total\tassoc_detach");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tbh2dis\tbh2h2p\tradasc\th3ph2p\th2phmh2h\tbh2h22hh2\th3phmh2hh\th3phm2h2\th32h2\teh3_h2h\th3ph2hp");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH_CH_C_H2\tH_CHP_CP_H2\tH_CH2_CH_H2\tH_CH3P_CH2P_H2\tH_OH_O_H2\tHminus_HCOP_CO_H2\tHminus_H3OP_H2O_H2\tHminus_H3OP_OH_H2_H");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tHP_CH2_CHP_H2\tHP_SiH_SiP_H2\tH2P_CH_CHP_H2\tH2P_CH2_CH2P_H2\tH2P_CO_COP_H2\tH2P_H2O_H2OP_H2\tH2P_O2_O2P_H2");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH2P_OH_OHP_H2\tH3P_C_CHP_H2\tH3P_CH_CH2P_H2\tH3P_CH2_CH3P_H2\tH3P_OH_H2OP_H2\tH3P_H2O_H3OP_H2\tH3P_CO_HCOP_H2");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH3P_O_OHP_H2\tH3P_SiH_SiH2P_H2\tH3P_SiO_SiOHP_H2\tH_CH3_CH2_H2\tH_CH4P_CH3P_H2\tH_CH5P_CH4P_H2\tH2P_CH4_CH3P_H2");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH2P_CH4_CH4P_H2\tH3P_CH3_CH4P_H2\tH3P_CH4_CH5P_H2\tHP_CH4_CH3P_H2\tHP_HNO_NOP_H2\tHP_HS_SP_H2\tH_HSP_SP_H2");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH3P_NH_NH2P_H2\tH3P_NH2_NH3P_H2\tH3P_NH3_NH4P_H2\tH3P_CN_HCNP_H2\tH3P_NO_HNOP_H2\tH3P_S_HSP_H2\tH3P_CS_HCSP_H2");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH3P_NO2_NOP_OH_H2\tH2P_NH_NHP_H2\tH2P_NH2_NH2P_H2\tH2P_NH3_NH3P_H2\tH2P_CN_CNP_H2\tH2P_HCN_HCNP_H2\tH2P_NO_NOP_H2");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH3P_Cl_HClP_H2\tH3P_HCl_H2ClP_H2\tH2P_C2_C2P_H2\tHminus_NH4P_NH3_H2\tH3P_HCN_HCNHP_H2");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tdestruction/creation\tav Einstein A\n");
	}
	else if( lgMatch("DEST",chCard) )
	{
		/* H2 destruction rates */
		strcpy( punch.chPunch[punch.npunch], "H2ds" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\ttot H2 rate create\ttot H2 rate destroy\ttot H- backwards\tSolomon H2g\tSolomon H2s\tphotodissoc H2s\tphotodissoc H2g");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\te- dissoc\trh2h2p\th2hph3p\tH0 dissoc\tCR\trheph2hpheh\theph2heh2p\thehph2h3phe\th3petc\tH2Ph3p");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\th2sh2g\th2h22hh2\th2sh2sh2g2h\th2sh2sh2s2h\tH2_CHP_CH2P_H\tH2_CH2P_CH3P_H\tH2_OHP_H2OP_H\tH2_H2OP_H3OP_H\tH2_COP_HCOP_H");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH2_OP_OHP_H\tH2_SiOP_SiOHP_H\tH2_C_CH_H\tH2_CP_CHP_H\tH2_CH_CH2_H\tH2_OH_H2O_H\tH2_O_OH_H");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\th2s_ch_ch2_h\th2s_o_oh_h\th2s_oh_h2o_h\th2s_c_ch_h\th2s_cp_chp_h\tH2_CH2_CH3_H\tH2_CH3_CH4_H");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH2_CH4P_CH5P_H\tH2s_CH2_CH3_H\tH2s_CH3_CH4_H\th2s_op_ohp_h\tH2_N_NH_H\tH2_NH_NH2_H\tH2_NH2_NH3_H\tH2_CN_HCN_H\tH2_NP_NHP_H");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH2_NHP_N_H3P\tH2_NHP_NH2P_H\tH2_NH2P_NH3P_H\tH2_NH3P_NH4P_H\tH2_CNP_HCNP_H\tH2_SP_HSP_H\tH2_CSP_HCSP_H");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tH2_ClP_HClP_H\tH2_HClP_H2ClP_H\tH2_HCNP_HCNHP_H");
		fprintf( punch.ipPnunit[punch.npunch], 
			"\tfrac H2g\tfrac H2s\n");
	}

	else if( lgMatch("HEAT",chCard) )
	{
		/* heating and cooling rates */
		strcpy( punch.chPunch[punch.npunch], "H2he" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#H2 depth\ttot Heat\tHeat(big)\tHeat(TH85)\tDissoc(Big)\tDissoc(TH85) \n" );
	}

	else if( lgMatch("LEVE",chCard) )
	{
		/* punch H2 level energies */
		strcpy( punch.chPunch[punch.npunch], "H2le" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#H2 energy(wn)\tstat wght\tv\tJ\n" );
	}

	else if( lgMatch("LINE",chCard) )
	{
		/* punch H2 lines - all in X */
		strcpy( punch.chPunch[punch.npunch], "H2ln" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#H2 line\tEhi\tVhi\tJhi\tElo\tVlo\tJlo\twl(mic)\twl(lab)\tlog flux\tI/Inorm\tExc Enr(K) \n" );
		/* first optional number changes the theshold of weakest line to print*/
		/* fe2thresh is intensity relative to normalization line,
		* normally Hbeta, and is set to zero in zero.c */

		/* threshold for faintest line to punch, default is 1e-4 of norm line */
		thresh_punline_h2 = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			/* this will be default relative intensity for faintest line to punch */
			thresh_punline_h2 = 1e-4f;
		}

		/* it is a log if negative */
		if( thresh_punline_h2 < 0. )
		{
			thresh_punline_h2 = (float)pow(10.f,thresh_punline_h2);
		}

		/* lines from how many electronic states?  default is one, just X, and is
		 * obtained with GROUND keyword.  ALL will produce all lines from all levels.
		 * else, if a number is present, will be the number.  if no number, no keyword,
		 * appear then just ground */
		if( lgMatch( "ELEC" , chCard ) ) 
		{
			if( lgMatch(" ALL",chCard) )
			{
				/* all electronic levels - when done, will set upper limit, the
				 * number of electronic levels actually computed, don't know this yet,
				 * so signify with negative number */
				h2.nElecLevelOutput = -1;
			}
			else if( lgMatch("GROU",chCard) )
			{
				/* just the ground electronic state */
				h2.nElecLevelOutput = 1;
			}
			else
			{
				h2.nElecLevelOutput = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
				if( lgEOL )
				{
					/* just the ground electronic state */
					h2.nElecLevelOutput = 1;
				}
			}
		}
	}

	else if( lgMatch(" PDR",chCard) )
	{
		/* creation and destruction processes */
		strcpy( punch.chPunch[punch.npunch], "H2pd" );
		fprintf( punch.ipPnunit[punch.npunch], "#H2 creation, destruction. \n" );
	}
	else if( lgMatch("POPU",chCard) )
	{
		/* punch H2_populations */
		strcpy( punch.chPunch[punch.npunch], "H2po" );

		/* this is an option to scan off highest vib and rot states 
		 * to punch pops - first is limit to vibration, then rotation 
		 * if no number is entered then 0 is set and all levels punched */
		/* now get vib lim */
		punch.punarg[0][punch.npunch] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* this is limit to rotation quantum index */
		punch.punarg[1][punch.npunch] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		if( lgMatch( "ZONE" , chCard ) )
		{
			/* punch v=0 pops for each zone, all along one line */
			punch.punarg[2][punch.npunch] = 0;
			fprintf( punch.ipPnunit[punch.npunch], "#depth\torth\tpar\te=1 rel pop\te=2 rel pop\tv=0 rel pops\n" );
		}
		else
		{
			/* will not do zone output, only output at the end of the calculation
			 * now check whether to punch triplets or a matrix for output -
			 * default is triplets, so only check for matrix */
			if( lgMatch( "MATR" , chCard ) )
			{
				/* matrix */
				punch.punarg[2][punch.npunch] = 1;
				fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\tpops\n" );
			}
			else
			{
				/* triplets */
				punch.punarg[2][punch.npunch] = -1;
				fprintf( punch.ipPnunit[punch.npunch], "#vib\trot\ts\tenergy(wn)\tpops/H2\told/H2\tpops/g\tdep coef\tFin(Col)\tFout(col)\tRCout\tRRout\tRCin\tRRin\n" );
			}
		}
	}

	else if( lgMatch("RATE",chCard) )
	{
		/* punch h2 rates - creation and destruction rates */
		strcpy( punch.chPunch[punch.npunch], "H2ra" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tN(H2)\tN(H2)/u(H2)\tA_V(star)\tn(Eval)"
			"\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tELWERT\tBigH2\telec->H2g\telec->H2s"
			"\tG(TH85)\tG(DB96)\tCR\tEleclife\tShield(BD96)\tShield(H2)\tBigh2/G0(spc)\ttot dest"
			"\tHeatH2Dish_TH85\tHeatH2Dexc_TH85\tHeatH2Dish_BigH2\tHeatH2Dexc_BigH2\thtot\n" );
	}
	else if( lgMatch("SOLO",chCard) )
	{
		/* rate of solomon process then fracs of exits from each v, J level */
		strcpy( punch.chPunch[punch.npunch], "H2so" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tSol tot\tpump/dissoc\tpump/dissoc BigH2\tavH2g\tavH2s\tH2g chem/big H2\tH2s chem/big H2\tfrac H2g BigH2\tfrac H2s BigH2\teHi\tvHi\tJHi\tvLo\tJLo\tfrac\twl(A)\n" );
	}
	else if( lgMatch("SPEC",chCard) )
	{
		/* special punch command*/
		strcpy( punch.chPunch[punch.npunch], "H2sp" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tspecial\n" );
	}
	else if( lgMatch("TEMP",chCard) )
	{
		/* various temperatures for neutral/molecular gas */
		strcpy( punch.chPunch[punch.npunch], "H2te" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tH2/H\tn(1/0)\tn(ortho/para)\tT(1/0)\tT(2/0)\tT(3/0)\tT(3/1)\tT(4/0)\tT(kin)\tT(21cm)\tT_sum(1/0)\tT_sum(2/0)\tT_sum(3/0)\tT_sum(3/1)\tT_sum(4/0) \n");
	}
	else if( lgMatch("THER",chCard) )
	{
		/* thermal heating cooling processes involving H2 */
		strcpy( punch.chPunch[punch.npunch], "H2th" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tH2/H\tn(1/0)\tn(ortho/para)\tT(1/0)\tT(2/0)\tT(3/0)\tT(3/1)\tT(4/0)\tT(kin)\tT(21cm)\tT_sum(1/0)\tT_sum(2/0)\tT_sum(3/0)\tT_sum(3/1)\tT_sum(4/0) \n");
	}
	else
	{
		fprintf( ioQQQ, 
			" There must be a second key; they are  RATE, LINE, COOL, COLUMN, _PDR, SOLOmon, TEMP, and POPUlations\n" );
		puts( "[Stop in ParsePunch]" );
		cdEXIT(EXIT_FAILURE);
	}
	return;
}


/*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
void H2_Prt_Zone(void)
{
	int iElecHi , iVibHi ;
#	ifdef DEBUG_FUN
	fputs( "<+>H2_Prt_line_tau()\n", debug_fp );
#	endif

	/* no print if H2 not turned on, or not computed for these conditions */
	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return;

	fprintf( ioQQQ, " H2 density   ");
	fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.H2_total));

	fprintf( ioQQQ, " orth/par");
	fprintf(ioQQQ,PrintEfmt("%9.2e", h2.ortho_density / SDIV( h2.para_density )));

	iElecHi = 0;
	iVibHi = 0;
	fprintf( ioQQQ, " v0 J=0,3");
	fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][0] / hmi.H2_total));
	fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][1] / hmi.H2_total));
	fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][2] / hmi.H2_total));
	fprintf(ioQQQ,PrintEfmt("%9.2e", H2_populations[iElecHi][iVibHi][3] / hmi.H2_total));

	fprintf( ioQQQ, " TOTv=0,3");
	fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total));
	fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total));
	fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total));
	fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total));
	fprintf( ioQQQ, "\n");

	return;
}

/*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */
void H2_Prt_column_density(	
	/* this is stream used for io, is stdout when called by final,
	 * is punch unit when punch output generated */
	 FILE *ioMEAN )

{
	int iVibHi ;
#	ifdef DEBUG_FUN
	fputs( "<+>H2_Prt_line_tau()\n", debug_fp );
#	endif

	/* no print if H2 not turned on, or not computed for these conditions */
	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return;

	fprintf( ioMEAN, " H2 total   ");
	fprintf(ioMEAN,"%7.3f", log10(SDIV(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])));

	fprintf( ioMEAN, " H2 orth   ");
	fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.ortho_colden)));

	fprintf( ioMEAN, " para");
	fprintf(ioMEAN,"%7.3f", log10(SDIV(h2.para_colden)));

	iVibHi = 0;
	fprintf( ioMEAN, " v0 J=0,3");
	fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][0])));
	fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][1])));
	fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][2])));
	fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][3])));

#	if 0
	fprintf( ioMEAN, "    v=0,3");
	fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][0] / hmi.H2_total));
	fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][1] / hmi.H2_total));
	fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][2] / hmi.H2_total));
	fprintf(ioMEAN,PrintEfmt("%9.2e", pops_per_vib[iElecHi][3] / hmi.H2_total));
	fprintf( ioMEAN, "\n");
#	endif

	return;
}


/*H2_ReadCollRates read collision rates */
void H2_ReadCollRates( long int nColl )
{
	/* the colliders are H, He, H2 ortho, H2 para, H+ */
	char cdDATAFILE[N_X_COLLIDER][FILENAME_PATH_LENGTH_2] = 
	{
		"H2_coll_H.dat" ,
		"H2_coll_He.dat" , 
		"H2_coll_H2ortho.dat" ,
		"H2_coll_H2para.dat",
		"H2_coll_Hp.dat" ,
		"H2_coll_He_Stancil.dat"  
	};
	FILE *ioDATA;
	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;
	long int i, n1, n2, n3;
	long int iVibHi , iVibLo , iRotHi , iRotLo;
	int lgEOL;

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

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , cdDATAFILE[nColl] );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , cdDATAFILE[nColl] );
	}

	if( nColl==5 )
	{
		/* special case, new data file from Oak Ridge project -
		 * call init routine and return - data is always read in when large H2 is included -
		 * but data are only used (for now, mid 2005) when command 
	 	 * atom H2 He OLD (Meudon) NEW (Stancil) and OFF given */
		/*H2_He_coll_init receives the name of the file that contrains the fitting coefficeints 
		 * of all transitions and read into 3d vectors. It outputs 'test.out' to test the arrays*/
		if( H2_He_coll_init( chFilename ) != 51018 )
		{
			fprintf(ioQQQ,"the H2 - He collision data file H2_coll_He_Stancil.dat does not have the correct magic number.\n");
			fprintf(ioQQQ,"I expected %i\n", 51018 );

			puts( "[Stop in H2_ReadCollRates]" );
			cdEXIT(EXIT_FAILURE);
		}

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

	/* now open the data file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadCollRates could not open %s\n", cdDATAFILE[nColl] );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " H2_ReadCollRates could not open %s\n", cdDATAFILE[nColl]);
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in H2_ReadCollRates]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read the first line and check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadCollRates could not read first line of %s\n", cdDATAFILE[nColl]);
		puts( "[Stop in H2_ReadCollRates]" );
		cdEXIT(EXIT_FAILURE);
	}
	i = 1;
	/* level 1 magic number */
	n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
	if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
	{
		fprintf( ioQQQ, 
			" H2_ReadCollRates: the version of %s is not the current version.\n", cdDATAFILE[nColl] );
		fprintf( ioQQQ, 
			" I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
			n1 , n2 , n3 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in H2_ReadCollRates]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read until not a comment */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		BadRead();

	while( chLine[0]=='#' )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
	}
	iRotHi = 1;
	while( iRotHi >= 0 )
	{
		double a[3];
		sscanf(chLine,"%li\t%li\t%li\t%li\t%le\t%le\t%le", 
			&iVibHi ,&iRotHi , &iVibLo , &iRotLo , &a[0],&a[1],&a[2] );
		/* negative iRotHi says end of data */
		if( iRotHi < 0 )
			continue;

		/* check that we actually included the levels in the model representation */
		ASSERT( iVibHi <= VIB_COLLID && 
		    iVibLo <= VIB_COLLID && 
			iRotHi <= nRot_hi[0][iVibHi] && 
			iRotLo <= nRot_hi[0][iVibLo]);

		/* some H collision rates have the same upper and lower indices - skip them */
		if( !( (iVibHi == iVibLo) && (iRotHi == iRotLo  )) )
		{
			/* this is downward transition - make sure that the energy difference is positive */
			ASSERT( (energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] ) > 0. );
			for( i=0; i<3; ++i )
			{
				CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][i] = (float)a[i];
			}

			/* this prints all levels with rates 
			fprintf(ioQQQ,"no\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n", 
				iVibHi,iRotHi,iVibLo,iRotLo,a[0],a[1],a[2]);*/
		}

		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
		while( chLine[0]=='#' )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				BadRead();
		}
	}
	fclose( ioDATA );

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

	return;
}
/* end of H2_ReadCollRates */

/*H2_ReadTransprob read transition probabilities */
void H2_ReadTransprob( long int nelec )
{
	char cdDATAFILE[N_H2_ELEC][FILENAME_PATH_LENGTH_2] = 
	{
		"H2_transprob_X.dat" ,
		"H2_transprob_B.dat" , 
		"H2_transprob_C_plus.dat" ,
		"H2_transprob_C_minus.dat" , 
		"H2_transprob_B_primed.dat" , 
		"H2_transprob_D_plus.dat" ,
		"H2_transprob_D_minus.dat" 
	};
	FILE *ioDATA;
	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;
	long int i, n1, n2, n3;
	long int iVibHi , iVibLo , iRotHi , iRotLo , iElecHi , iElecLo;
	int lgEOL;

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , cdDATAFILE[nelec] );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , cdDATAFILE[nelec] );
	}

	/* now open the data file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadTransprob could not open %s\n", cdDATAFILE[nelec] );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " H2_ReadTransprob could not open %s\n", cdDATAFILE[nelec]);
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in H2_ReadTransprob]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read the first line and check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadTransprob could not read first line of %s\n", cdDATAFILE[nelec]);
		puts( "[Stop in H2_ReadTransprob]" );
		cdEXIT(EXIT_FAILURE);
	}
	i = 1;
	/* level 1 magic number */
	n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
	if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
	{
		fprintf( ioQQQ, 
			" H2_ReadTransprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
		fprintf( ioQQQ, 
			" I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
			n1 , n2 , n3 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in H2_ReadTransprob]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read until not a comment */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		BadRead();

	while( chLine[0]=='#' )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
	}
	iVibHi = 1;
	while( iVibHi >= 0 )
	{
		double Aul;
		sscanf(chLine,"%li\t%li\t%li\t%li\t%li\t%li\t%le", 
			&iElecHi , &iVibHi ,&iRotHi , &iElecLo , &iVibLo , &iRotLo , &Aul );
		ASSERT( iElecHi == nelec );
		/* negative iVibHi says end of data */
		if( iVibHi < 0 )
			continue;

		/* check that we actually included the levels in the model representation */
		if( iVibHi <= nVib_hi[iElecHi] && 
		    iVibLo <= nVib_hi[iElecLo] && 
			iRotHi <= nRot_hi[iElecHi][iVibHi] && 
			iRotLo <= nRot_hi[iElecLo][iVibLo])
		{
			double ener = energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo];

			/* lines were initialized with Aul of zero, here overwrite with correct value */
			H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul = (float)Aul;
			/* say that this line exists */
			lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] = TRUE;
			/* prints transitions with negative energies  -  should not happen */
			if( ener <= 0. )
			{
				fprintf(ioQQQ,"negative energy H2 transition\t%li\t%li\t%li\t%li\t%.2e\t%.2e\n", 
					iVibHi,iVibLo,iRotHi,iRotLo,Aul,
					H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN);
				ShowMe();
			}
		}
#		if 0
		/* this prints all levels with As but without energies */
		else
		{
			fprintf(ioQQQ,"no\t%li\t%li\t%li\t%li\t%.2e\n", 
				iVibHi,iVibLo,iRotHi,iRotLo,Aul);
		}
#		endif

		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
		while( chLine[0]=='#' )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				BadRead();
		}
	}
	fclose( ioDATA );
	return;
}

#if 0
/*H2_Read_Cosmicray_distribution read distribution function for H2 population following cosmic ray collisional excitation */
void H2_Read_Cosmicray_distribution(void)
{
	/*>>refer	H2	cr excit	Tine, S., Lepp, S., Gredel, R., & Dalgarno, A. 1997, ApJ, 481, 282 */
	FILE *ioDATA;
	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;
	long int i, n1, n2, n3, iVib , iRot;
	long neut_frac;
	int lgEOL;

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

	/* check on path if file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "H2_CosmicRay_collision.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "H2_CosmicRay_collision.dat" );
	}

	/* now open the data file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " H2_Read_CosmicRay_distribution could not open %s\n", chFilename );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " H2_Read_CosmicRay_distribution could not open %s\n", chFilename );
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in H2_Read_Cosmicray_distribution]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read the first line and check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " H2_Read_Cosmicray_distribution could not read first line of %s\n", "H2_Cosmic_collision.dat");
		puts( "[Stop in H2_Read_Cosmicray_distribution]" );
		cdEXIT(EXIT_FAILURE);
	}

	i = 1;
	/* level 1 magic number */
	n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of H2_Cosmic_collision.dat 01 21 03 */
	if( ( n1 != 1 ) || ( n2 != 21 ) || ( n3 != 3 ) )
	{
		fprintf( ioQQQ, 
			" H2_Read_Cosmicray_distribution: the version of %s is not the current version.\n", "H2_Cosmic_collision.dat" );
		fprintf( ioQQQ, 
			" I expected to find the number 1 21 3 and got %li %li %li instead.\n" ,
			n1 , n2 , n3 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in H2_Read_Cosmicray_distribution]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read until not a comment */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		BadRead();

	while( chLine[0]=='#' )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
	}

	iRot = 1;
	iVib = 1;
	neut_frac = 0;
	while( iVib >= 0 )
	{
		long int j_minus_ji ;
		double a[10] ;
	
		sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf", 
			&iVib ,&j_minus_ji , &a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9] 
			);
		/* negative iVib says end of data */
		if( iVib < 0 )
			continue;

		/* cr_rate[CR_X][CR_VIB][CR_J][CR_EXIT];*/
		/* check that we actually included the levels in the model representation */
		ASSERT( iVib < CR_VIB );
		ASSERT( j_minus_ji == -2 || j_minus_ji == +2 || j_minus_ji == 0 );
		ASSERT( neut_frac < CR_X );

		/* now make i_minus_ji an array index */
		j_minus_ji = 1 + j_minus_ji/2;
		ASSERT( j_minus_ji>=0 && j_minus_ji<=2 );
		
		/* option to add gaussian random mole */
		for( iRot=0; iRot<CR_J; ++iRot )
		{
			cr_rate[neut_frac][iVib][iRot][j_minus_ji] = (float)a[iRot];
		}
		if( mole.lgH2_NOISECOSMIC )
		{
			float r;
			r = (float)RandGauss( mole.xMeanNoise , mole.xSTDNoise );

			for( iRot=0; iRot<CR_J; ++iRot )
			{
				cr_rate[neut_frac][iVib][iRot][j_minus_ji] *= (float)pow(10.,(double)r);
			}
		}

		if( CR_PRINT )
		{
			fprintf(ioQQQ,"cr rate\t%li\t%li", iVib , j_minus_ji ); 
			for( iRot=0; iRot<CR_J; ++iRot )
			{ 
				fprintf(ioQQQ,"\t%.3e", cr_rate[neut_frac][iVib][iRot][j_minus_ji] );
			} 
			fprintf(ioQQQ,"\n" );
		}

		/* now get next line */
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
		while( chLine[0]=='#' )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				BadRead();
		}
	}
	fclose( ioDATA );

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

	return;
}
#endif

/*H2_ReadEnergies read energies for all electronic levels */
void H2_ReadEnergies( long int nelec )
{
	char cdDATAFILE[N_H2_ELEC][FILENAME_PATH_LENGTH_2] = 
	{
		"H2_energy_X.dat" ,
		"H2_energy_B.dat" , 
		"H2_energy_C_plus.dat" ,
		"H2_energy_C_minus.dat" , 
		"H2_energy_B_primed.dat" , 
		"H2_energy_D_plus.dat" ,
		"H2_energy_D_minus.dat" 
	};
	FILE *ioDATA;
	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;
	long int i, n1, n2, n3, iVib , iRot;
	int lgEOL;

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , cdDATAFILE[nelec] );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , cdDATAFILE[nelec] );
	}

	/* now open the data file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadEnergies could not open %s\n", cdDATAFILE[nelec] );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " H2_ReadEnergies could not open %s\n", cdDATAFILE[nelec]);
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in H2_ReadEnergies]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read the first line and check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadEnergies could not read first line of %s\n", cdDATAFILE[nelec]);
		puts( "[Stop in H2_ReadEnergies]" );
		cdEXIT(EXIT_FAILURE);
	}
	i = 1;
	/* level 1 magic number */
	n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
	if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
	{
		fprintf( ioQQQ, 
			" H2_ReadEnergies: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
		fprintf( ioQQQ, 
			" I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
			n1 , n2 , n3 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in H2_ReadEnergies]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read until not a comment */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		BadRead();

	while( chLine[0]=='#' )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
	}

	/* this will count the number of levels within each elec state */
	nLevels_per_elec[nelec] = 0;

	for( iVib=0; iVib<=nVib_hi[nelec]; ++iVib )
	{
		for( iRot=Jlowest[nelec]; iRot<=nRot_hi[nelec][iVib]; ++iRot )
		{
			i = 1;
			sscanf(chLine,"%li\t%li\t%le", &n1 , &n2 , &energy_wn[nelec][iVib][iRot] );
			ASSERT( n1 == iVib );
			ASSERT( n2 == iRot );
#			if 0
			/* in atomic units, or 1 Hartree, or two rydbergs */
			if( nelec == 0 )
			{
				/* only do this for Phillip Stancil's file */
				/* corrections are to get lowest rotation level to have energy of zero */
				energy_wn[0][iVib][iRot] = -( energy_wn[0][iVib][iRot]- 3.6118114E+04 );
			}
#			endif
			ASSERT( energy_wn[nelec][iVib][iRot]> 0. || (nelec==0 && iVib==0 && iRot==0 ) );
			/* increment number of levels within this elec state */
			++nLevels_per_elec[nelec];

			/* now start reading next line */
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				BadRead();
			while( chLine[0]=='#' )
			{
				if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
					BadRead();
			}
		}
	}
	fclose( ioDATA );
	return;
}

/*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
void H2_ReadDissprob( long int nelec )
{
	char cdDATAFILE[N_H2_ELEC][FILENAME_PATH_LENGTH_2] = 
	{
		"H2_dissprob_X.dat" ,/* this does not exist and nelec == 0 is not valid */
		"H2_dissprob_B.dat" , 
		"H2_dissprob_C_plus.dat" ,
		"H2_dissprob_C_minus.dat" , 
		"H2_dissprob_B_primed.dat" , 
		"H2_dissprob_D_plus.dat" ,
		"H2_dissprob_D_minus.dat" 
	};
	FILE *ioDATA;
	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;
	long int i, n1, n2, n3, iVib , iRot;
	int lgEOL;

	ASSERT( nelec > 0 );

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , cdDATAFILE[nelec] );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , cdDATAFILE[nelec] );
	}

	/* now open the data file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadDissprob could not open %s\n", cdDATAFILE[nelec] );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " H2_ReadDissprob could not open %s\n", cdDATAFILE[nelec]);
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in H2_ReadDissprob]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read the first line and check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " H2_ReadDissprob could not read first line of %s\n", cdDATAFILE[nelec]);
		puts( "[Stop in H2_ReadDissprob]" );
		cdEXIT(EXIT_FAILURE);
	}
	i = 1;
	/* level 1 magic number */
	n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
	if( ( n1 != 3 ) || ( n2 != 2 ) || ( n3 != 11 ) )
	{
		fprintf( ioQQQ, 
			" H2_ReadDissprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
		fprintf( ioQQQ, 
			" I expected to find the number 3 2 11 and got %li %li %li instead.\n" ,
			n1 , n2 , n3 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in H2_ReadDissprob]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read until not a comment */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		BadRead();

	while( chLine[0]=='#' )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
	}

	for( iVib=0; iVib<=nVib_hi[nelec]; ++iVib )
	{
		for( iRot=Jlowest[nelec]; iRot<=nRot_hi[nelec][iVib]; ++iRot )
		{
			double a, b;
			i = 1;
			sscanf(chLine,"%li\t%li\t%le\t%le", 
				&n1 , &n2 , 
				/* dissociation probability */
				&a ,
				/* dissociation kinetic energy - eV not ergs */
				&b);

			/* these have to agree if data file is valid */
			ASSERT( n1 == iVib );
			ASSERT( n2 == iRot );

			/* dissociation probability */
			H2_dissprob[nelec][iVib][iRot] = (float)a;
			/* dissociation kinetic energy - eV not ergs */
			H2_disske[nelec][iVib][iRot] = (float)b;

			/* now get next line */
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				BadRead();
			while( chLine[0]=='#' )
			{
				if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
					BadRead();
			}
		}
	}
	fclose( ioDATA );
}


/*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */
void H2_Read_hminus_distribution(void)
{
	FILE *ioDATA;
	char chLine[FILENAME_PATH_LENGTH_2] , 
		chFilename[FILENAME_PATH_LENGTH_2] ;
	long int i, n1, n2, n3, iVib , iRot;
	int lgEOL;
	double sumrate[nTE_HMINUS];
	/* set true for lots of printout */
#	define H2HMINUS_PRT	FALSE

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "H2_hminus_deposit.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "H2_hminus_deposit.dat" );
	}

	/* now open the data file */
	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " H2_Read_hminus_distribution could not open %s\n", "H2_hminus_deposit.dat" );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " H2_Read_hminus_distribution could not open %s\n", "H2_hminus_deposit.dat");
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in H2_Read_hminus_distribution]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read the first line and check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " H2_Read_hminus_distribution could not read first line of %s\n", "H2_hminus_deposit.dat");
		puts( "[Stop in H2_Read_hminus_distribution]" );
		cdEXIT(EXIT_FAILURE);
	}

	i = 1;
	/* level 1 magic number */
	n1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	n3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* magic number
	 * the following is the set of numbers that appear at the start of H2_hminus_deposit.dat 01 08 10 */
	if( ( n1 != 2 ) || ( n2 != 10 ) || ( n3 != 17 ) )
	{
		fprintf( ioQQQ, 
			" H2_Read_hminus_distribution: the version of %s is not the current version.\n", "H2_hminus_deposit.dat" );
		fprintf( ioQQQ, 
			" I expected to find the number 2 10 17 and got %li %li %li instead.\n" ,
			n1 , n2 , n3 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in H2_Read_hminus_distribution]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read until not a comment */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		BadRead();

	while( chLine[0]=='#' )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
	}

	/* convert temps to log */
	for(i=0; i<nTE_HMINUS; ++i )
	{
		H2_te_hminus[i] = (float)log10(H2_te_hminus[i]);
		sumrate[i] = 0.;
	}

	iRot = 1;
	iVib = 1;
	while( iVib >= 0 )
	{
		/* set true to print rates */
		
		double a[nTE_HMINUS] , ener;
		sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf", 
			&iVib ,&iRot , &ener, &a[0],&a[1],&a[2] , &a[3],&a[4],&a[5] ,&a[6] 
			);
		/* negative iVib says end of data */
		if( iVib < 0 )
			continue;

		/* check that we actually included the levels in the model representation */
		ASSERT( iVib <= nVib_hi[0] && 
			iRot <= nRot_hi[0][iVib] );

		if( H2HMINUS_PRT )
			fprintf(ioQQQ,"hminusss\t%li\t%li", iVib , iRot );
		for( i=0; i<nTE_HMINUS; ++i )
		{
			H2_X_hminus_formation_distribution[i][iVib][iRot] = (float)pow(10.,-a[i]);
			sumrate[i] += H2_X_hminus_formation_distribution[i][iVib][iRot];
			if( H2HMINUS_PRT )
				fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
		}
		if( H2HMINUS_PRT )
			fprintf(ioQQQ,"\n" );

		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			BadRead();
		while( chLine[0]=='#' )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
				BadRead();
		}
	}
	fclose( ioDATA );

	if( H2HMINUS_PRT )
	{
		/* print total rate */
		fprintf(ioQQQ," total H- formation rate ");
		/* convert temps to log */
		for(i=0; i<nTE_HMINUS; ++i )
		{
			fprintf(ioQQQ,"\t%.3e" , sumrate[i]);
		}
		fprintf(ioQQQ,"\n" );
	}

	/* convert to dimensionless factors that add to unity */
	for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
	{
		for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
		{
			for(i=0; i<nTE_HMINUS; ++i )
			{
				H2_X_hminus_formation_distribution[i][iVib][iRot] /= (float)sumrate[i];
			}
		}
	}

	if( H2HMINUS_PRT )
	{
		/* print total rate */
		fprintf(ioQQQ,"  H- distribution function ");
		for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
		{
			for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
			{
				fprintf(ioQQQ,"%li\t%li", iVib , iRot );
				for(i=0; i<nTE_HMINUS; ++i )
				{
					fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
				}
				fprintf(ioQQQ,"\n" );
			}
		}
	}

	return;
}

/* ===================================================================== */
/*H2_Punch_line_data punch line data for H2 molecule */
void H2_Punch_line_data(
	/* io unit for punch */
	FILE* ioPUN ,
	/* punch all levels if true, only subset if false */
	int lgDoAll )
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;

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

	if( !h2.lgH2ON )
		return;

	if( lgDoAll )
	{
		fprintf( ioQQQ, 
			" H2_Punch_line_data ALL option not implemented in H2_Punch_line_data yet 1\n" );
		puts( "[Stop in H2_Punch_line_data]" );
		cdEXIT(EXIT_FAILURE);
	}
	else
	{

		/* punch line date, looping over all possible lines */
		for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
		{
			for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
			{
				for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
				{
					/* now the lower levels */
					/* NB - X is the only lower level considered here, since we are only 
					* concerned with excited electronic levels as a photodissociation process
					* code exists to relax this assumption - simply change following to iElecHi */
					long int lim_elec_lo = 0;
					for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
					{
						/* want to include all vib states in lower level if different elec level,
						 * but only lower vib levels if same elec level */
						long int nv = nVib_hi[iElecLo];
						if( iElecLo==iElecHi )
							nv = iVibHi;
						for( iVibLo=0; iVibLo<=nv; ++iVibLo )
						{
							long nr = nRot_hi[iElecLo][iVibLo];
							if( iElecLo==iElecHi && iVibHi==iVibLo )
								nr = iRotHi-1;

							for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
							{
								/* >>chng 03 feb 14, from !=0 to >0 */
								if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
								{
									/*TODO	1	add logic to deduce cs */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cs = 0.;
									/* print quantum indices */
									fprintf(ioPUN,"%2li %2li %2li %2li %2li %2li ",
										iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,iRotLo );
									Punch1LineData( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , ioPUN);
								}
							}
						}
					}
				}
			}
		}
		fprintf( ioPUN , "\n");
	}

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

/*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from PunchLineStuff */
void H2_PunchLineStuff( FILE * io , float xLimit  , long index)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;

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

	if( !h2.lgH2ON )
	{
#		ifdef DEBUG_FUN
		fputs( " <->H2_PunchLineStuff()\n", debug_fp );
#		endif
		return;
	}

	/* loop over all possible lines */
	for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
	{
		for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
		{
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{
				/* now the lower levels */
				/* NB - X is the only lower level considered here, since we are only 
				* concerned with excited electronic levels as a photodissociation process
				* code exists to relax this assumption - simply change following to iElecHi */
				long int lim_elec_lo = 0;
				for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
				{
					/* want to include all vib states in lower level if different elec level,
					* but only lower vib levels if same elec level */
					long int nv = nVib_hi[iElecLo];
					if( iElecLo==iElecHi )
						nv = iVibHi;
					for( iVibLo=0; iVibLo<=nv; ++iVibLo )
					{
						long nr = nRot_hi[iElecLo][iVibLo];
						if( iElecLo==iElecHi && iVibHi==iVibLo )
							nr = iRotHi-1;

						for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
						{
							/*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul != 0. )*/
							/* >>chng 03 feb 14, from !=0 to >0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								pun1Line( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , io , xLimit  , index , 1.);
							}
						}
					}
				}
			}
		}
	}

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

	return;

}


/*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
void H2_Prt_line_tau(void)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;

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

	if( !h2.lgH2ON )
		return;

	/* loop over all possible lines */
	for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
	{
		for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
		{
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{
				/* now the lower levels */
				/* NB - X is the only lower level considered here, since we are only 
				* concerned with excited electronic levels as a photodissociation process
				* code exists to relax this assumption - simply change following to iElecHi */
				long int lim_elec_lo = 0;
				for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
				{
					/* want to include all vib states in lower level if different elec level,
					* but only lower vib levels if same elec level */
					long int nv = nVib_hi[iElecLo];
					if( iElecLo==iElecHi )
						nv = iVibHi;
					for( iVibLo=0; iVibLo<=nv; ++iVibLo )
					{
						long nr = nRot_hi[iElecLo][iVibLo];
						if( iElecLo==iElecHi && iVibHi==iVibLo )
							nr = iRotHi-1;

						for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
						{
							/* >>chng 03 feb 14, from !=0 to >0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								prme(" c",&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
							}
						}
					}
				}
			}
		}
	}

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

	return;
}


/*chMolBranch returns a char with the spectroscopic branch of a transition */
static char chMolBranch( long iRotHi , long int iRotLo )
{
	/* these are the spectroscopic branches */
	char chBranch[5] = {'O','P','Q','R','S'};
	/* this is the index within the chBranch array */
	int ip = 2 + (iRotHi - iRotLo);
	if( ip<0 || ip>=5 )
	{
		fprintf(ioQQQ," chMolBranch called with insane iRotHi=%li iRotLo=%li ip=%i\n",
			iRotHi , iRotLo , ip );
		ip = 0;
	}

	return( chBranch[ip] );
}

/*H2_PunchDo punch some properties of the large H2 molecule */
void H2_PunchDo( FILE* io ,  char chJOB[] , const char chTime[] , long int ipPun )
{
	long int iVibHi , iElecHi , iRotHi , iVibLo , iElecLo , iRotLo,
		ip;
	long int ilo , iRot , iVib;
	long int LimVib , LimRot;

	/* which job are we supposed to do? This routine is active even when H2 is not turned on
	 * so do not test on h2.lgH2ON initially */

	/* H2_populations in computed last zone - 
	 * give all of molecule in either matrix or triplet format */
	if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") == 0) &&
		(punch.punarg[2][ipPun] != 0) )
	{
		/* >>chng 04 feb 19, do not punch if H2 not yet evaluated */
		if( h2.lgH2ON  && hmi.lgBigH2_evaluated )
		{
			iVibHi= 0;
			iRotHi = 0;
			iElecHi=0;
			/* the limit to the number of vibration levels punched -
			* default is all, but first two numbers on punch h2 pops command
			* reset limit */
			/* this is limit to vibration */
			if( punch.punarg[0][ipPun] > 0 )
			{
				LimVib = (long)punch.punarg[0][ipPun];
			}
			else
			{
				LimVib = nVib_hi[iElecHi];
			}

			/* first punch the current ortho, para, and total H2 density */
			fprintf(io,"%i\t%i\t%.3e\tortho\n", 
				103 , 
				103 ,
				h2.ortho_density );
			fprintf(io,"%i\t%i\t%.3e\tpara\n", 
				101 , 
				101 ,
				h2.para_density );
			fprintf(io,"%i\t%i\t%.3e\ttotal\n", 
				0 , 
				0 ,
				hmi.H2_total );

			/* now punch the actual H2_populations, first part both matrix and triplets */
			for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
			{
				/* this is limit to rotation quantum index */
				if( punch.punarg[1][ipPun] > 0 )
				{
					LimRot = (long)punch.punarg[1][ipPun];
				}
				else
				{
					LimRot = nRot_hi[iElecHi][iVibHi];
				}
				if( punch.punarg[2][ipPun] > 0 )
				{
					long int i;
					/* this option punch matrix */
					if( iVibHi == 0 )
					{
						fprintf(io,"vib\\rot");
						/* this is first vib, so make row of rot numbs */
						for( i=0; i<=LimRot; ++i )
						{
							fprintf(io,"\t%li",i);
						}
						fprintf(io,"\n");
					}
					fprintf(io,"%li",iVibHi );
					for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
					{
						fprintf(io,"\t%.3e", 
							H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total );
					}
					fprintf(io,"\n" );
				}
				else if( punch.punarg[2][ipPun] < 0 )
				{
					/* this option punch triplets - the default */
					for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
					{
						fprintf(io,"%li\t%li\t%c\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
							/* upper vibration and rotation quantum numbers */
							iVibHi , iRotHi ,
							/* an 'O' or 'P' for ortho or para */
							chlgPara[H2_lgOrtho[iElecHi][iVibHi][iRotHi]],
							/* the level excitation energy in wavenumbers */
							energy_wn[iElecHi][iVibHi][iRotHi],
							/* actual population relative to total H2 */
							H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ,
							/* old level H2_populations for comparison */
							H2_old_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total ,
							/* H2_populations per h2 and per statistical weight */
							H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total/H2_stat[iElecHi][iVibHi][iRotHi] ,
							/* lte departure coefficient */
							/* >>chng 05 jan 26, missing factor of H2 abundance LTE is norm to unity, not tot abund */
							H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi]*hmi.H2_total ) ,
							/* fraction of exits that were collisional */
							H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[0][iVibHi][iRotHi]) ,
							/* fraction of entrys that were collisional */
							H2_col_rate_in[iVibHi][iRotHi]/SDIV(H2_col_rate_in[iVibHi][iRotHi]+H2_rad_rate_in[iVibHi][iRotHi]),
							/* collisions out */
							H2_col_rate_out[iVibHi][iRotHi],
							/* radiation out */
							H2_rad_rate_out[0][iVibHi][iRotHi] ,
							/* radiation out */
							H2_col_rate_in[iVibHi][iRotHi],
							/* radiation in */
							H2_rad_rate_in[iVibHi][iRotHi]
							) ;
					}
				}
			}
		}
	}
	/* punch H2 populations for each zone 
	 * H2_populations of v=0 for each zone */
	else if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") != 0) &&
		(punch.punarg[2][ipPun] == 0) )
	{
		/* >>chng 04 feb 19, do not punch if h2 not yet evaluated */
		if( h2.lgH2ON  && hmi.lgBigH2_evaluated )
		{
			fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth_mid_zone , 
				h2.ortho_density , h2.para_density);
			/* rel pops of first two excited electronic states */
			fprintf(io,"\t%.3e\t%.3e", 
				pops_per_elec[1] , pops_per_elec[2]);
			iElecHi = 0;
			iVibHi = 0;
			/* this is limit to vibration quantum index */
			if( punch.punarg[0][ipPun] > 0 )
			{
				LimVib = (long)punch.punarg[1][ipPun];
			}
			else
			{
				LimVib = nRot_hi[iElecHi][iVibHi];
			}
			LimVib = MIN2( LimVib , nVib_hi[iElecHi] );
			/* this is limit to rotation quantum index */
			if( punch.punarg[1][ipPun] > 0 )
			{
				LimRot = (long)punch.punarg[1][ipPun];
			}
			else
			{
				LimRot = nRot_hi[iElecHi][iVibHi];
			}
			for( iVibHi = 0; iVibHi<=LimVib; ++iVibHi )
			{
				long int LimRotVib = MIN2( LimRot , nRot_hi[iElecHi][iVibHi] );
				for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRotVib; ++iRotHi )
				{
					fprintf(io,"\t%.3e", 
						H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total );
				}
				fprintf(io,"\t");
			}
			fprintf(io,"\n");
		}
	}

	/* punch column densities */
	else if( (strcmp( chJOB , "H2cl" ) == 0) && (strcmp(chTime,"LAST") == 0) )
	{
		iVibHi= 0;
		iRotHi = 0;
		iElecHi=0;
		/* the limit to the number of vibration levels punched -
		 * default is all, but first two numbers on punch h2 pops command
		 * reset limit */
		/* this is limit to vibration */
		if( punch.punarg[0][ipPun] > 0 )
		{
			LimVib = (long)punch.punarg[0][ipPun];
		}
		else
		{
			LimVib = nVib_hi[iElecHi];
		}

		/* first punch ortho and para H2_populations */
		fprintf(io,"%i\t%i\t%.3e\tortho\n", 
			103 , 
			103 ,
			h2.ortho_colden );
		fprintf(io,"%i\t%i\t%.3e\tpara\n", 
			101 , 
			101 ,
			h2.para_colden );
		/* total H2 column density */
		fprintf(io,"%i\t%i\t%.3e\ttotal\n", 
			0 , 
			0 ,
			colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]);

		/* punch level column densities */
		for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
		{
		if( h2.lgH2ON )
		{
			/* this is limit to rotation quantum index */
			if( punch.punarg[1][ipPun] > 0 )
			{
				LimRot = (long)punch.punarg[1][ipPun];
			}
			else
			{
				LimRot = nRot_hi[iElecHi][iVibHi];
			}
			if( punch.punarg[2][ipPun] > 0 )
			{
				long int i;
				/* punch matrix */
				if( iVibHi == 0 )
				{
					fprintf(io,"vib\\rot");
					/* this is first vib, so make row of rot numbs */
					for( i=0; i<=LimRot; ++i )
					{
						fprintf(io,"\t%li",i);
					}
					fprintf(io,"\n");
				}
				fprintf(io,"%li",iVibHi );
				for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
				{
					fprintf(io,"\t%.3e", 
						H2_X_colden[iVibHi][iRotHi]/hmi.H2_total );
				}
				fprintf(io,"\n" );
			}
			else
			{
				/* punch triplets - the default */
				for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
				{
					fprintf(io,"%li\t%li\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\n", 
						iVibHi , 
						iRotHi ,
						/* energy relative to 0,0, T1CM converts wavenumber to K */
						energy_wn[iElecHi][iVibHi][iRotHi]*T1CM,
						/* these are column densities for actual molecule */
						H2_X_colden[iVibHi][iRotHi] ,
						H2_X_colden[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi] ,
						/* these are same column densities but for LTE populations */
						H2_X_colden_LTE[iVibHi][iRotHi] ,
						H2_X_colden_LTE[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi]);
				}
			}
		}
		}
	}
	else if( (strcmp(chJOB , "H2pd" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* punch pdr 
		 * output some pdr information (densities, rates) for each zone */
		fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
			/* depth in cm */
			radius.depth_mid_zone ,
			/* the computed ortho and para densities */
			h2.ortho_density , 
			h2.para_density ,
			/* the Lyman Werner band dissociation, Tielens & Hollenbach */
			hmi.H2_Solomon_dissoc_rate_TH85_H2g , 
			/* the Lyman Werner band dissociation, Bertoldi & Draine */
			hmi.H2_Solomon_dissoc_rate_BD96_H2g,
			/* the Lyman Werner band dissociation, big H2 mole */
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g);
	}
	else if( (strcmp(chJOB , "H2co" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* punch H2 cooling - do heating cooling for each zone old new H2 */
		fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
			/* depth in cm */
			radius.depth_mid_zone ,
			/* total cooling, equal to total heating */
			thermal.ctot , 
			/* H2 destruction by Solomon process, TH85 rate */
			hmi.H2_Solomon_dissoc_rate_TH85_H2g,
			/* H2 destruction by Solomon process, big H2 model rate */
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g +
				hmi.H2_Solomon_dissoc_rate_BigH2_H2s,
			/* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
			hmi.HeatH2Dish_TH85,
			/* heating due to dissociation of elec excited states */
			hmi.HeatH2Dish_BigH2 , 
			/* cooling (usually neg and so heating) due to collisions within X */
			hmi.HeatH2Dexc_TH85,
			hmi.HeatH2Dexc_BigH2 
			);

	}
	else if( (strcmp(chJOB , "H2cr" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* PUNCH H2 CREATION - show H2 creation processes for each zone */
		/* >>chng 05 jul 15, TE, punch all H2 creation processes, unit cm-3 s-1 */
		fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e", 
			/* depth in cm */
			radius.depth_mid_zone ,
			/* creation cm-3 s-1, destruction rate, s-1 */
			hmi.H2_rate_create, 
			hmi.H2_rate_destroy, 
			gv.rate_h2_form_grains_used_total * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 
			hmi.assoc_detach * hmi.Hmolec[ipMH]*hmi.Hmolec[ipMHm] / hmi.H2_rate_create, 
			hmi.bh2dis * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 
			hmi.bh2h2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 
			hmi.radasc * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create, 
			hmi.h3ph2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create, 
			hmi.h2phmh2h * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
			hmi.bh2h22hh2 * 2 * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2g] / hmi.H2_rate_create,
			hmi.h3phmh2hh * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
			hmi.h3phm2h2 * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
			hmi.h32h2 * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			hmi.eh3_h2h * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			hmi.h3ph2hp * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create
			);
		fprintf(io,"\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e", 
		/*chemical network*/
			/*light elements*/
			co.H_CH_C_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H_CHP_CP_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H_CH2_CH_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H_CH3P_CH2P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H_OH_O_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.Hminus_HCOP_CO_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
			co.Hminus_H3OP_H2O_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
			co.Hminus_H3OP_OH_H2_H * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
			co.HP_CH2_CHP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
			co.HP_SiH_SiP_H2* hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
			co.H2P_CH_CHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_CH2_CH2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_CO_COP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_H2O_H2OP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_O2_O2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_OH_OHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H3P_C_CHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_CH_CH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_CH2_CH3P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_OH_H2OP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_H2O_H3OP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_CO_HCOP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_O_OHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_SiH_SiH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_SiO_SiOHP_H2	 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H_CH3_CH2_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H_CH4P_CH3P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H_CH5P_CH4P_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H2P_CH4_CH3P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_CH4_CH4P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H3P_CH3_CH4P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_CH4_CH5P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.HP_CH4_CH3P_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
			/* heavy elements */
			co.HP_HNO_NOP_H2  * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
			co.HP_HS_SP_H2 * hmi.Hmolec[ipMHp] / hmi.H2_rate_create,
			co.H_HSP_SP_H2 * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_create,
			co.H3P_NH_NH2P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_NH2_NH3P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_NH3_NH4P_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_CN_HCNP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_NO_HNOP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_S_HSP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_CS_HCSP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_NO2_NOP_OH_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H2P_NH_NHP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_NH2_NH2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_NH3_NH3P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_CN_CNP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_HCN_HCNP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.H2P_NO_NOP_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create, 
			co.H3P_Cl_HClP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H3P_HCl_H2ClP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create,
			co.H2P_C2_C2P_H2 * hmi.Hmolec[ipMH2p] / hmi.H2_rate_create,
			co.Hminus_NH4P_NH3_H2 * hmi.Hmolec[ipMHm] / hmi.H2_rate_create,
			co.H3P_HCN_HCNHP_H2 * hmi.Hmolec[ipMH3p] / hmi.H2_rate_create
			);
		fprintf(io,"\t%.3e\t%.3e\n",
			hmi.H2_rate_destroy  * hmi.H2_total / hmi.H2_rate_create,
			hmi.h2s_sp_decay
			);
	}
	else if( (strcmp(chJOB , "H2ds" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* punch H2 destruction - show H2 destruction processes for each zone 
		 * >>chng 05 nov 17, TE, added the new reaction H2s + O+ -> OH+ + H 
		 * >>chng 05 oct 04, TE, remove eh2hhm(was double) and include dissociation by electrons, H2g/H2s + e -> 2H  
		 * >>chng 05 jul 15, TE, punch all H2 destruction rates, weighted by fractional abundance of H2g, H2s, unit s-1 */
		fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e", 
			/* depth in cm */
			radius.depth_mid_zone ,
			hmi.H2_rate_create, 
			/* destruction rate, s-1 */
			hmi.H2_rate_destroy , 
			/* H2 + e -> H- + H0 */
			(hmi.assoc_detach_backwards_grnd * hmi.Hmolec[ipMH2g] + hmi.assoc_detach_backwards_exct * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy  / hmi.H2_total, 
			/*photons*/
			hmi.H2_Solomon_dissoc_rate_used_H2g / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 
			hmi.H2_Solomon_dissoc_rate_used_H2s / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			hmi.H2_photodissoc_used_H2s / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total, 
			hmi.H2_photodissoc_used_H2g / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 
			/*electrons*/
			(hmi.h2ge2h * hmi.Hmolec[ipMH2g] + hmi.h2se2h * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy  / hmi.H2_total,
			/*H+*/
			hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total, 
			hmi.h2hph3p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			/*H0*/
			(hmi.rh2dis * hmi.Hmolec[ipMH2g] + hmi.h2sh * hmi.Hmolec[ipMH2s]) * dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_destroy  / hmi.H2_total,
			/*CR*/
			(hmi.CR_reac_H2g * hmi.Hmolec[ipMH2g] + hmi.CR_reac_H2s * hmi.Hmolec[ipMH2s]) / hmi.H2_rate_destroy / hmi.H2_total,
			/*He+,HeH+*/
			hmi.rheph2hpheh*dense.xIonDense[ipHELIUM][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			hmi.heph2heh2p*dense.xIonDense[ipHELIUM][1] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			hmi.hehph2h3phe*hmi.Hmolec[ipMHeHp] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			/*H3+*/
			hmi.h3petc*hmi.Hmolec[ipMH3p] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			/*H2+*/
			hmi.h2ph3p*hmi.Hmolec[ipMH2p] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			/*H2s+H2g -> H2g + 2H*/
			hmi.h2sh2g * hmi.Hmolec[ipMH2g] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			/*H2g+H2g -> H2g + 2H*/
			hmi.h2h22hh2*2*hmi.Hmolec[ipMH2g] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			/*H2s+H2s -> H2g + 2H*/
			hmi.h2sh2sh2g2h*2*hmi.Hmolec[ipMH2s] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			/*H2s+H2s -> H2s + 2H*/
			hmi.h2sh2sh2s2h*2*hmi.Hmolec[ipMH2s] / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			/*chemical network*/
				/*light elements*/
			co.H2_CHP_CH2P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CH2P_CH3P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_OHP_H2OP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_H2OP_H3OP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_COP_HCOP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_OP_OHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_SiOP_SiOHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_C_CH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CP_CHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CH_CH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_OH_H2O_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_O_OH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.h2s_ch_ch2_h /  hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			co.h2s_o_oh_h /  hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			co.h2s_oh_h2o_h /  hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			co.h2s_c_ch_h /  hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,	
			co.h2s_cp_chp_h /  hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			co.H2_CH2_CH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CH3_CH4_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CH4P_CH5P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2s_CH2_CH3_H /  hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			co.H2s_CH3_CH4_H /  hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			co.h2s_op_ohp_h / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2s] / hmi.H2_total,
			/*heavy elements*/
			co.H2_N_NH_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_NH_NH2_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_NH2_NH3_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CN_HCN_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_NP_NHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_NHP_N_H3P / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_NHP_NH2P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_NH2P_NH3P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_NH3P_NH4P_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CNP_HCNP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_SP_HSP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_CSP_HCSP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_ClP_HClP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_HClP_H2ClP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total,
			co.H2_HCNP_HCNHP_H / hmi.H2_rate_destroy * hmi.Hmolec[ipMH2g] / hmi.H2_total
			);
		fprintf(io,"\t%.4e\t%.4e\n",
			/*H2g/Htot, H2s/Htot chemical network and from big molecule model*/
			hmi.Hmolec[ipMH2g] / hmi.H2_total,
			hmi.Hmolec[ipMH2s] / hmi.H2_total
			);

	}
	else if( (strcmp(chJOB , "H2le" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		for( ilo=0; ilo < nLevels_per_elec[0]; ilo++ )
		{
			ip = H2_ipX_ener_sort[ilo];
			iRot = ipRot_H2_energy_sort[ip];
			iVib = ipVib_H2_energy_sort[ip];

			/* punch H2 level energies */
			fprintf(io,"%.2f\t%li\t%li\t%li\n", 
				energy_wn[0][iVib][iRot],
				(long)H2_stat[0][iVib][iRot],
				iVib , iRot );
		}

	}
	else if( (strcmp(chJOB , "H2ra" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* punch h2 rates - some rates and lifetimes */
		double sumpop = 0. , sumlife = 0.;

		/* this block, find lifetime against photo excitation into excited elec states */
		iElecLo = 0;
		iVibLo = 0;
		if( h2.lgH2ON && hmi.lgBigH2_evaluated )
		{
			for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
			{
				for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
				{
					for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
					{
						for( iRotLo=Jlowest[iElecLo]; iRotLo<=nRot_hi[iElecLo][iVibLo]; ++iRotLo )
						{
							/* >>chng 03 feb 14, from !=0 to >0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								sumlife +=
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump *
									H2_populations[iElecLo][iVibLo][iRotLo];
								sumpop +=
									H2_populations[iElecLo][iVibLo][iRotLo];
							}
						}
					}
				}
			}
		}

		/* continue output from punch h2 rates command */
		/* find photoexcitation rates from v=0 */
		/* PDR information for each zone */
		fprintf(io,
			"%.5e\t%.3e\t%.3e\t%.3e\t%li", 
			/* depth in cm */
			radius.depth_mid_zone ,
			/* the column density (cm^-2) in H2 */
			colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
			/* this is a special form of column density - should be proportional
			 * to total shielding */
			colden.coldenH2_ov_vel ,
			/* visual extinction due to dust alone, of point source (star)*/
			rfield.extin_mag_V_point,
			/* number of large molecule evaluations in this zone */
			h2.nCallH2_this_zone );
		fprintf(io,
			"\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
			/* total H2 fraction */
			hmi.H2_total/dense.gas_phase[ipHYDROGEN] ,
			/* chemistry renorm factor */
			H2_renorm_chemistry,
			/* rate H2 forms on grains */
			gv.rate_h2_form_grains_used_total , 
			/* rate H2 forms by H minus route */
			hmi.Hmolec[ipMHm]*1.35e-9,
			/* H2 destruction by Solomon process, TH85 rate */
			hmi.H2_Solomon_dissoc_rate_TH85_H2g + hmi.H2_Solomon_dissoc_rate_TH85_H2s,
			/* H2 destruction by Solomon process, Bertoldi & Draine rate */
			hmi.H2_Solomon_dissoc_rate_BD96_H2g + hmi.H2_Solomon_dissoc_rate_BD96_H2s,
			/* H2 destruction by Solomon process, Elwert et al. in preparation */
			hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + hmi.H2_Solomon_dissoc_rate_ELWERT_H2g,
			/* H2 destruction by Solomon process, big H2 model rate */
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g + hmi.H2_Solomon_dissoc_rate_BigH2_H2s,
			/* rate s-1 H2 elec excit states decay into H2g */
			hmi.H2_Solomon_elec_decay_H2g ,
			/* rate s-1 H2 elec excit states decay into H2s */
			hmi.H2_Solomon_elec_decay_H2s 
			);
		fprintf(io,
			"\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e", 
			/* The TH85 estimate of the radiation field relative to the Habing value */
			hmi.UV_Cont_rel2_Habing_TH85_depth,
			/* The DB96 estimate of the radiation field relative to the Habing value */
			hmi.UV_Cont_rel2_Draine_DB96_depth,
			/* cosmic ray ionization rate */
			secondaries.csupra[ipHYDROGEN][0]*0.93,
			sumlife/SDIV( sumpop ) ,
			hmi.H2_Solomon_dissoc_rate_BD96_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth) ,
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth),
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.UV_Cont_rel2_Habing_spec_depth),
			hmi.H2_rate_destroy);
		fprintf(io,
			"\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			hmi.HeatH2Dish_TH85,
			hmi.HeatH2Dexc_TH85,
			hmi.HeatH2Dish_BigH2, 
			hmi.HeatH2Dexc_BigH2,
			thermal.htot);
	}
	/* punch h2 solomon */
	else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* remember as many as NSOL lines contributing to total Solomon process */
#		define NSOL 100
		double sum, one;
		long int jlosave[NSOL] , ivlosave[NSOL],
			iehisave[NSOL] ,jhisave[NSOL] , ivhisave[NSOL],
			nsave,
			ipOrdered[NSOL];
		int lgFail;
		int i;
		float fsave[NSOL], wlsave[NSOL];
		/* solomon process, and where it came from */
		fprintf(io,"%.5e\t%.3e", 
			/* depth in cm */
			radius.depth_mid_zone ,
			/* H2 destruction by Solomon process, big H2 model rate */
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g +
				hmi.H2_Solomon_dissoc_rate_BigH2_H2s);
		sum = 0.;
		iElecLo = 0;
		/* find sum of all radiative exits from X into excited electronic states */
		if( h2.lgH2ON && hmi.lgBigH2_evaluated )
		{
			for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
			{
				for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
				{
					for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
					{
						long int nv = nVib_hi[iElecLo];
						for( iVibLo=0; iVibLo<=nv; ++iVibLo )
						{
							long nr = nRot_hi[iElecLo][iVibLo];
							for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
							{
								if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
								{
									one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo *
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump;
									sum += one;
								}
							}
						}
					}
				}
			}

			/* make sure it is safe to div by sum */
			sum = SDIV( sum );
			nsave = 0;
			/* now loop back over X and print all those which contribute more than FRAC of the total */
#			define FRAC	0.01
			for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
			{
				for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
				{
					for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
					{
						long int nv = nVib_hi[iElecLo];
						for( iVibLo=0; iVibLo<=nv; ++iVibLo )
						{
							long nr = nRot_hi[iElecLo][iVibLo];
							for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
							{
								if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
								{
									one = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo *
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump;
									if( one/sum > FRAC && nsave<NSOL)
									{
										fsave[nsave] = (float)(one/sum);
										jlosave[nsave] = iRotLo;
										ivlosave[nsave] = iVibLo;
										jhisave[nsave] = iRotHi;
										ivhisave[nsave] = iVibHi;
										iehisave[nsave] = iElecHi;
										wlsave[nsave] = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng;
										++nsave;
										/*fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f", 
											iElecHi,iVibHi,iRotHi,iVibLo , iRotLo , one/sum );*/
									}
								}
							}
						}
					}
				}
			}/* iElecHi */
			/* now sort these into decreasing order */

			/* now sort by decreasing importance */
			/*spsort netlib routine to sort array returning sorted indices */
			spsort(
				/* input array to be sorted */
				fsave, 
				/* number of values in x */
				nsave, 
				/* permutation output array */
				ipOrdered, 
				/* flag saying what to do - 1 sorts into increasing order, not changing
				* the original routine */
				-1, 
				/* error condition, should be 0 */
				&lgFail);

			/* print ratio of pumps to dissociations - this is 9:1 in TH85 */
			/*>>chng 05 jul 20, TE, punch average energy in H2s and renormalization factors for H2g and H2s */
			/* >>chng 05 sep 16, TE, chng denominator to do g and s with proper dissoc rates */
			fprintf(io,"\t%.3f\t%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e",
				/* this is sum of photons and CRs */
				(sum + secondaries.csupra[ipHYDROGEN][0]*2.02f)/SDIV((hmi.H2_Solomon_dissoc_rate_BigH2_H2g * hmi.Hmolec[ipMH2g] +
					hmi.H2_Solomon_dissoc_rate_BigH2_H2s * hmi.Hmolec[ipMH2s]) ), 
				/* this is sum of photons and CRs */
				(sum + secondaries.csupra[ipHYDROGEN][0]*2.02f) /SDIV((hmi.H2_Solomon_dissoc_rate_BigH2_H2g *  hmi.H2g_BigH2 +
					hmi.H2_Solomon_dissoc_rate_BigH2_H2s * hmi.H2s_BigH2) ),
				hmi.H2_BigH2_H2g_av, hmi.H2_BigH2_H2s_av, 
				hmi.H2_chem_BigH2_H2g, hmi.H2_chem_BigH2_H2s,
				hmi.H2g_BigH2/SDIV(hmi.H2_total_BigH2), hmi.H2s_BigH2/SDIV(hmi.H2_total_BigH2)
				);
			for( i=0; i<nsave; ++i )
			{
				ip = ipOrdered[i];
				/*lint -e644 not init */
				fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f\t%.3f", 
					iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[ip] , wlsave[ip] );
				/*lint +e644 not init */
			}
			fprintf(io,"\n"); 
		}
		/*fprintf(io,"DEBUG tau\t%.3e\t%.3f\n",
			H2Lines[1][0][1][0][0][0].TauIn, H2Lines[1][0][1][0][0][0].WLAng); */
#		undef NSOL
	}

	else if( (strcmp(chJOB , "H2te" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* punch h2 temperatures */
		double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40;
		double T10,T20,T30,T31,T40;
		/* subscript"sum" denotes integrated quantities */
		double T10_sum,T20_sum,T30_sum,T31_sum,T40_sum;
		double pop_ratio10_sum,pop_ratio20_sum,pop_ratio30_sum,pop_ratio31_sum,pop_ratio40_sum;
		if( h2.lgH2ON && h2.nCallH2_this_zone )
		{
			double energyK = T1CM*(energy_wn[0][0][1] - energy_wn[0][0][0]);
			/* the ratio of H2_populations of J=1 to 0 */
			pop_ratio10 = H2_populations[0][0][1]/SDIV(H2_populations[0][0][0]);
			pop_ratio10_sum = H2_X_colden[0][1]/SDIV(H2_X_colden[0][0]);
			/* the corresponding temperature */
			T10 = -170.5/log(SDIV(pop_ratio10) * H2_stat[0][0][0]/H2_stat[0][0][1]);
			T10_sum = -170.5/log(SDIV(pop_ratio10_sum) * H2_stat[0][0][0]/H2_stat[0][0][1]);

			energyK = T1CM*(energy_wn[0][0][2] - energy_wn[0][0][0]);
			pop_ratio20 = H2_populations[0][0][2]/SDIV(H2_populations[0][0][0]);
			T20 = -energyK/log(SDIV(pop_ratio20) * H2_stat[0][0][0]/H2_stat[0][0][2]);

			pop_ratio20_sum = H2_X_colden[0][2]/SDIV(H2_X_colden[0][0]);
			T20_sum = -energyK/log(SDIV(pop_ratio20_sum) * H2_stat[0][0][0]/H2_stat[0][0][2]);

			energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][0]);
			pop_ratio30 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][0]);
			T30 = -energyK/log(SDIV(pop_ratio30) * H2_stat[0][0][0]/H2_stat[0][0][3]);

			pop_ratio30_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][0]);
			T30_sum = -energyK/log(SDIV(pop_ratio30_sum) * H2_stat[0][0][0]/H2_stat[0][0][3]);

			energyK = T1CM*(energy_wn[0][0][3] - energy_wn[0][0][1]);
			pop_ratio31 = H2_populations[0][0][3]/SDIV(H2_populations[0][0][1]);
			T31 = -energyK/log(SDIV(pop_ratio31) * H2_stat[0][0][1]/H2_stat[0][0][3]);

			pop_ratio31_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][1]);
			T31_sum = -energyK/log(SDIV(pop_ratio31_sum) * H2_stat[0][0][1]/H2_stat[0][0][3]);

			energyK = T1CM*(energy_wn[0][0][4] - energy_wn[0][0][0]);
			pop_ratio40 = H2_populations[0][0][4]/SDIV(H2_populations[0][0][0]);
			T40 = -energyK/log(SDIV(pop_ratio40) * H2_stat[0][0][0]/H2_stat[0][0][4]);

			pop_ratio40_sum = H2_X_colden[0][4]/SDIV(H2_X_colden[0][0]);
			T40_sum = -energyK/log(SDIV(pop_ratio40_sum) * H2_stat[0][0][0]/H2_stat[0][0][4]);
		}
		else
		{
			pop_ratio10 = 0.;
			pop_ratio10_sum = 0.;
			T10 = 0.;
			T20 = 0.;
			T30 = 0.;
			T31 = 0.;
			T40 = 0.;
			T10_sum = 0.;
			T20_sum = 0.;
			T30_sum = 0.;
			T31_sum = 0.;
			T40_sum = 0.;
		}

		/* various temperatures for neutral/molecular gas */
		fprintf( io, 
			"%.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\n" ,
			/* depth in cm */
			radius.depth_mid_zone ,
			/* total H2 fraction */
			hmi.H2_total/dense.gas_phase[ipHYDROGEN] ,
			/* ratio of H2_populations of 1 to 0 only */
			pop_ratio10 ,
			/* sum of all ortho and para */
			h2.ortho_density / SDIV(h2.para_density),
			T10,T20,T30,T31,T40,
			phycon.te ,
			hyperfine.Tspin21cm,T10_sum,T20_sum,T30_sum,T31_sum,T40_sum  );
	}
	else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) )
	{
		/* punch H2 lines - output the full emission-line spectrum */
		double thresh ;
		double renorm;
		if( h2.lgH2ON )
		{
			/* get the normalization line */
			if( LineSv[LineSave.ipNormWavL].sumlin > SMALLFLOAT )
			{
				renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].sumlin;
			}
			else
			{
				renorm = 1.;
			}
			if( renorm > SMALLFLOAT )
			{
				/* this is threshold for faintest line, normally 0, set with 
				* number on punch verner command */
				thresh = thresh_punline_h2/(float)renorm;
			}
			else
			{
				thresh = 0.f;
			}

			/* punch H2 line intensities at end of iteration 
			 * h2.nElecLevelOutput is elec level with 1 for ground, so this loop is < h2.nElecLevelOutput */
			for( iElecHi=0; iElecHi < h2.nElecLevelOutput; ++iElecHi )
			{
				for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
				{
					for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
					{
						/* now the lower levels */
						/* NB - X is the only lower level considered here, since we are only 
						* concerned with excited electronic levels as a photodissociation process
						* code exists to relax this assumption - simply change following to iElecHi */
						long int lim_elec_lo = 0;
						for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
						{
							long int nv = nVib_hi[iElecLo];
							if( iElecLo==iElecHi )
								nv = iVibHi;
							for( iVibLo=0; iVibLo<=nv; ++iVibLo )
							{
								long nr = nRot_hi[iElecLo][iVibLo];
								if( iElecLo==iElecHi && iVibHi==iVibLo )
									nr = iRotHi-1;
								for( iRotLo=Jlowest[iElecLo]; iRotLo<nr; ++iRotLo )
								{
									if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > thresh )
									{
										/* air wavelength in microns */
										/* WLAng contains correction for index of refraction of air */
										double wl = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng/1e4;
										/*ASSERT( abs(iRotHi-iRotLo)<=2 );*/

										fprintf(io, "%li-%li %c(%li)", 
											iVibHi , 
											iVibLo ,
											chMolBranch( iRotHi , iRotLo ) ,
											iRotLo );
										fprintf( io, "\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld", 
											iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo);
										/* WLAng contains correction for index of refraction of air */
										fprintf( io, "\t%.7f\t", wl );
										/*prt_wl print floating wavelength in Angstroms, in output format */
										prt_wl( io , H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng );
										fprintf( io, "\t%.3f\t%.3e", 
											log10(MAX2(1e-37,H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo])) + radius.Conv2PrtInten, 
											H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]*renorm );
										fprintf( io, "\t%.3f", energy_wn[iElecHi][iVibHi][iRotHi]*T1CM );
										fprintf( io, "\n");
									}
								}
							}
						}
					}
				}
			}
		}
	}
	else if( (strcmp(chJOB , "H2sp" ) == 0)  )
	{
		iVib = 0;
		iRot = 0;
#		if 0
		/* punch h2 special */
		fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			radius.depth_mid_zone , 
			H2_populations[0][iVib][iRot] ,
			radius.depth_mid_zone * H2_populations[0][iVib][iRot] ,
			H2_rad_rate_out[0][iVib][iRot] ,
			H2_rad_rate_in[iVib][iRot] ,
			H2_col_rate_out_ old[iVib][iRot] ,
			H2_col_rate_in_ old[iVib][iRot] );
#		endif
		fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			radius.depth_mid_zone , 
			H2_populations[0][iVib][iRot] ,
			H2Lines[1][1][1][0][iVib][iRot].pump,
			H2Lines[1][1][1][0][iVib][iRot].TauIn,
			H2Lines[1][1][1][0][iVib][iRot].TauCon);;
	}

	return;
}
 /*cdH2_Line determines intensity and luminosity of and H2 line.  The first
  * six arguments give the upper and lower quantum designation of the levels.
  * The function returns 1 if we found the line, 
  * and FALSE==0 if we did not find the line because ohoto-para transition
  * or upper level has lower energy than lower level  */
long int cdH2_Line(
	  /* indices for the upper level */
	  long int iElecHi, 
	  long int iVibHi ,
	  long int iRotHi ,
	  /* indices for lower level */
	  long int iElecLo, 
	  long int iVibLo ,
	  long int iRotLo ,
	  /* linear intensity relative to normalization line*/
	  double *relint, 
	  /* log of luminosity or intensity of line */
	  double *absint )
{

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

	/* these will be return values if we can't find the line */
	*relint = 0.;
	*absint = 0.;

	/* for now both electronic levels must be zero */
	if( iElecHi!=0 || iElecLo!=0 )
	{
#		ifdef DEBUG_FUN
		fputs( " <->cdH2_Line()\n", debug_fp );
#		endif
		return 0;
	}

	/* check that energy of first level is higher than energy of second level */
	if( energy_wn[iElecHi][iVibHi][iRotHi] < energy_wn[iElecLo][iVibHi][iRotHi] )
	{
#		ifdef DEBUG_FUN
		fputs( " <->cdH2_Line()\n", debug_fp );
#		endif
		return 0;
	}

	/* check that ortho-para does not change */
	if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] - H2_lgOrtho[iElecLo][iVibLo][iRotLo] != 0 )
	{
#		ifdef DEBUG_FUN
		fputs( " <->cdH2_Line()\n", debug_fp );
#		endif
		return 0;
	}

	/* exit if lines does not exist */
	if( !lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
	{
#		ifdef DEBUG_FUN
		fputs( " <->cdH2_Line()\n", debug_fp );
#		endif
		return 0;
	}

	/* does the normalization line have a positive intensity*/
	if( LineSv[LineSave.ipNormWavL].sumlin > 0. )
	{
		*relint = H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]/
			LineSv[LineSave.ipNormWavL].sumlin * LineSave.ScaleNormLine;
	}
	else
	{
		*relint = 0.;
	}

	/* return log of line intensity if it is positive */
	if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > 0. )
	{
		*absint = log10(H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]) + 
			radius.Conv2PrtInten;
	}
	else
	{
		/* line intensity is actually zero, return small number */
		*absint = -37.;
	}

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

	/* this indicates sucess */
	return 1;
}

/*lint +e662 possible access of out-of-bounds pointer */

