/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*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_Read_Cosmicray_distribution read distribution function for H2 population following cosmic ray collisional excitation */
/*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 */
#include "cddefines.h" 
#include "physconst.h" 
#include "punch.h" 
#include "hmi.h"
#include "path.h"
#include "prt.h"
#include "mole.h"
#include "secondaries.h"
#include "grainvar.h"
#include "phycon.h"
#include "rfield.h"
#include "hyperfine.h"
#include "thermal.h"
#include "lines.h"
#include "dense.h"
#include "radius.h"
#include "colden.h"
#include "taulines.h"
#include "h2_priv.h"
#include "h2.h"

/* 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 */
char chlgPara[2]={'P','O'};

/*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("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("LINE",chCard) )
	{
		/* punch H2 lines - all in X */
		strcpy( punch.chPunch[punch.npunch], "H2ln" );
		fprintf( punch.ipPnunit[punch.npunch], "#H2 line\tVhi\tJhi\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);
		}
	}
	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);

		/* 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\tpops\n" );
		}
		else 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\tv=0 rel pops\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) )
	{
		/* creation and destruction rates */
		strcpy( punch.chPunch[punch.npunch], "H2ra" );
		fprintf( punch.ipPnunit[punch.npunch], 
			"#depth\tA_V(star)\tn(Eval)\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tBigH2\tG(TH85)\tG(DB96)\tCR\tN(H2)\tEleclife \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\teHi\tvHi\tJHi\tvLo\tJLo\tfrac...\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) \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(hmi.H2_total)));

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

	fprintf( ioMEAN, " para");
	fprintf(ioMEAN,"%7.3f", log10(SDIV(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"  
	};
	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;

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

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


/*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 );
		
		/*lint -e662 possible creation of out of bounds pointer */
		/* 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,r);
			}
		}
		/*lint +e662 possible creation of out of bounds pointer */

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

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

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

#	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[] , char chTime[] , long int ipPun )
{
	long int iVibHi , iElecHi , iRotHi , iVibLo , iElecLo , iRotLo;
	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 last zone */
	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 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 */
							H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi] ) ,
							/* fraction of exits that were collisional */
							H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[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[iVibHi][iRotHi] ,
							/* radiation out */
							H2_col_rate_in[iVibHi][iRotHi],
							/* radiation in */
							H2_rad_rate_in[iVibHi][iRotHi]
							
							) ;
					}
				}
			}
		}
	}
	/* PUNCH H2 POPULATIONS 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 )
		{
			iElecHi = 0;
			iVibHi = 0;
			/* 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];
			}
			fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth , 
				h2.ortho_density , h2.para_density);
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
			{
				fprintf(io,"\t%.3e", 
					H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total );
			}
			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 ,
			ortho_colden );
		fprintf(io,"%i\t%i\t%.3e\tpara\n", 
			101 , 
			101 ,
			para_colden );
		/* total H2 column density */
		fprintf(io,"%i\t%i\t%.3e\ttotal\n", 
			0 , 
			0 ,
			colden.colden[ipCOLH2] );

		/* 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 ,
			/* 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 , 
			/* the Lyman Werner band dissociation, Bertoldi & Draine */
			hmi.H2_Solomon_dissoc_rate_BD96,
			/* the Lyman Werner band dissociation, big H2 mole */
			hmi.H2_Solomon_dissoc_rate_BigH2_total);
	}
	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 ,
			/* total cooling, equal to total heating */
			thermal.ctot , 
			/* H2 destruction by Solomon process, TH85 rate */
			hmi.H2_Solomon_dissoc_rate_TH85,
			/* H2 destruction by Solomon process, big H2 model rate */
			hmi.H2_Solomon_dissoc_rate_BigH2_total,
			/* 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 , "H2he" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
		/* PUNCH H2 HEATING - show heating dur to big H2 and TH85 */
		fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
			/* depth in cm */
			radius.depth ,
			/* total heating */
			thermal.htot , 
			/* heating by big H2 model */
			hmi.HeatH2Dexc_BigH2,
			/* heating by TH85 H2 model */
			hmi.HeatH2Dexc_TH85 ,
			hmi.HeatH2Dish_BigH2,
			hmi.HeatH2Dish_TH85 
			);

	}
	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%li\t%.2e\t%.3e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
			/* depth in cm */
			radius.depth ,
			/* 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,
			/* 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,
			/* H2 destruction by Solomon process, Bertoldi & Draine rate */
			hmi.H2_Solomon_dissoc_rate_BD96,
			/* H2 destruction by Solomon process, big H2 model rate */
			hmi.H2_Solomon_dissoc_rate_BigH2_total,
			/* 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,
			/* the column density in H2 */
			colden.colden[ipCOLH2],
			sumlife/SDIV( sumpop ) );
	}
	/* punch h2 solomon */
	else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) )
	{
#		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];
		/* solomon process, and where it came from */
		fprintf(io,"%.5e\t%.3e", 
			/* depth in cm */
			radius.depth ,
			/* H2 destruction by Solomon process, big H2 model rate */
			hmi.H2_Solomon_dissoc_rate_BigH2_total);
		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 xx of tot */
#			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;
										++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);
			for( i=0; i<nsave; ++i )
			{
				long int ip = ipOrdered[i];
				/*lint -e644 not init */
				fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f", 
					iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[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) )
	{
		double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40;
		double T10,T20,T30,T31,T40;
		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]);
			/* the corresponding temperature */
			T10 = -170.5/log(pop_ratio10 * 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(pop_ratio20 * 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(pop_ratio30 * 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(pop_ratio31 * 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(pop_ratio40 * H2_stat[0][0][0]/H2_stat[0][0][4]);
		}
		else
		{
			pop_ratio10 = 0.;
			T10 = 0.;
			T20 = 0.;
			T30 = 0.;
			T31 = 0.;
			T40 = 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\n" ,
			/* depth in cm */
			radius.depth ,
			/* 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 );
	}
	else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) )
	{
		/* the full emission line spectrum */
		double thresh ;
		double renorm;
		if( h2.lgH2ON )
		{
			/* get the normalization line */
			if( LineSv[LineSave.ipNormWavL].sumlin > 0. )
			{
				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;
			}

			iElecLo = 0;
			iElecHi = 0;
			/* punch lines at end of iteration */
			for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
			{
				for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
				{
					for( iVibLo=0; iVibLo<=iVibHi; ++iVibLo )
					{
						long nr = nRot_hi[iElecLo][iVibLo];
						if( iVibHi==iVibLo )
							/* max because cannot malloc 0 bytes */
							nr = MAX2(1,iRotHi-1);
						for( iRotLo=0; iRotLo<nr; ++iRotLo )
						{
							if( H2_SaveLine[iVibHi][iRotHi][iVibLo][iRotLo] > thresh )
							{
								/* air wavelength in microns */
								/* WLAng contains correction for index of refraction of air */
								double wl = H2Lines[0][iVibHi][iRotHi][0][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", 
									iVibHi , iRotHi , 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[0][iVibHi][iRotHi][0][iVibLo][iRotLo].WLAng );
								fprintf( io, "\t%.3f\t%.3e", 
									log10(MAX2(1e-37,H2_SaveLine[iVibHi][iRotHi][iVibLo][iRotLo])) + radius.Conv2PrtInten, 
									H2_SaveLine[iVibHi][iRotHi][iVibLo][iRotLo]*renorm );
								fprintf( io, "\t%.3f", energy_wn[0][iVibHi][iRotHi]*T1CM );
								fprintf( io, "\n");
							}
						}
					}
				}
			}
		}
	}
	else if( (strcmp(chJOB , "H2sp" ) == 0)  )
	{
		int iVib = 0;
		int 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 , 
			H2_populations[0][iVib][iRot] ,
			radius.depth * H2_populations[0][iVib][iRot] ,
			H2_rad_rate_out[iVib][iRot] ,
			H2_rad_rate_in[iVib][iRot] ,
			H2_col_rate_out[iVib][iRot] ,
			H2_col_rate_in[iVib][iRot] );
#		endif
		fprintf(io,"%.4e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			radius.depth , 
			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;
}

