/* 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_Create create H2 molecules, called by ContCreatePointers after continuum mesh is set up */
/*H2_ReadEnergies read energies for all electronic levels */
/*H2_ReadTransprob read transition probabilities */
/*H2_ReadCollRates read transition probabilities */
/*H2_ContPoint set the ipCont struc element for the H2 molecule, called by ContCreatePointers */
/*H2_Accel radiative acceleration due to H2 */
/*H2_RadPress rad pre due to h2 lines called in PresTotCurrent */
/*H2_Zero zero out vars in the large H2 molecule, called from zero */
/*H2_InterEnergy internal energy of H2 called in PresTotCurrent */
/*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from PunchLineStuff */
/*H2_RT_diffuse do emission from H2 - called from RT_diffuse */
/*H2_itrzn - average number of H2 pop evaluations per zone */
/*H2_RTMake do RT for H2 - called from RT_line_all */
/*H2_collid_rates - set H2 collision rates */
/*H2_RT_tau_inc increment optical depth for the H2 molecule, called from RT_tau_inc */
/*H2_LineZero initialize optical depths in H2, called from RT_tau_init */
/*H2_ParsePunch parse the punch h2 command */
/*H2_RT_tau_reset the large H2 molecule, called from RT_tau_reset */
/*H2_Reset called to reset variables that are needed after an iteration */
/*H2_Colden maintain H2 column densities within X */
/*H2_LevelPops do level H2_populations for H2, called by Hydrogenic */
/*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
/*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
/*H2_PunchDo punch some properties of the large H2 molecule */
/*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
/*H2_Prt_column_density print H2 column density, called from prtcolumns */
/*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
/*H2_Level_low_matrix evaluate CO rotation cooling */
/*H2_cooling evaluate cooling and heating due to H2 molecule */
/*H2_vib_dist evaluates the vibration distribution for H2 formed on grains */
/*cdH2_colden return column density in H2, negative -1 if cannot find state,
 * header is cddrive */
/*H2_DR choose next zone thickness based on H2 big molecule */
/*TODO	2	put in excited molecular dissociation from v >=4 as in hazy 2 */
/* turn this flag on to do minimal debug print of pops */
#define	PRT_POPS	FALSE
/* this is limit to number of loops over H2 pops before giving up */
#define	LIM_H2_POP_LOOP	100
#include "cddefines.h" 
#include "physconst.h" 
#include "taulines.h" 
#include "lines.h" 
#include "atoms.h" 
#include "converge.h" 
#include "secondaries.h" 
#include "trace.h" 
#include "hmi.h" 
#include "rt.h" 
#include "radius.h" 
#include "ipoint.h" 
#include "phycon.h" 
#include "thermal.h" 
#include "dense.h" 
#include "rfield.h" 
#include "lines_service.h" 
#include "mole.h"
#include "h2_priv.h"
#include "h2.h"

/* this is the number of electronic levels - in h2.h */
/* #define N_H2_ELEC	7*/

/* the number of times the H2 molecules has been called in this iteration.  For the
 * very first call we will use lte for the level H2_populations, for later calls
 * use the last solution */
static long int nCallH2_this_iteration;

/* this counts how many times we go through the H2 level H2_populations loop */
static long int loop_h2_pops;

float H2_te_hminus[nTE_HMINUS] = {10.,30.,100.,300.,1000.,3000.,10000.};

/* the total population in each elec state */
static float pops_per_elec[N_H2_ELEC];

/* this will contain a vector for collisions within the X ground elec state,
 * CollRateFit[coll_type][vib_up][rot_up][vib_lo][rot_lo][3] */
static float collider_density[N_X_COLLIDER];
static float collider_density_total;

/* the order of the electronic states is
 * X, B, C+, C-, B', D+, and D- */
/* this will be the number of vibration levels within each elec */
/* number of vib states within electronic states from
 * >>refer	H2	energies	Abgrall, */
long int nVib_hi[N_H2_ELEC] = {14 , 37 , 13 , 13, 9, 2 , 2};

/* this integer is added to rotation quantum number J for the test of whether
 * a particular J state is ortho or para - the state is ortho if J+below is odd,
 * and para if J+below is even */
int H2_nRot_add_ortho_para[N_H2_ELEC] = {0 , 1 , 1 , 0, 1, 1 , 0};

/* this gives the first rotational state for each electronic state - J=0 does
 * not exist when Lambda = 1 */
long int Jlowest[N_H2_ELEC] = {0 , 0 , 1  , 1 , 0 , 1 , 1 };

#if 0
static float SolomonRateMax;
static long nRot_SolomonRateMax , nVib_SolomonRateMax ,
	nRotHi_SolomonRateMax , nVibHi_SolomonRateMax , nEkcHi_SolomonRateMax;
#endif

/* dissociation energies (cm-1) for each electronic state, from 
 * >>refer	H2	energies	Sharp, T. E., 1971, Atomic Data, 2, 119 */
/* >>chng 02 oct 08, improved energies */
double H2_DissocEnergies[N_H2_ELEC] = 
{ 36118.11, 118375.6, 118375.6, 118375.6, 118375.6,133608.6,133608.6 };
/* original values
 { 36113., 118372., 118372., 118372., 118372.,0.,0. };*/

/* number of rotation levels within each elec - vib */
/*lint -e785 too few init for aggregate */
long int nRot_hi[N_H2_ELEC][50]=
	/* ground, X */
	{ {31, 30, 28, 27, 25, 
		23, 22, 20, 18, 16, 
		14, 12, 10,  7,  3 } ,
	/* B */
	{25,25,25,25,25,25,25,25, 25,25,
	 25,25,25,25,25,25,25,25, 25,25,
	 25,25,25,25,25,25,25,25, 23,21,
	 19,17,15,15,11,9,7, 7},
	/* C plus */
	{ 25, 25, 25, 25, 24, 23, 21, 19, 17, 14, 12, 10, 6, 2 },
	/* C minus (the same) */
	{ 25, 25, 25, 25, 24, 23, 21, 19, 17, 15, 13, 10, 7, 2 },
	/* B primed */
	{19,17, 14, 12, 9, 8, 7, 7, 4, 1 },
	/* D plus */
	{13, 10, 5},
	/* D minus */
	{25 , 25 ,25 } }
	;
/*lint +e785 too few init for aggregate */

static long int nH2_pops , nH2_zone;
/*H2_itrzn - average number of H2 pop evaluations per zone */
double H2_itrzn( void )
{
	if( h2.lgH2ON && nH2_zone>0 )
	{
		return( (double)nH2_pops / (double)nH2_zone );
	}
	else
	{
		return 0.;
	}
}

/*H2_Reset called by IterRestart to reset variables that are needed after an iteration */
void H2_Reset( void )
{

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

	if(mole.lgH2_TRACE) 
		fprintf(ioQQQ,
		"\n***************H2_Reset called, resetting nCallH2_this_iteration, zone %.2f iteration %li\n", 
		fnzone,
		iteration );

	/* number of times large molecules evaluated in this iteration,
	 * is FALSE if never evaluated, on next evaluation will start with lte populations */
	nCallH2_this_iteration = 0;

	/* these remember the largest and smallest factors needed to
	 * renormalize the H2 chemistry */
	h2.renorm_max = 1.;
	h2.renorm_min = 1.;

	/* counters used by H2_itrzn to find number of calls of h2 per zone */
	nH2_pops  = 0;
	nH2_zone = 0;

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

	return;

}

/*H2_init - called to initialize things from cdInit */
void H2_Init(void)
{

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

	/* the number of electronic quantum states to include.
	 * To do both Lyman and Werner bands want nelec = 3,
	 * default is to do all bands included */
	mole.n_h2_elec_states = N_H2_ELEC;
	h2.nCallH2_this_zone = 0;

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

	return;

}


/*H2_Zero zero out vars in the large H2 molecule, called from zero 
 * before any commands are parsed 
 * NB - this routine is called before space allocated - must not zero out
 * any allocated arrays */
void H2_Zero( void )
{

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

	/* this is the smallest ratio of H2/H where we will bother with the large H2 molecule
	 * this value was chosen so that large mole used at very start of TH85 standard pdr,
	 * NB - this appears in headinfo and must be updated there if changed here */
	/* >>chng 03 jun 02, from 1e-6 to 1e-8 - in orion veil large H2 turned on half way
	 * across, and solomon process was very fast since all lines optically thin.  correct
	 * result has some shielding, so reset to lower value so that H2 comes on sooner. */
	mole.H2_to_H_limit = 1e-8;

	h2.lgH2ON = FALSE;
	mole.lgH2_TRACE = FALSE;

	/* counters used by H2_itrzn to find number of calls of h2 per zone */
	nH2_pops  = 0;
	nH2_zone = 0;

	/* option to scramble collision data */
	mole.lgH2_NOISE = FALSE;
    mole.lgH2_NOISECOSMIC = FALSE; 

	/* option to turn off or on gbar collisions of the collision rate,
	 * default is to have it on */
	/* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
	mole.lgColl_gbar = TRUE;

	/* include collision rates that come from real calculations,
	 * off with atom h2 collisions off command */
	mole.lgColl_deexec_Calc = TRUE;
	mole.lgColl_dissoc_coll = TRUE;

	/* these remember the largest and smallest factors needed to
	 * renormalize the H2 chemistry */
	h2.renorm_max = 1.;
	h2.renorm_min = 1.;

	nCallH2_this_iteration = 0;
	h2.ortho_density = 0.;
	h2.para_density = 0.;

	hmi.H2_Solomon_dissoc_rate_BigH2_total = 0.;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
	hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
	hmi.H2_photodissoc_BigH2 = 0.;

	/* say that H2 has never been computed */
	hmi.lgBigH2_evaluated = FALSE;

	hmi.lgH2_Thermal_BigH2 = TRUE;
	hmi.lgH2_Chemistry_BigH2 = TRUE;

	hmi.H2_star_BigH2 = 0.;

	if( !lgH2_READ_DATA )
	{
		/* the number of electronic levels in the H2 molecule,
		 * to just do the Lyman and Werner bands set to 3 -
		 * reset with atom h2 levels command,
		 * default is all levels with data */
		mole.n_h2_elec_states = N_H2_ELEC;
	}

	/* the number of levels used in the matrix solution
	 * of the level H2_populations - set with atom h2 matrix nlevel,
	 * >>chng 04 oct 05, make default 30 levels */
	nXLevelsMatrix = 30;

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

	return;

}

/* set the ipCont struc element for the H2 molecule, called by ContCreatePointers */
void H2_ContPoint( void )
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;

	if( !h2.lgH2ON )
		return;

	/* set array index for line energy within continuum array */
	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, change to > 0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont = 
									(int)ipLineEnergy(
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD , 
									"H2  " , 0 );
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipFine = 
									ipFineCont(
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD );
							}
						}
					}
				}
			}
		}
	}
	return;
}

/* ===================================================================== */
/* radiative acceleration due to H2 called in rt_line_driving */
double H2_Accel(void)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;
	double h2_drive;

	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return(0.);

	/* this routine computes the line driven radiative acceleration
	 * due to H2 molecule*/

	h2_drive = 0.;
	/* 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. )
							{
								h2_drive += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump*
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg*
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc;
							}
						}
					}
				}
			}
		}
	}
	return h2_drive;
}

/* ===================================================================== */
/* rad pre due to h2 lines called in PresTotCurrent */
double H2_RadPress(void)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;
	double press;

	/* will be used to check on size of opacity, was capped at this value */
	float smallfloat=SMALLFLOAT*10.f;

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

	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return(0.);

	press = 0.;
	/* 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. )
							{
								if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi > smallfloat &&
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc > smallfloat )
								{
									double RadPres1 = 5.551e-2*(
										powi(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN/1.e6,4))*
										(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi/
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi)/
										(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo/
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo)*
										RT_LineWidth(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]);

										press += RadPres1;
								}
							}
						}
					}
				}
			}
		}
	}

	if(mole.lgH2_TRACE) 
		fprintf(ioQQQ,
		"  H2_RadPress returns, radiation pressure is %.2e\n", 
		press );

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

	return press;
}

/* ===================================================================== */
/* internal energy of H2 called in PresTotCurrent */
double H2_InterEnergy(void)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;
	double energy;

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

	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return(0.);

	energy = 0.;
	/* 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. )
							{
								energy += 
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi* 
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg;
							}
						}
					}
				}
			}
		}
	}

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

	return energy;
}

/*H2_RT_diffuse do emission from H2 - called from RT_diffuse */
void H2_RT_diffuse(void)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;

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

	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return;

	/* loop over all possible lines */
	/* NB - this loop does not include the electronic lines */
	for( iElecHi=0; iElecHi<1; ++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. )
							{
								outline( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
#								if 0
								/* pointer to line energy in continuum array */
								ip = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont;
								rfield.outlin[ip-1] += (float)(
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots*
									VolFac*opac.tmn[ip-1]*
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot);
								rfield.reflin[ip-1] += (float)(
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots * ref);
#								endif
							}
						}
					}
				}
			}
		}
	}

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

	return;
}

/* do RT for H2 - called from RT_line_all */
void H2_RTMake( int lgDoEsc , int lgUpdateFineOpac )
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo ;

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

	if( !h2.lgH2ON )
		return;
	/* this routine drives calls to make RT relations for H2 molecule */
	/* 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, change test from !=0 to >0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								/* >>chng 03 jun 18, added 4th parameter in call to this routine - says to not
								 * include self-shielding of line across this zone.  This introduces a dr dependent
								 * variation in the line pumping rate, which made H2 abundance fluctuate due to
								 * Solomon process having slight dr-caused mole. */
								RT_line_one( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , lgDoEsc , lgUpdateFineOpac,FALSE);
							}
						}
					}
				}
			}
		}
	}
	/* this is to take mean of upward transition probability */
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			double sumpop = 0.;
			double sumpopA = 0.;
			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 )
					{
						/* 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 */
						iElecLo = 0;
						/* want to include all vib states in lower level if different elec level,
						* but only lower vib levels if same elec level */
						for( iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
						{
							long nr = nRot_hi[iElecLo][iVibLo];
							for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
							{
								/* >>chng 03 feb 14, change test from !=0 to >0 */
								if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
								{
									sumpop += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo;
									sumpopA += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo*
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul;
								}
							}
						}
					}
				}
			}
			fprintf(ioQQQ,"sumpop = %.3e sumpopA= %.3e A=%.3e\n", sumpop, sumpopA, 
				sumpopA/SDIV(sumpop) );
		}
	}
	return;
}

/* increment optical depth for the H2 molecule, called from RT_tau_inc which is called  by cloudy,
 * one time per zone */
void H2_RT_tau_inc(void)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;

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

	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return;

	/* remember largest and smallest chemistry renorm factor -
	 * if both networks are parallel will be unity,
	 * but only do this after we have stable solution */
	if( nzone > 0 && nCallH2_this_iteration>2 )
	{
		h2.renorm_max = MAX2( H2_renorm_chemistry , h2.renorm_max );
		h2.renorm_min = MIN2( H2_renorm_chemistry , h2.renorm_min );
	}
	/*fprintf(ioQQQ,"DEBUG popopc\t%.3e",H2Lines[1][0][0][0][0][1].TauCon);*/

	/* 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. )
							{
								/*double save = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc;
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc = 0.;bro ken ( );*/
								RT_line_one_tauinc( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ,
									-9 , iRotHi , iVibLo , iRotLo);
								/*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc = save;*/
							}
						}/* iRotLo loop */
					}/* iVibLo loop */
				}/* iElecLo loop */
			}/* iRotHi loop */
		}/* iVibHi loop */
	}/* iElecHi loop */
	/*fprintf(ioQQQ,"\t%.3e\n",H2Lines[1][0][0][0][0][1].TauCon);*/

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

	return;
}


/* initialize optical depths in H2, called from RT_tau_init */
void H2_LineZero( void )
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;

	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 )
						{
							/*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. )
							{
								/* >>chng 03 feb 14, use EmLineZero rather than explicit sets */
								EmLineZero( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
							}
						}
					}
				}
			}
		}
	}
	return;
}

/* the large H2 molecule, called from RT_tau_reset */
void H2_RT_tau_reset( void )
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;

	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, change test from !=0 to >0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								/* inward optical depth */
								RT_line_one_tau_reset( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , 0.5);
							}
						}
					}
				}
			}
		}
	}
}

/* this is fraction of population that is within levels done with matrix */
static double frac_matrix;

/*H2_Level_low_matrix evaluate lower populations within X */
static void H2_Level_low_matrix(
	/* total abundance within matrix */
	float abundance )
{

	/* will need to MALLOC space for these but only on first call */
	static double **data, 
	  **dest, 
	  /* pump[low][high] is rate (s^-1) from lower to upper level */
	  **pump,
	  **CollRate_levn,
	  *pops,
	  *create,
	  *destroy,
	  *depart,
	  /* statistical weight */
	  *stat_levn ,
	  /* excitation energies in kelvin */
	  *excit;
	long int ip, ip1;
	static int lgFirst=TRUE;
	long int i,
		j,
		ilo , 
		ihi,
		iElec,
		iElecHi,
		iVib,
		iRot,
		iVibHi,
		iRotHi;
	int lgDeBug,lgNegPop,
		lgOrthoPara_possible;
	double rot_cooling , dCoolDT ;
	double pop_e, factor, sum_pops;
	static long int ndimMalloced = 0;
	double rateout , ratein , energy_ilo;

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

	/* option to not use the matrix */
	if( nXLevelsMatrix <= 1 )
		return;

	if( lgFirst )
	{
		/* check that not more levels than there are in X */
		if( nXLevelsMatrix > nLevels_per_elec[0] )
		{
			/* number is greater than number of levels within X */
			fprintf( ioQQQ, 
				" The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n",
				nLevels_per_elec[0]);
			puts( "[Stop in ParseAtomH2]" );
			cdEXIT(EXIT_FAILURE);
		}
		/* will never do this again */
		lgFirst = FALSE;
		/* remember how much space we malloced in case ever called with more needed */
		ndimMalloced = nXLevelsMatrix;
		/* allocate the 1D arrays*/
		if( (excit = (double *)MALLOC( sizeof(double)*(size_t)(nXLevelsMatrix) )) == NULL )
			BadMalloc();
		if( (stat_levn = (double *)MALLOC( sizeof(double)*(size_t)(nXLevelsMatrix) )) == NULL )
			BadMalloc();
		if( (pops = (double *)MALLOC( sizeof(double)*(size_t)(nXLevelsMatrix) )) == NULL )
			BadMalloc();
		if( (create = (double *)MALLOC( sizeof(double)*(size_t)(nXLevelsMatrix) )) == NULL )
			BadMalloc();
		if( (destroy = (double *)MALLOC( sizeof(double)*(size_t)(nXLevelsMatrix) )) == NULL )
			BadMalloc();
		if( (depart = (double *)MALLOC( sizeof(double)*(size_t)(nXLevelsMatrix) )) == NULL )
			BadMalloc();
		/* create space for the 2D arrays */
		if( (pump = ((double **)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (CollRate_levn = ((double **)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (dest = ((double **)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double *)))) == NULL )
			BadMalloc();
		if( (data = ((double **)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double *)))) == NULL )
			BadMalloc();
		for( i=0; i<(nXLevelsMatrix); ++i )
		{
			if( (pump[i] = ((double *)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (CollRate_levn[i] = ((double *)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (dest[i] = ((double *)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double )))) == NULL )
				BadMalloc();
			if( (data[i] = ((double *)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double )))) == NULL )
				BadMalloc();
			/*if( (ipdest[i] = ((long int *)MALLOC((size_t)(nXLevelsMatrix)*sizeof(double )))) == NULL )
				BadMalloc();*/
		}

		for( j=0; j < nXLevelsMatrix; j++ )
		{
			stat_levn[j]=0;
			excit[j] =0;
		}
		/* the statistical weights of the levels
		 * and excitation potentials of each level relative to ground */
		for( j=0; j < nXLevelsMatrix; j++ )
		{
			/* obtain the proper indices for the upper level */
			ip = H2_ipX_ener_sort[j];
			iVib = ipVib_H2_energy_sort[ip];
			iRot = ipRot_H2_energy_sort[ip];

			/* statistical weights for each level */
			stat_levn[j] = H2_stat[0][iVib][iRot];
			/* excitation energy of each level relative to ground, in K */
			excit[j] = energy_wn[0][iVib][iRot]*T1CM;
			/*fprintf(ioQQQ,"DEBUG H2energy\t%li\t%li\t%.4e\n",
				iVib , 
				iRot ,
				excit[j] );*/
		}

		for( j=0; j < nXLevelsMatrix-1; j++ )
		{
			/*fprintf(ioQQQ,"DEBUG H2ener2 %li %.4e %.4e\n",
				j, excit[j+1] , excit[j] );*/
			/* make sure that the energies are ok */
			ASSERT( excit[j+1] > excit[j] );
		}

	}
	/* end mallocing space and creating constant terms */

	/* this is test for call with too many rotation levels to handle - logic needs
	 * for largest model atom to be called first */
	if( nXLevelsMatrix > ndimMalloced )
	{
		fprintf(ioQQQ," H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n");
		puts( "[Stop in H2_Level_low_matrix]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* when no o-p collisions are possible some downward collisions rates are also zero,
	 * test against 1e-10 since this density multiplies rate coef and must fit into a float 
	 * all this does is defeat an assert - no physics is affected */
	lgOrthoPara_possible = TRUE;
	if( collider_density[0] + collider_density[4] <1e-10 )
		lgOrthoPara_possible = FALSE;

	/* all elements are used, and must be set to zero */
	for( i=0; i < nXLevelsMatrix; i++ )
	{
		create[i] = 0.;
		destroy[i] = 0.;
		pops[i] = 0.;
		depart[i] = 0;
		for( j=0; j < nXLevelsMatrix; j++ )
		{
			data[j][i] = 0.;
			dest[j][i] = 0.;
			pump[j][i] = 0.;
			CollRate_levn[j][i] = 0.;
		}
	}

	/* this is the total population in all excited electronic states */
	pop_e = 0;
	for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
	{
		pop_e += pops_per_elec[iElecHi];
	}

	/* find all radiative interactions within matrix, and between
	 * matrix and upper X and excited electronic states */
	iElec = 0;
	for( ilo=0; ilo < nXLevelsMatrix; ilo++ )
	{
		ip = H2_ipX_ener_sort[ilo];
		iRot = ipRot_H2_energy_sort[ip];
		iVib = ipVib_H2_energy_sort[ip];
		rateout = 0.;
		ratein = 0.;
		/* this loop does radiative decays from upper states inside matrix, 
		 * and upward pumps within matrix region into this lower level */
		for( ihi=ilo+1; ihi<nXLevelsMatrix; ++ihi )
		{
			ip = H2_ipX_ener_sort[ihi];
			iRotHi = ipRot_H2_energy_sort[ip];
			iVibHi = ipVib_H2_energy_sort[ip];
			/* general case - but line may not actually exist */
			if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
			{
				if( H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Aul > 0. )
				{
					/* NB - the dest prob is included in the total and the dest is set to zero
					* since we want to only count one ots rate, in main calling routine,
					* and do not want matrix solver below to include it */
					data[ilo][ihi] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Aul*(
						H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Pesc + 
						H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Pdest +
						/* >>chng 04 sep 08, do not include pump here since done in atom leveln, gs */
						H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Pelec_esc)/*+H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].pump *
							H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].gLo/
							H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].gHi*/;
					dest[ilo][ihi] = 0.;
					pump[ilo][ihi] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].pump;
				}
			}
		}

		iElecHi = 0;
		iElec = 0;
		/* now do all levels within X, which are above nXLevelsMatrix,
		 * the highest level inside the matrix */
		for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
		{
			ip = H2_ipX_ener_sort[ihi];
			iRotHi = ipRot_H2_energy_sort[ip];
			iVibHi = ipVib_H2_energy_sort[ip];
			if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
			{
				if( H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul > 0. )
				{
					/* these will enter as net creation terms in creation vector, with
					 * units cm-3 s-1
					 * radiative transitions from above the matrix within X */
					ratein +=
						H2_populations[iElecHi][iVibHi][iRotHi] *
						(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul*
						(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pesc + 
						H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pelec_esc + 
						H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pdest)+H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump *
							H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gLo/
							H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gHi);
					/* rate out has units s-1 - destroys current lower level */
					rateout +=
						H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump;
				}
			}
		}

		/* Solomon process, this sum dos all pump and decays from all electronic excited states */
		/* pop_e is zero the first time this is called in an iteration, and might
		 * become small if there is no uv light - the lower level is [0][iVib][iRot] */
		if( pop_e > SMALLFLOAT )
		{
			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 )
					{
						/* the rate electrons enter this state from excited elec states */
						/* >>chng 03 feb 14, from !=0 to >0 */
						if( H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul > 0. )
						{
							/* these will enter as net creation terms in creation vector, with
							 * units cm-3 s-1 
							 * radiative transitions to excited elec states */
							ratein +=
								H2_populations[iElecHi][iVibHi][iRotHi] *
								(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul*
								(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pesc + 
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pelec_esc + 
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pdest)+H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump *
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gLo/
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gHi);
							/* rate out has units s-1 - destroys current lower level */
							rateout +=
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump + 
								/* >>chng 04 sep 08, include destryctio by secondaries gs */
								hmi.lgLeidenCRHack*secondaries.x12tot /1e3f;

							/* these are general entries and exits going into vector */
							/*if( iRot==5 )
								fprintf(ioQQQ," rad 1\t%li\t%li\t%e\t%e\n",iVibHi,iRotHi,
								H2_populations[iElecHi][iVibHi][iRotHi] *
								(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul*
								(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pesc + 
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pelec_esc + 
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pdest)+H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump *
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gLo/
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gHi),
								rateout );*/
						}
					}
				}
			}/* end loop over all electronic excited states */
		}/* end test for significant population in electronic excited states */
		create[ilo] += ratein;
		destroy[ilo] += rateout;

		/* rates H2 is created from grains and H- units cm-3 s-1, evaluated in mole_H2_form */
		create[ilo] += H2_X_formation[iVib][iRot];
		/* >>chng 04 sep 14, cosmic ray or secondary photodissociation into mainly H2+ ,GS */
		destroy[ilo] += secondaries.csupra[ipHYDROGEN][0]*2.02f * mole.lgColl_deexec_Calc;
		/* >>chng 04 sep 14, continuum photodissociation into 2H, GS */ 
		if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
		{ 
			destroy[ilo] += rfield.flux_accum[H2_ipPhoto[iVib][iRot]-1]*0.25e-18f;
		}
		/* >>chng 04 sep 14, collisional dissociation into continuum of X, GS */
		energy_ilo = H2_DissocEnergies[0] - energy_wn[0][iVib][iRot];
		ASSERT( energy_ilo > 0. );
		destroy[ilo] += collider_density_total* 
			1e-14f * (float)sexp(energy_ilo/phycon.te_wn) * mole.lgColl_dissoc_coll;

		/* >>chng 04 sep 15,excitation within X due to cosmic rays,
		 * csupra is total ionization rate, cr_rate is branching ratio, GS */

			if( iVib < CR_VIB && !hmi.lgLeidenCRHack)
			{
				/* cr_rate[CR_X][CR_VIB][CR_J][CR_EXIT];*/
				if( iRot> 1 && iRot < CR_J )
				{
					/* collision J to J-2 */
					destroy[ilo] += 
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][0] * mole.lgColl_deexec_Calc ;
					
					ipEnergySort[iVib][iRot-2]=ip;
                    ip1 = H2_ipX_ener_sort[ip];
			

					create[ip1] += H2_old_populations[0][iVib][iRot]*
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][0] * mole.lgColl_deexec_Calc ;
				}

				if( iRot<CR_J-2 && iRot+2 <= nRot_hi[iElec][iVib])
				{   
					ipEnergySort[iVib][iRot+2]=ip;
                    ip1 = H2_ipX_ener_sort[ip];

					destroy[ip1] += 
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][2] * mole.lgColl_deexec_Calc ;
					create[ip1] += H2_old_populations[0][iVib][iRot]*
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][2] * mole.lgColl_deexec_Calc ;
				}
			}
             /* >>chng 04 sep 15, ortho para conversion, and deexcitation, while on grain surface,
			 * H2 + grain => H2 + grain, do not want to add rate to 0,0 ; GS */
			if( iVib!= 0 || iRot!=0 )
			{
				/* this is the ground v J = 1->0 transition */
				if( iVib == 0 && iRot == 1 )
				{
					destroy[ilo] += 
						/* H2 grain interactions
						* H2 ortho - para conversion on grain surface,
						* rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
						(float)hmi.rate_h2_allX_2_J0_grains * mole.lgColl_deexec_Calc ;
				}
				/* add both rates if not 0,0 or 0,1 */
				else if( iVib >0  || iRot > 1 )
				{
					destroy[ilo] += 
						/* H2 grain interactions
						* H2 ortho - para conversion on grain surface,
						* rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
						(float)(hmi.rate_h2_allX_2_J0_grains + 
						/* rate (s-1) all v,J levels go to 0, regardless of nuclear spin */
						hmi.rate_h2_ortho_para_conserve) * mole.lgColl_deexec_Calc ;

				}
			}

	}

	/* now evaluate the collision rates */
	lgDeBug = FALSE;
	/* loop over lower levels within matrix */
	for( ilo=0; ilo < nXLevelsMatrix; ilo++ )
	{    
		ip = H2_ipX_ener_sort[ilo];
		iRot = ipRot_H2_energy_sort[ip];
		iVib = ipVib_H2_energy_sort[ip];
		if(lgDeBug)fprintf(ioQQQ,"DEBUG H2_Level_low_matrix, ilo=%li",ilo);
		for( ihi=ilo+1; ihi < nXLevelsMatrix; ihi++ )
		{
			int nColl;
			ip = H2_ipX_ener_sort[ihi];
			iRotHi = ipRot_H2_energy_sort[ip];
			iVibHi = ipVib_H2_energy_sort[ip];
			/* first do downward deexcitation rate */
			CollRate_levn[ihi][ilo] = 0.;
			/* >>chng 04 sep 14, do all levels */
			for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
			{
				CollRate_levn[ihi][ilo] +=
					H2_CollRate[nColl][iVibHi][iRotHi][iVib][iRot]*collider_density[nColl];
			}
			ASSERT( !lgOrthoPara_possible || !mole.lgColl_deexec_Calc || !mole.lgColl_gbar ||
				CollRate_levn[ihi][ilo]>SMALLFLOAT );

			/*create[ilo] +=CollRate_levn[ihi][ilo]*H2_populations[0][iVibHi][iRotHi];*/
			if(lgDeBug)fprintf(ioQQQ,"\t%.1e",CollRate_levn[ihi][ilo]);

			/* now get upward excitation rate */
			CollRate_levn[ilo][ihi] = CollRate_levn[ihi][ilo]*
				H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
				H2_stat[0][iVibHi][iRotHi] / 
				H2_stat[0][iVib][iRot];


			/*destroy[ilo] += CollRate_levn[ilo][ihi];*/
		}
		if(lgDeBug)fprintf(ioQQQ,"\n");
		/* now do all collisions for levels within X, which are above nXLevelsMatrix,
		 * the highest level inside the matrix */
		iElecHi = 0;
        
		for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
		{
			int nColl;
			ip = H2_ipX_ener_sort[ihi];
			iRotHi = ipRot_H2_energy_sort[ip];
			iVibHi = ipVib_H2_energy_sort[ip];
			ratein = 0;
			rateout = 0;
			/* first do downward deexcitation rate */
			/* >>chng 04 sep 14, do all levels */
			for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
			{
				ratein +=
					H2_CollRate[nColl][iVibHi][iRotHi][iVib][iRot]*collider_density[nColl];
			}
			if(lgDeBug)fprintf(ioQQQ,"\t%.1e",ratein);
			ASSERT( !lgOrthoPara_possible || !mole.lgColl_deexec_Calc || !mole.lgColl_gbar ||
				ratein>SMALLFLOAT );

			/* now get upward excitation rate */
			rateout = ratein *
				H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
				H2_stat[0][iVibHi][iRotHi]/H2_stat[0][iVib][iRot];
			/* >>chng 04 aug 19, use main stat weight vector */
				/*H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].gHi / 
				H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].gLo;*/
			/*if( iRot==5 )
			{
				fprintf(ioQQQ,"DEBUG he coll 1\t%li\t%li\t%e\t%e\n",iVibHi,iRotHi,
				ratein,H2_populations[iElecHi][iVibHi][iRotHi] );
			}*/

			/* these are general entries and exits going into vector */
			create[ilo] += ratein*H2_populations[iElecHi][iVibHi][iRotHi];
			destroy[ilo] += rateout;
		}
	}

	lgDeBug = FALSE;
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			fprintf(ioQQQ,"DEBUG H2 matexcit");
			for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
			{
				fprintf(ioQQQ,"\t%li",ilo );
			}
			fprintf(ioQQQ,"\n");
			for(ihi=0; ihi<nXLevelsMatrix;++ihi)
			{
				fprintf(ioQQQ,"\t%.2e",excit[ihi] );
			}
            fprintf(ioQQQ,"\n");
			for(ihi=0; ihi<nXLevelsMatrix;++ihi)
			{
				fprintf(ioQQQ,"\t%.2e",stat_levn[ihi] );
			}
			fprintf(ioQQQ,"\n");

			fprintf(ioQQQ,"data [n][]\\[][n] = Aul*Pesc\n");
			for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
			{
				fprintf(ioQQQ,"\t%li",ilo );
			}
			fprintf(ioQQQ,"\n");
			for(ihi=0; ihi<nXLevelsMatrix;++ihi)
			{
				fprintf(ioQQQ,"%li", ihi);
				for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
				{
					fprintf(ioQQQ,"\t%.2e",data[ihi][ilo] );
				}
				fprintf(ioQQQ,"\n");
			}

			fprintf(ioQQQ,"pump [n][]\\[][n]\n");
			for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
			{
				fprintf(ioQQQ,"\t%li",ilo );
			}
			fprintf(ioQQQ,"\n");
			for(ihi=0; ihi<nXLevelsMatrix;++ihi)
			{
				fprintf(ioQQQ,"%li", ihi);
				for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
				{
					fprintf(ioQQQ,"\t%.2e",pump[ihi][ilo] );
				}
				fprintf(ioQQQ,"\n");
			}

			fprintf(ioQQQ,"CollRate_levn [n][]\\[][n]\n");
			for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
			{
				fprintf(ioQQQ,"\t%li",ilo );
			}
			fprintf(ioQQQ,"\n");
			for(ihi=0; ihi<nXLevelsMatrix;++ihi)
			{
				fprintf(ioQQQ,"%li", ihi);
				for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
				{
					fprintf(ioQQQ,"\t%.2e",CollRate_levn[ihi][ilo] );
				}
				fprintf(ioQQQ,"\n");
			}
			fprintf(ioQQQ,"SOURCE");
			for(ihi=0; ihi<nXLevelsMatrix;++ihi)
			{
				fprintf(ioQQQ,"\t%.2e",create[ihi]);
			}
			fprintf(ioQQQ,"\nSINK");
			for(ihi=0; ihi<nXLevelsMatrix;++ihi)
			{
				fprintf(ioQQQ,"\t%.2e",destroy[ihi]);
			}
			fprintf(ioQQQ,"\n");
		}
	}

	atom_levelN(
		/* number of levels */
		nXLevelsMatrix,
		abundance,
		stat_levn,
		excit,
		pops,
		depart,
		/* data[ilo][iho] are net A's, from up to low, A * esc prob
		 * data[ihi][ilo] is col str or collision rate from ihi to ilo */
		&data,
		&dest,
		&pump,
		&CollRate_levn,
		create,
		destroy,
		/* say that we have evaluated the collision rates already */
		TRUE,
		/*&ipdest,*/
		&rot_cooling,
		&dCoolDT,
		" H2 ",
		/* lgNegPop positive if negative pops occured, negative if too cold */
		&lgNegPop,
		lgDeBug );/* option to print suff - set to true for debug printout */

	/* this will be sum of H2_populations of levels within correct bounds of matrix,
	 * not increased bounds we set above */
	sum_pops = 0.;
	/* copy populations back into 3D vector, also keep track of sum of populations */
	for( i=0; i< nXLevelsMatrix; ++i )
	{
		sum_pops += pops[i];
	}

	/* now correct H2_populations of levels that should have been in the matrix */
	factor =  abundance / SDIV(sum_pops );
	/*>>chng 04 aug 19, try using soln from matrix without renorm */
	factor = 1.;

	if(mole.lgH2_TRACE) 
		fprintf(ioQQQ," ratio of correct matrix pops to total is %.4e\n", factor );
	/*fprintf(ioQQQ,"DEBUG levpops");*/
	for( i=0; i< nXLevelsMatrix; ++i )
	{
		ip = H2_ipX_ener_sort[i];
		iRot = ipRot_H2_energy_sort[ip];
		iVib = ipVib_H2_energy_sort[ip];
		/* correct all H2_populations, since same scale factor is off for all */
		pops[i] *= factor;
		/* >>chng 04 sep 08, retail old pops, then update new, gs */
		/*H2_populations[0][iVib][iRot] = (float)pops[i];
		H2_old_populations[0][iVib][iRot] = (float)pops[i];*/
		H2_old_populations[0][iVib][iRot] = H2_populations[0][iVib][iRot];
		H2_populations[0][iVib][iRot] = (float)pops[i];
		/*fprintf(ioQQQ,"\t%.3e",pops[i]/sum_pops);*/
	}
	/*fprintf(ioQQQ,"\n",pops[i]/sum_pops);*/

	if(mole.lgH2_TRACE) 
	{
		/* print pops that came out of matrix */
		fprintf(ioQQQ,"\tH2_Level_lowJ hmi.H2_total: %.3e matrix rel pops\n",hmi.H2_total);
		fprintf(ioQQQ,"v\tJ\tpop\n");
		for( i=0; i<nXLevelsMatrix; ++i )
		{
			ip = H2_ipX_ener_sort[i];
			iRot = ipRot_H2_energy_sort[ip];
			iVib = ipVib_H2_energy_sort[ip];
			fprintf(ioQQQ,"%3li\t%3li\t%.3e\n",
				iVib , iRot , H2_populations[0][iVib][iRot]/hmi.H2_total);
		}
	}

	if( lgNegPop > 0 )
	{
		fprintf(ioQQQ,"H2_Level_low_matrix called atom_levelN which returned negative H2_populations.\n");
	}

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

/*H2_collid_rates - set H2 collision rates */
void H2_collid_rates( void )
{
	long int numb_coll_trans = 0;
	double excit;
	double t = phycon.te/1000. + 1.;
	double t2 = POW2(t);
	long int iElecHi , iElecLo , ipHi , iVibHi , iRotHi , 
		ipLo , iVibLo , iRotLo , nColl;

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


	iElecHi = 0;
	iElecLo = 0;
	if(mole.lgH2_TRACE) 
		fprintf(ioQQQ,"H2 set collision rates\n");
	/* set true to print all collision rates then quit */
#		define PRT_COLL	FALSE
	/* loop over all possible bound-bound collisional changes within X 
		* and set collision rates, which only depend on Te
		* will go through array in energy order since coll trans do not
		* correspond to a line */
	H2_coll_dissoc_rate_coef[0][0] = 0.;
	for( ipHi=1; ipHi<nLevels_per_elec[0]; ++ipHi )
	{
		double energy;

		/* obtain the proper indices for the upper level */
		long int ip = H2_ipX_ener_sort[ipHi];
		iVibHi = ipVib_H2_energy_sort[ip];
		iRotHi = ipRot_H2_energy_sort[ip];

		/* this is a guess of the collisional dissociation rate coeffient -
			* will be multiplied by the sum of all colliders */
		energy = H2_DissocEnergies[0] - energy_wn[0][iVibHi][iRotHi];
		ASSERT( energy > 0. );
		/* we made this up - boltzmann factor times rough coefficient */
		H2_coll_dissoc_rate_coef[iVibHi][iRotHi] = 
			1e-14f * (float)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll;
		/*fprintf(ioQQQ,"coll_dissoc_rateee\t%li\t%li\t%.3e\t%.3e\n",
			iVibHi,iRotHi,energy, H2_coll_dissoc_rate_coef[iVibHi][iRotHi] );*/

		for( ipLo=0; ipLo<ipHi; ++ipLo )
		{
			/* these are fits to the existing collision data */
			double gbarcoll[N_X_COLLIDER][3] = 
			{   
				{-9.9265 , -0.1048 , 0.456  },
				{-8.281  , -0.1303 , 0.4931 },
				{-10.0357, -0.0243 , 0.67   },
				{-8.6213 , -0.1004 , 0.5291 },
				{-9.2719 , -0.0001 , 1.0391 } };

			ip = H2_ipX_ener_sort[ipLo];
			iVibLo = ipVib_H2_energy_sort[ip];
			iRotLo = ipRot_H2_energy_sort[ip];

			/* >>chng 04 sep 14, do all levels */
			/* do not do the very large v levels - space not even allocated */
			/*if( iVibLo>VIB_COLLID || iVibHi>VIB_COLLID )
				continue;*/

			ASSERT( energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] > 0.);

			/* in following the colliders are H, He, H2(ortho), H2(para), and H+ */
			/* fits were read in from the following files: "H2_coll_H.dat" ,
				* "H2_coll_He.dat" , "H2_coll_H2ortho.dat" ,"H2_coll_H2para.dat",
				* "H2_coll_Hp.dat" */

			/* keep track of number of different collision routes */
			++numb_coll_trans;
			/* this is sum over all different colliders, except last one which is special,
				* linear rather than log formula for that one */
			for( nColl=0; nColl<N_X_COLLIDER-1; ++nColl )
			{
				/* this branch - real collision rate coefficients exist, use them */
				if( CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][0]!= 0 )
				{
					double r = CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][0] + 
							    CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][1]/t + 
							    CollRateFit[nColl][iVibHi][iRotHi][iVibLo][iRotLo][2]/t2;
					H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] = 
						(float)pow(10.,r)*mole.lgColl_deexec_Calc;
					if( PRT_COLL )
						fprintf(ioQQQ,"col fit\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
							nColl,
							iVibHi,iRotHi,iVibLo,iRotLo,
							energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
							H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
				}
				/* this is option to use guess of collision rate coefficient - but only if this is 
					* a downward transition that does not mix ortho and para */
				/* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
				else if( mole.lgColl_gbar  && 
					(H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) )
				{
					/* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient,
					* and x is the energy in wavenumbers */
					double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
					/* do not let energy diff be smaller than 100 wn, the smallest
						* diff for which we fit the rate coefficients */
					ediff = MAX2(100., ediff );
					H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] = 
						(float)pow(10. ,
						gbarcoll[nColl][0] + gbarcoll[nColl][1] * 
						pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;

					if( PRT_COLL )
						fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
							nColl+10,
							iVibHi,iRotHi,iVibLo,iRotLo,
							energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
							H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
				}
			}
			/* last one is special, with a different fit - this is the collision of H2 with
				* protons - of this group, cause ortho - para conversion */
			/* >>refer	H2	coll Hp	Gerlich, D., 1990, J. Chem. Phys., 92, 2377-2388 */
			if( CollRateFit[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo][1] != 0 )
			{
				H2_CollRate[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo] = 
					CollRateFit[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo][0] * 1e-10f *
					/* sec fit coef was dE in milli eV */
					(float)sexp( CollRateFit[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo][1]/1000./phycon.te_eV)*mole.lgColl_deexec_Calc;
				if( PRT_COLL )
					fprintf(ioQQQ,"col fit\t%i\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
						N_X_COLLIDER-1,
						iVibHi,iRotHi,iVibLo,iRotLo,
						energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
						H2_CollRate[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo] );
			}
			/* this is option to use guess of rate coefficient for ortho-para
				* conversion by collision with protons */
			/* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
			else if( mole.lgColl_gbar )
			{
				/* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient,
					* and x is the energy in wavenumbers */
				double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
				ediff = MAX2(100., ediff );
				H2_CollRate[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo] = 
					(float)pow(10. ,
					gbarcoll[N_X_COLLIDER-1][0] + gbarcoll[N_X_COLLIDER-1][1] * 
					pow(ediff ,gbarcoll[N_X_COLLIDER-1][2])	)*mole.lgColl_deexec_Calc;

				if( PRT_COLL )
					fprintf(ioQQQ,"col gbr\t%i\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
						N_X_COLLIDER-1+10,
						iVibHi,iRotHi,iVibLo,iRotLo,
						energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
						H2_CollRate[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo] );
			}
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC )
				{
					fprintf(ioQQQ,"bugcoll\tiVibHi\t%li\tiRotHi\t%li\tiVibLo\t%li\tiRotLo\t%li\tcoll\t%.2e\n",
						iVibHi,iRotHi,iVibLo,iRotLo,
						H2_CollRate[N_X_COLLIDER-1][iVibHi][iRotHi][iVibLo][iRotLo] );
				}
			}
		}
	}

	/* at this stage the collision rates that came in from the large data files
		* have been entered into the H2_CollRate array.  Now add on three extra collision
		* terms, the ortho para atomic H collision rates from
		* >>>refer	H2	collision	Sun, Y., & Dalgarno, A., 1994, ApJ, 427, 1053-1056 
		*/
	nColl = 0;
	iElecHi = 0;
	iElecLo = 0;
	iVibHi = 0;
	iVibLo = 0;

	/* >>chng 02 nov 13, the sun and dalgarno rates diverge to + inf below this temp */
	if( phycon.te >= 100. )
	{
		double excit1;
		/* this is the J=1-0 d
		ownward collision rate */
		iRotLo = 0;
		iRotHi = 1;
		excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/phycon.te);
		excit = sexp( -(POW2(5.30-460./phycon.te)-21.2) )*1e-13;

		H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo] = (float)(
			excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
			H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi / 
			/* >>chng 02 nov 13, from 2nd to first */
			SDIV(excit1) )*mole.lgColl_deexec_Calc;
			/*sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/phycon.te))*mole.lgColl_deexec_Calc;*/

		if( PRT_COLL )
			fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
				nColl,
				iVibHi,iRotHi,iVibLo,iRotLo,
				energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
				H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );

		/* this is the J=3-0 downward collision rate */
		iRotLo = 0;
		iRotHi = 3;
		excit = sexp( -(POW2(6.36-373./phycon.te)-34.5) )*1e-13;
		excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/phycon.te);
		H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo] = (float)(
			excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
			H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi / 
			SDIV(excit1) )*mole.lgColl_deexec_Calc;

		if( PRT_COLL )
			fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
				nColl,
				iVibHi,iRotHi,iVibLo,iRotLo,
				energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
				H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );

		/* this is the downward J=2-1 collision rate */
		iRotLo = 1;
		iRotHi = 2;
		excit = sexp( -(POW2(5.35-454./phycon.te)-23.1 ) )*1e-13;
		excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/phycon.te);
		H2_CollRate[0][iVibHi][iRotHi][iVibLo][iRotLo] = (float)(
			excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
			H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi / 
			SDIV(excit1) )*mole.lgColl_deexec_Calc;

		if( PRT_COLL )
			fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
				nColl,
				iVibHi,iRotHi,iVibLo,iRotLo,
				energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
				H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo] );
	}
	else
	{
		H2_CollRate[0][iVibHi][1][iVibLo][0] = 0.;
		H2_CollRate[0][iVibHi][3][iVibLo][0] = 0.;
		H2_CollRate[0][iVibHi][2][iVibLo][1] = 0.;
	}

	if( mole.lgH2_TRACE )
		fprintf(ioQQQ,
		" collision rates updated for new temp, number of trans is %li\n",
		numb_coll_trans);

	return;

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

/* do level H2_populations for H2, called by Hydrogenic after ionization and H chemistry
 * has been recomputed */
void H2_LevelPops( void )
{
	static float TeUsedBoltz=-1.f;
	static float TeUsedColl=-1.f;
	float H2_renorm_conserve=0. , 
		H2_renorm_conserve_init=0. , 
		sumold;
	static double part_fun = 0.;
	double oldrate;
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
	long int nColl,
		i;
	static int lgPopsSetZero=FALSE;
	long int n_pop_oscil = 0;
	int kase=0;
	int lgConv_h2_soln,
		lgPopsConv,
		lgHeatConv,
		lgSolomonConv,
		lgOrthoParaRatioConv;
	double quant_old=-1.,
		quant_new=-1.;

	double H2s_pump;
	double ratio_H2s_pump;

	int lgH2_pops_oscil=FALSE,
		lgH2_pops_ever_oscil=FALSE;
	long int nEner,
		ipHi, ipLo;
	long int iElec , iVib , iRot,ip;
	float sum_pops_matrix;
	float colldown,
		collup;
	/* old and older ortho - para ratios, used to determine whether soln is converged */
	static double ortho_para_old=0. , ortho_para_older=0. , ortho_para_current=0.;
	float frac_new_oscil=0.25;
	/* flags indicating whether non-positive level H2_populations occurred
	int lgNegPop , lgZeroPop; */
	/* keep track of changes in population */
	double PopChgMax=0. , PopChgMaxOld=0.;
	long int iRotMaxChng , iVibMaxChng;
	/* reason not converged */
	char chReason[100];

	/* H2 not on, so space not allocated and return,
	 * also return if calculation has been declared a failure */
	if( !h2.lgH2ON || conv.lgAbort )
		return;

	if(mole.lgH2_TRACE) 
		fprintf(ioQQQ,
		"\n***************H2_LevelPops call %li this iteration, zone is %.2f, H2/H is %.2e\n", 
		nCallH2_this_iteration,
		fnzone,
		hmi.H2_total/dense.gas_phase[ipHYDROGEN]
		);

	/* if H2 fraction is small then just zero out H2_populations and cooling, and return,
	 * but, if H2 has ever been done, redo irregarless of abundance -
	 * if large H2 is ever evaluated then mole.H2_to_H_limit is ignored */
	if( (!hmi.lgBigH2_evaluated && hmi.H2_total/dense.gas_phase[ipHYDROGEN] < mole.H2_to_H_limit )
		|| hmi.H2_total < 1e-20 )
	{
		/* will not compute the molecule */
		if(mole.lgH2_TRACE) 
			fprintf(ioQQQ,
			"  H2_LevelPops pops too small, not computing, H2/H is %.2e and mole.H2_to_H_limit is %.2e.",
			hmi.H2_total/dense.gas_phase[ipHYDROGEN] ,
			mole.H2_to_H_limit);
		/* no need to zero if already zeroed */
		if( lgPopsSetZero )
		{
			if(mole.lgH2_TRACE) 
				fprintf(ioQQQ,
				" pops already zero, just return\n");
			return;
		}
		if(mole.lgH2_TRACE) 
			fprintf(ioQQQ,
			" zero out pops.\n");
		/* set pops to zero and set flag saying not to do this again until pops set non-zero */
		lgPopsSetZero = TRUE;
		if(mole.lgH2_TRACE) fprintf(ioQQQ," zero and return\n");
		/* zero everything out - loop over all possible lines */
		for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
		{
			pops_per_elec[iElecHi] = 0.;
			for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
			{
				pops_per_vib[iElecHi][iVibHi] = 0.;
				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, change !=0 to >0 */
								if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
								{
									/* population of lower level with correction for stim emission */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc = 0.;

									/* population of lower level */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo = 0.;

									/* population of upper level */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi = 0.;

									/* following two heat exchange excitation, deexcitation */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cool = 0.;
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].heat = 0.;

									/* intensity of line */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity = 0.;

									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots = 0.;
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ots = 0.;
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot = 0.;
								}
							}
						}
					}
				}
			}
		}
		hmi.H2_star_BigH2 = 0.;
		hmi.H2_photodissoc_BigH2 = 0.;
		hmi.HeatH2Dish_BigH2 = 0.;
		hmi.HeatH2Dexc_BigH2 = 0.;
		hmi.deriv_HeatH2Dexc_BigH2 = 0.;
		hmi.H2_Solomon_dissoc_rate_BigH2_total = 0.;
		hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
		hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
		hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
		/* end of zero abundance branch */
		return;
	}

	/* we will set the pops to something, remember this */
	lgPopsSetZero = FALSE;

	/* check whether we need to update the H2_Boltzmann factors, LTE level H2_populations,
	 * and partition function.  lte level pops normalized by partition function,
	 * so sum of pops is unity */

	/* say that H2 has been computed, ignore previous limit to abund
	 * in future - this is to prevert oscillations as model is engaged */
	hmi.lgBigH2_evaluated = TRUE;

	/* do we need to update the boltzmann factors? */
	/*lint -e777 float test equality */
	if( phycon.te != TeUsedBoltz )
	/*lint +e777 float test equality */
	{
		part_fun = 0.;
		TeUsedBoltz = phycon.te;
		/* loop over all levels setting H2_Boltzmann and deriving partition function */
		for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
		{
			for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
			{
				for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
				{
					H2_Boltzmann[iElec][iVib][iRot] = 
						sexp( energy_wn[iElec][iVib][iRot] / phycon.te_wn );
					/* sum the partition funciton - boltzmann factor times statistical weight */
					part_fun += H2_Boltzmann[iElec][iVib][iRot] * H2_stat[iElec][iVib][iRot];
					ASSERT( part_fun > 0 );
				}
			}
		}
		/* have partition function, set H2_populations_LTE (populations for unit H2 density) */
		for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
		{
			for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
			{
				for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
				{
					/* these are the H2 LTE populations for a unit H2 density -
					 * these populations will sum up to unity */
					H2_populations_LTE[iElec][iVib][iRot] = 
						H2_Boltzmann[iElec][iVib][iRot] * 
						H2_stat[iElec][iVib][iRot] / part_fun;
					/*if( iElec==0 && iVib < 2)
						fprintf(ioQQQ,"DEBUG lte pop\t%i\t%i\t%e\n",
						iVib,iRot,H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total);*/
				}
			}
		}
		if(mole.lgH2_TRACE) 
			fprintf(ioQQQ,
			"H2 set H2_Boltzmann factors, T=%.2f, partition function is %.2f\n",
			phycon.te,
			part_fun);
	}
	/* end loop setting H2_Boltzmann factors, partition function, and lte H2_populations */

	/* check whether we need to update the collision rates */
	if( fabs(1. - TeUsedColl / phycon.te ) > 0.05  )
	{
		H2_collid_rates();
		TeUsedColl = phycon.te;
	}
	if( PRT_COLL )
		exit(98);
	/* end loop setting collision rates */

	/* set the H2_populations when this is the first call to this routine on 
	 * current iteration- will use LTE H2_populations  */
	if( nCallH2_this_iteration==0 )
	{
		/* very first call so use lte H2_populations */
		if(mole.lgH2_TRACE) fprintf(ioQQQ,"H2 1st call - using lte level pops\n");
		for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
		{
			pops_per_elec[iElec] = 0.;
			for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
			{
				pops_per_vib[iElec][iVib] = 0.;
				for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
				{
					/* LTE populations are for unit H2 density, so need to multiply
					 * by total H2 density */
					H2_old_populations[iElec][iVib][iRot] = 
						(float)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total;
					H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot];
				}
			}
		}
		/* first guess at ortho and para densities */
		h2.ortho_density = 0.75*hmi.H2_total;
		h2.para_density = 0.25*hmi.H2_total;
		ortho_para_current = h2.ortho_density / SDIV( h2.para_density );
		ortho_para_older = ortho_para_current;
		ortho_para_old = ortho_para_current;
		/* this is the fraction of the H2 pops that are within the levels done with a matrix */
		frac_matrix = 1.;
	}
	/* end loop setting very first LTE H2_populations */

#	define	FRAC_OLD	(0.5f)
	/* find sum of all H2_populations in X */
	iElec = 0;
	pops_per_elec[0] = 0.;
	for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
	{
		pops_per_vib[0][iVib] = 0.;
		for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
		{
			pops_per_elec[0] += H2_populations[iElec][iVib][iRot];
			pops_per_vib[0][iVib] += H2_populations[iElec][iVib][iRot];
		}
	}
	/* now renorm the old populations to the correct current H2 density, 
	 * At this point pops_per_elec[0] (pop of X)
	 * is the result of the previous evaluation of the H2 population, 
	 * following is the ratio of the current chemistry soln H2 to the previous H2 */
	H2_renorm_chemistry = hmi.H2_total/pops_per_elec[0];
	if(mole.lgH2_TRACE) 
		fprintf(ioQQQ,
			"H2 H2_renorm_chemistry is %.4e, hmi.H2_total is %.4e pops_per_elec[0] is %.4e\n",
			H2_renorm_chemistry ,
			hmi.H2_total,
			pops_per_elec[0]);

	/* renormalize all level populations for the current chemical solution */
	for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
	{
		for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
		{
			H2_populations[iElec][iVib][iRot] *= H2_renorm_chemistry;
			H2_old_populations[iElec][iVib][iRot] = H2_populations[iElec][iVib][iRot];
		}
	}
	ASSERT( fabs(h2.ortho_density+h2.para_density-hmi.H2_total)/hmi.H2_total < 0.001 );

	if(mole.lgH2_TRACE)
		fprintf(ioQQQ,
		" H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n",
		pops_per_elec[0],
		hmi.H2_total);

	/* these are the colliders that will be considered as depopulating agents */
	/* the colliders are H, He, H2 ortho, H2 para, H+ */
	/* atomic hydrogen */
	collider_density[0] = dense.xIonDense[ipHYDROGEN][0];
	/* atomic helium */
	collider_density[1] = dense.xIonDense[ipHELIUM][0];
	/* all ortho h2 */
	collider_density[2] = (float)h2.ortho_density;
	/* all para H2 */
	collider_density[3] = (float)h2.para_density;
	/* protons - ionized hydrogen */
	collider_density[4] = dense.xIonDense[ipHYDROGEN][1];
	/* 02 oct 13 add this */
	/* assume that H3+ has precisely same effects as proton */
	collider_density[4] += hmi.Hmolec[ipMH3p];

	/* this is total density of all collideres, is only used for collisional dissociation */
	collider_density_total = 0.;
	for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
	{
		collider_density_total += collider_density[nColl];
	}
	/* >>chng 02 oct 24, add electrons since have roughly same rate as others
	 * for collisional dissociation */
	collider_density_total += (float)dense.eden;

	if(mole.lgH2_TRACE)
	{
		fprintf(ioQQQ," Collider densities are:");
		for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
		{
			fprintf(ioQQQ,"\t%.3e", collider_density[nColl]);
		}
		fprintf(ioQQQ,"\n");
	}

	/* update state specific rates in H2_X_formation (cm-3 s-1) that H2 forms from grains and H- */
	mole_H2_form();

#	if 0
	/* rate of entry into X from H- and formation on grain surfaces 
	 * will one of several distribution functions derived elsewhere
	 * first zero out formation rates and rates others collide into particular level */
	for( iVib = 0; iVib <= nVib_hi[0] ; ++iVib )
	{
		for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
		{
			/* this will be the rate formation (s-1) of H2 due to
			 * both formation on grain surfaces and the H minus route,
			 * also H2+ + H => H2 + H+ into one vJ level */
			/* units cm-3 s-1 */
			H2_X_formation[iVib][iRot] = 0.;
			/* this is the total collisional rate into a vJ level within X due to
			 * all collisional processes, including cosmic rays and thermal collisions */
			/* cm-3 s-1 */
			H2_col_rate_in[iVib][iRot] = 0.;
		}
	}

	/* loop over all grain types, finding total formation rate into each ro-vib level,
	 * also keeps trace of total formation into H2 ground and star, as defined by Tielens & Hollenbach,
	 * these are used in the H molecular network */
	rate = gv.rate_h2_form_grains_used_total;
	hmi.H2_forms_grain = 0.;
	hmi.H2star_forms_grain = 0.;
	/* >>chng 02 oct 08, resolved grain types */
	for( nd=0; nd<gv.nBin; ++nd )
	{
		int ipH2 = (int)gv.which_H2distr[gv.bin[nd]->matType];
		for( iVib = 0; iVib <= nVib_hi[0] ; ++iVib )
		{
			for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
			{
				/* >>chng 02 nov 14, changed indexing into H2_X_grain_formation_distribution and gv.bin, PvH */
				float one = 
					/* H2_X_grain_formation_distribution is normalized to a summed total of unity */
					H2_X_grain_formation_distribution[ipH2][iVib][iRot] * 
					/* units of following are s-1 */
					(float)gv.bin[nd]->rate_h2_form_grains_used;
				/* final units are s-1 */
				/* units cm-3 s-1 */
				/* >>chng 04 may 05, added atomic hydrogen density, units cm-3 s-1 */
				H2_X_formation[iVib][iRot] += one*dense.xIonDense[ipHYDROGEN][0];

				/* save rates for formation into "H2" and "H2*" in the chemical
				 * network - it resolves the H2 into two species, as in 
				 * Hollenbach / Tielens work - these rates will be used in the
				 * chemistry solver to get H2 and H2* densities */
				if( energy_wn[0][iVib][iRot] < ENERGY_H2_STAR )
				{
					hmi.H2_forms_grain += one;
				}
				else
				{
					hmi.H2star_forms_grain += one;
				}
			}
		}
	}

	/* formation of H2 in excited states from H- H minus */
	/* >>chng 02 oct 17, temp dependent fit to rate, updated reference,
	 * about 40% larger than before */
	/* rate in has units cm-3 s-1 */
	rate = hmi.Hmolec[ipMHm] * dense.xIonDense[ipHYDROGEN][0] * hmi.assoc_detach;
	/*rate = hmi.hminus*1.35e-9f;*/
	/* convert to dimensionless factors that add to unity */
	/* >>chng 02 oct 17, use proper distribution function */
	hmi.H2star_forms_hminus = 0.;
	hmi.H2_forms_hminus = 0.;
	oldrate = 0.;
	/* which temperature point to use? */
	if( phycon.alogte<=1. )
	{
		ipT = 0;
		frac_lo = 1.;
		frac_hi = 0.;
	}
	else if( phycon.alogte>=4. )
	{
		ipT = nTE_HMINUS-2;
		frac_lo = 0.;
		frac_hi = 1.;
	}
	else
	{
		/* find the temp */
		for( ipT=0; ipT<nTE_HMINUS-1; ++ipT )
		{
			if( H2_te_hminus[ipT+1]>phycon.alogte )
				break;
		}
		frac_hi = (phycon.alogte-H2_te_hminus[ipT])/(H2_te_hminus[ipT+1]-H2_te_hminus[ipT]);
		frac_lo = 1.-frac_hi;
	}

	/* we now know how to interpolate, now fill in H- formation sites */
	for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
	{
		for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
		{
			/* the temperature-interpolated distribution function, adds up to unity, 
			 * dimensionless */
			double rate_interp =
				frac_lo*H2_X_hminus_formation_distribution[ipT][iVib][iRot] +
				frac_hi*H2_X_hminus_formation_distribution[ipT+1][iVib][iRot];

			/* above rate was set, had dimensions cm-3 s-1 */
			float one = (float)(rate * rate_interp);

			/* units cm-3 s-1 */
			H2_X_formation[iVib][iRot] += one;

			oldrate += rate_interp;

			/* save rates to pass back into molecule network */
			if( energy_wn[0][iVib][iRot] < ENERGY_H2_STAR )
			{
				hmi.H2_forms_hminus += one;
			}
			else
			{
				hmi.H2star_forms_hminus += one;
			}
		}
	}
	/* confirm that shape function is normalized correctly */
	ASSERT( fabs(1.-oldrate)<1e-4 );

	/* >>chng 03 feb 10, add this population process */
	/* H2+ + H => H2 + H+,
	 * >>refer	H2	population	Krstic, P.S., preprint 
	 * all goes into v=4 but no J information, assume into J = 0 */
	/* >>chng 04 may 05, add density at end */
	rate = hmi.bh2h2p * hmi.Hmolec[ipMH2p] * dense.xIonDense[ipHYDROGEN][0];
	iVib = 4;
	iRot = 0;
	/* units cm-3 s-1 */
	H2_X_formation[iVib][iRot] += (float)rate;
#	endif

	/* this flag will say whether H2 H2_populations have converged,
	 * by comparing old and new values */
	lgConv_h2_soln = FALSE;
	/* this will count number of passes around following loop */
	loop_h2_pops = 0;
#	if 0
	SolomonRateMax = -1.f;
	nRot_SolomonRateMax = -1;
	nVib_SolomonRateMax = -1;
	nRotHi_SolomonRateMax = -1;
	nVibHi_SolomonRateMax = -1;
	nEkcHi_SolomonRateMax = -1;
#	endif
	{
		static long int nzoneEval=-1;
		if( nzone != nzoneEval )
		{
			nzoneEval = nzone;
			/* this is number of zones with H2 solution in this iteration */
			++nH2_zone;
		}
	}

	/* begin - start level population solution
	 * first do electronic excited states, Lyman, Werner, etc
	 * using old solution for X
	 * then do matrix if used, then solve for pops of rest of X */
	/* >>chng 04 apr 06, sub number of oscillations from limit - don't waste loops 
	 * if solution is unstable */
	while( loop_h2_pops < LIM_H2_POP_LOOP-n_pop_oscil && !lgConv_h2_soln )
	{
		float rate_in , rate_out, rate_tot;
		static float old_HeatH2Dexc_BigH2=0. , HeatChangeOld=0. , HeatChange=0.;

		/* this is number of trips around loop this time */
		++loop_h2_pops;
		/* this is number of times through this loop in entire iteration */
		++nH2_pops;

		/* beging solution for electronic excited states
		 * loop over all possible pumping routes to excited electronic states
		 * to get radiative excitation and dissociation rates */
		hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
		H2s_pump = 0;
		ratio_H2s_pump = 0;
		rate_tot = 0;

		for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
		{
			/* this will be total population in each electronic state */
			pops_per_elec[iElecHi] = 0.;

			if(mole.lgH2_TRACE) fprintf(ioQQQ," Pop(e=%li):",iElecHi);
			for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
			{
				pops_per_vib[iElecHi][iVibHi] = 0.;
				/* =======================INSIDE POPULATIONS EXEC ELEC CONVERGE LOOP =====================*/
				for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
				{
					/* Solomon process done here,
					 * sum of all rates into and out of these upper levels */
					/* all inward rates have units cm-3 s-1 */
					rate_in = 0.;
					/* this term is spontaneous dissociation of excited elec states into 
					 * the X continuum */
					/* all outward rates have units s-1 */
					rate_out = H2_dissprob[iElecHi][iVibHi][iRotHi];

					/* now loop over all levels within X */
					iElecLo=0;
					for( iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
					{
						for( iRotLo=Jlowest[iElecLo]; iRotLo<=nRot_hi[iElecLo][iVibLo]; ++iRotLo )
						{
							/* >>chng 03 feb 14, change test from !=0 to > 0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								/* solv elec exc state, 
								 * rate lower level in X go to elec excited state, cm-3 s-1 */
								/* >>chng 04 may 7, use straight lya rate, rather than lya / 1e3 */
								/* >>chng 04 may 25, put lya cr excit back to x12 /1000 since others
								 * in leiden did not include it */
								float rate_one =
									H2_old_populations[iElecLo][iVibLo][iRotLo]* 
									(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump 
									/* add in collisional excitation by suprathermals */
									+hmi.lgLeidenCRHack*secondaries.x12tot /1e3f );
								/*TODO	0	get better relationship between cr excit of Lya rate */

								/* this is a permitted electronic transition, must presurve nuclear spin */
								ASSERT( H2_lgOrtho[iElecHi][iVibHi][iRotHi]-H2_lgOrtho[iElecLo][iVibLo][iRotLo]==0 );

								/* this is the rate [cm-3 s-1] electrons move into the upper level from X */
								rate_in += rate_one;

								/* this is total X -> exc elec state rate, cm-3 s-1 */
								rate_tot += rate_one;

								/* excitation rate for solomon process - this currently has units
								 * cm-3 s-1 but will be divided by total pop of X and become s-1 */
								/* >>chng 04 may 25, Gargi Shaw found this bug - had been total sum */
								/*hmi.H2_H2g_to_H2s_rate_BigH2 += rate_one/SDIV(hmi.H2_total);*/
								/* this has unit s-1 and will be used in hmole_step.c for H2g->H2s */
								if( energy_wn[0][iVibLo][iRotLo] > ENERGY_H2_STAR )
								{
									hmi.H2_H2g_to_H2s_rate_BigH2 += (H2_old_populations[iElecHi][iVibHi][iRotHi]*
										(H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Aul*
										/* escape and destruction */
										(H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pesc + 
										H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pelec_esc + 
										H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].Pdest) +
										/* induced emission down */
										H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].pump *
										H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gLo/
										H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].gHi))/SDIV(hmi.H2_total); /*rate_one/SDIV(hmi.H2_total);*/
								
									/* pumping from H2s, unit cm-3s-1*/
									H2s_pump += H2_old_populations[0][iVibLo][iRotLo]*(H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo].pump);
									/* Fraction of all X pumps that are from H2s, unitless*/
									ratio_H2s_pump = H2s_pump/SDIV(rate_tot);
								}

								/* solv elec excit states,
								 * this is the rate [s-1] electrons leave the excit elec upper level
								 * and decay into X */
								/* >>chng 04 apr 22, big bug, pump rate multiplied the A */
								rate_out +=
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*
									/* escape and destruction */
									 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc + 
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc + 
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest) +
									 /* induced emission down */
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump *
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo/
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi;

								ASSERT( (rate_in >=0.) && rate_out >= 0. );
							}
						}
					}

					/* population [cm-3] of the electronic excited state this only includes 
					 * radiative processes between X and excited electronic states - 
					 * collisions are neglected
					 * X is done below and includes collisions */
					if( 0 )
					{
						H2_populations[iElecHi][iVibHi][iRotHi] = 
							rate_in / SDIV( rate_out);
					}
					else
					{
						/* take mean of old and new pops to damp noise */
						H2_populations[iElecHi][iVibHi][iRotHi] = (H2_populations[iElecHi][iVibHi][iRotHi]+
							rate_in /SDIV( rate_out) ) / 2.f;
					}
					ASSERT( H2_populations[iElecHi][iVibHi][iRotHi] >= 0. && 
						H2_populations[iElecHi][iVibHi][iRotHi] <= hmi.H2_total );

					/* this is total pop in this vib state */
					pops_per_vib[iElecHi][iVibHi] += H2_populations[iElecHi][iVibHi][iRotHi];

					/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
				}
				/* end excit vib pops loop */
				if(mole.lgH2_TRACE) fprintf(ioQQQ,"\t%.2e",pops_per_vib[iElecHi][iVibHi]/hmi.H2_total);

				/* total pop in each elec state */
				pops_per_elec[iElecHi] += pops_per_vib[iElecHi][iVibHi];
			}
			/* end excited elec pops loop */
			if(mole.lgH2_TRACE) fprintf(ioQQQ,"\n");
		}

		/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
		/* now do lowest levels H2_populations with matrix, 
		 * these should be collisionally dominated */
		/*fprintf(ioQQQ,"DEBUG h2 matr\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n" ,
			H2_populations[0][0][5] ,
			H2_populations[0][0][6] ,
			H2_populations[0][0][7] ,
			H2_populations[0][0][8] ,
			H2_populations[0][1][0],
			H2_populations[0][1][1] ,
			H2_populations[0][1][2] ,
			H2_populations[0][1][3],
			H2_populations[0][0][9] );*/
		if( nXLevelsMatrix )
		{
			H2_Level_low_matrix(
				/* the total abundance - frac_matrix is fraction of pop that was in these
				 * levels the last time this was done */
				hmi.H2_total * (float)frac_matrix );
			/*{ static int neval=0;++neval;
				if( neval > 10 )
				{
					fprintf(ioQQQ,"stopping h2\n");
					exit(99);
				}
			}*/
		}

		/* above set pops of excited electronic levels and low levels within X - 
		 * now do something about excited levels within the X state */
		/* nLevels_per_elec is number of levels within elec 0 - so nEner is one
		 * beyond end of array here - but will be decremented at start of loop */
		nEner = nLevels_per_elec[0];
		iElec = 0;

		/* these will do convergence check */
		PopChgMaxOld = PopChgMax;
		PopChgMax = 0.;
		iElecHi = 0;
		iRotMaxChng =-1;
		iVibMaxChng = -1;

		/* =======================INSIDE X POPULATIONS CONVERGE LOOP =====================*/
		/* begin solving for X
		 * this is the main loop that determines H2_populations within X */
		/* units of all rates in are cm-3 s-1, all rates out are s-2  */
		while( (--nEner) >= nXLevelsMatrix )
		{
			/* array of energy sorted indices within X */
			ip = H2_ipX_ener_sort[nEner];
			iVib = ipVib_H2_energy_sort[ip];
			iRot = ipRot_H2_energy_sort[ip];

			/* count formation from grains and H- as a collisional formation process */
			/* cm-3 s-1, evaluated in mole_H2_form */
			H2_col_rate_in[iVib][iRot] = H2_X_formation[iVib][iRot];

			/* this represents collisional dissociation into continuum of X,
			 * rates are just guesses */
			H2_col_rate_out[iVib][iRot] = collider_density_total *
				H2_coll_dissoc_rate_coef[iVib][iRot] * mole.lgColl_deexec_Calc ;

			/* this is cosmic ray or secondary photodissociation into mainly H2+ */
			H2_col_rate_out[iVib][iRot] += secondaries.csupra[ipHYDROGEN][0]*2.02f * mole.lgColl_deexec_Calc;

			/* will become rate (cm-3 s-1) other levels have radiative transitions to here */
			H2_rad_rate_in[iVib][iRot] = 0.;

			/* >>chng 04 may 25, add this test on energy, had always
			 * added this term, which is small.
			 * this is continuum photodissociation into 2H */
			if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
			{ 
				H2_rad_rate_out[iVib][iRot] = rfield.flux_accum[H2_ipPhoto[iVib][iRot]-1]*0.25e-18f;
			}
			else
			{
				H2_rad_rate_out[iVib][iRot] = 0;
			}
			/* ortho para conversion, and deexcitation, while on grain surface,
			 * H2 + grain => H2 + grain, do not want to add rate to 0,0 */
			if( iVib!= 0 || iRot!=0 )
			{
				/* this is the ground v J = 1->0 transition */
				if( iVib == 0 && iRot == 1 )
				{
					H2_col_rate_out[iVib][iRot] += 
						/* H2 grain interactions
						* H2 ortho - para conversion on grain surface,
						* rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
						(float)hmi.rate_h2_allX_2_J0_grains * mole.lgColl_deexec_Calc ;
				}
				/* add both rates if not 0,0 or 0,1 */
				else if( iVib >0  || iRot > 1 )
				{
					H2_col_rate_out[iVib][iRot] += 
						/* H2 grain interactions
						* H2 ortho - para conversion on grain surface,
						* rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
						(float)(hmi.rate_h2_allX_2_J0_grains + 
						/* rate (s-1) all v,J levels go to 0, regardless of nuclear spin */
						hmi.rate_h2_ortho_para_conserve) * mole.lgColl_deexec_Calc ;

				}
			}

			/* excitation within X due to cosmic rays,
			 * csupra is total ionization rate, cr_rate is branching ratio */
			if( iVib < CR_VIB && !hmi.lgLeidenCRHack)
			{
				/* cr_rate[CR_X][CR_VIB][CR_J][CR_EXIT];*/
				if( iRot> 1 && iRot < CR_J )
				{
					/* collision J to J-2 */
					H2_col_rate_out[iVib][iRot] += 
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][0] * mole.lgColl_deexec_Calc ;
					H2_col_rate_in[iVib][iRot-2] += H2_old_populations[0][iVib][iRot]*
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][0] * mole.lgColl_deexec_Calc ;
				}

				if( iRot<CR_J-2 && iRot+2 <= nRot_hi[iElec][iVib])
				{
					H2_col_rate_out[iVib][iRot+2] += 
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][2] * mole.lgColl_deexec_Calc ;
					H2_col_rate_in[iVib][iRot+2] += H2_old_populations[0][iVib][iRot]*
						secondaries.csupra[ipHYDROGEN][0] * cr_rate[0][iVib][iRot][2] * mole.lgColl_deexec_Calc ;
				}
			}

			/* this sum is for X going into all electronic excited states */
			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 )
					{

						/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
						/* solv X, the rate electrons enter this state from excited elec states */
						/* >>chng 03 feb 14, change test from != 0 to > 0 */
						if( H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul > 0. )
						{
							float rateone;
							/* the spontaneous rate down, escape and destruction, this
							 * ignores stilulated emission - very good approximation */
							/* >>chng 04 apr 22, big bug, pump rate multiplied the A */

							/* >>chng 04 may 25, uxe s12tot over 1000, had been x12 */
							/*TODO this is factor for lya excit  do something about this */
							/* cosmic ray or supra thermal collisions into excited states */
							H2_col_rate_out[iVib][iRot] += hmi.lgLeidenCRHack*secondaries.x12tot/1e3f * mole.lgColl_deexec_Calc;

							/* this is a permitted electronic transition, must presurve nuclear spin */
							ASSERT( H2_lgOrtho[iElecHi][iVibHi][iRotHi]-H2_lgOrtho[iElec][iVib][iRot]==0 );

							/* elec excit state decay into X */
							rateone =
								H2_old_populations[iElecHi][iVibHi][iRotHi]*(
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul*
								/* escape and destruction */
								(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pesc + 
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pelec_esc + 
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pdest) + 
								/* induced emission down */
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump *
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gLo/
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].gHi);
							ASSERT( rateone >= 0. );
							/* units cm-3 s-1 */
							H2_rad_rate_in[iVib][iRot] += rateone;

							/* pump up units - cm-1 */
							H2_rad_rate_out[iVib][iRot] +=
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].pump;
							{
								/*@-redef@*/
								enum {DEBUG_LOC=FALSE};
								/*@+redef@*/
								if( DEBUG_LOC && iRot==1 && iVib==0 )
								{
									fprintf(ioQQQ,"DEBUG E vibhi\t%li\trothi\t%li\trate_in\t%.2e\tAul\t%.2e\tpop\t%.2e\n",
										iVibHi,iRotHi, rateone ,
										H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul,
										H2_old_populations[iElecHi][iVibHi][iRotHi] );
								}
							}

							/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
						}
					}
				}
			}

			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				/* this debug gives total rates in and out after the excited elec states
				 * have been dealt with */
				if( DEBUG_LOC && iVib == 4 && iRot==0 )
				{
					fprintf(ioQQQ,"DEBUG eleccc vib\t%li\trot\t%li\tin\t%.2e\tout\t%.2e\trel pop\t%.2e\n",
						iVib,iRot, H2_rad_rate_in[iVib][iRot] ,
						H2_rad_rate_out[iVib][iRot],
						H2_old_populations[iElec][iVib][iRot] );
				}
			}

			/* now sum over states within X which are higher than current state */
			iElecHi = 0;
			for( ipHi = nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi )
			{
				ip = H2_ipX_ener_sort[ipHi];
				iVibHi = ipVib_H2_energy_sort[ip];
				iRotHi = ipRot_H2_energy_sort[ip];
				/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
				/* the rate we enter this state from more highly excited states within X
				 * by radiative decays, which have delta J = 0 or 2 */
				/* note test on vib is needed - iVibHi<iVib, energy order ok and space not allocated */
				if( (fabs(iRotHi-iRot)==2 || iRotHi==iRot )&& (iVib<=iVibHi) )
				/* >>chng 03 feb 14, change to test on Aul */
				/*if( H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul>0. )*/
				{
					float rateone;
					rateone =
						H2_old_populations[iElecHi][iVibHi][iRotHi]*
						 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul*
						(H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pesc + 
						 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pelec_esc + 
						 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Pdest);
					ASSERT( rateone >=0 );

					/* units cm-3 s-1 */
					H2_rad_rate_in[iVib][iRot] += rateone;

					{
						/*@-redef@*/
						enum {DEBUG_LOC=FALSE};
						/*@+redef@*/
						if( DEBUG_LOC && iRot==7 && iVib==3 )
						{
							fprintf(ioQQQ,"DEBUG bug X vibhi\t%li\trothi\t%li\trate_in\t%.2e\tAul\t%.2e\tpop\t%.2e\tsum\n",
								iVibHi,iRotHi, rateone ,
								H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Aul,
								H2_old_populations[iElecHi][iVibHi][iRotHi] );
						}
					}
				}
				/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
				/* collisional interactions with upper levels within X */
				/* VIB_COLLID is the highest vib level that has coll data */
				if( iVibHi <= VIB_COLLID && iVib <= VIB_COLLID)
				{
					/*SPEEDUP - move to before sweep, in routine that does collision rates,
					 * define rates up and down for all transitions*/
					colldown = 0.;
					for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
					{
						/* downward collision rate, units s-1 */
						float colldown1 = 
							H2_CollRate[nColl][iVibHi][iRotHi][iVib][iRot]*collider_density[nColl] * 
							mole.lgColl_deexec_Calc ;
						ASSERT( colldown1>=0. );
						colldown += colldown1;
					}
					/* rate in from upper level, units cm-3 s-1 */
					H2_col_rate_in[iVib][iRot] += H2_old_populations[iElecHi][iVibHi][iRotHi]*colldown;
					/* convert real rate back into electron cs in case other parts of code
					 * need this for reference */
					if( iVibHi >= iVib )
					{
						H2Lines[iElec][iVibHi][iRotHi][iElec][iVib][iRot].ColUL = colldown;
						LineConvRate2CS( &H2Lines[iElec][iVibHi][iRotHi][iElec][iVib][iRot] , 
							H2Lines[iElec][iVibHi][iRotHi][iElec][iVib][iRot].ColUL);
					}
					/* ortho para conversion on grain surfaces, also ortho-ortho and para-para */
					/*TODO	2	is this right - two rates that add together?  or is
					 * J=1 to 0 only added for 0,1 to 0,0? */
					if( (iVib ==0) && (iRot <=1) &&
						(H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVib][iRot])==0 )
					{
						/* this is same nuclear spin */
						H2_col_rate_in[iVib][iRot] += H2_old_populations[iElecHi][iVibHi][iRotHi]*
							(float)hmi.rate_h2_ortho_para_conserve;
					}
					else if( iVib==0 && iRot== 0 )
					{
						/* this is everything down to v=0 && J=0 */
						H2_col_rate_in[iVib][iRot] += H2_old_populations[iElecHi][iVibHi][iRotHi]*
							(float)hmi.rate_h2_allX_2_J0_grains;
					}

					/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
					/* upward collisions out of this level - probably very small
					 * since dE >> kT - H2_Boltzmann is e(- nu/kT ) */
					collup = (float)(colldown *	
						H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElec][iVib][iRot] *
						H2_Boltzmann[iElecHi][iVibHi][iRotHi] /
						SDIV( H2_Boltzmann[iElec][iVib][iRot] ) );
					H2_col_rate_out[iVib][iRot] += collup;
					{
						/*@-redef@*/
						enum {DEBUG_LOC=FALSE};
						/*@+redef@*/
						if( DEBUG_LOC && iRot==1 && iVib==0 && iRotHi==3 && iVibHi==0 )
						{
							fprintf(ioQQQ,"DEBUG  C1 vibhi\t%li\trothi\t%li\trate_in\t%.2e\tpop\t%.2e\tsum\t%.2e\n",
								iVibHi,iRotHi, colldown , collup ,
								H2_col_rate_in[iVib][iRot] );
						}
						if( DEBUG_LOC && iRotHi==1 && iVibHi==0 && iRot==0 && iVib==0 )
						{
							fprintf(ioQQQ,"DEBUG  C1 vibhi\t%li\trothi\t%li\trate_in\t%.2e\tpop\t%.2e\tsum\t%.2e\n",
								iVibHi,iRotHi, 
								colldown , 
								collup ,
								H2_col_rate_in[iVib][iRot] );
						}
					}
				}

				/* this branch no coll data, but will make some up on the fly,
				* buy only if g-bar is on and not spin changing */
				/* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
				else if( mole.lgColl_gbar && 
					(H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVib][iRot]==0) )
				{
					/* "typical" rates for H, He, H2ortho, H2para, H+ */
					float csguess[N_X_COLLIDER]={1e-15f,1e-17f,1e-18f,1e-18f,3e-10f};

					colldown = 0.;
					for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
					{
						colldown +=
							csguess[nColl]*collider_density[nColl];
					}
					H2_col_rate_in[iVib][iRot] += H2_old_populations[iElecHi][iVibHi][iRotHi]*colldown;

					/* >>chng 04 may 6, back collision had not been included! */
					collup = (float)(colldown *	
						H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElec][iVib][iRot] *
						H2_Boltzmann[iElecHi][iVibHi][iRotHi] /
						SDIV( H2_Boltzmann[iElec][iVib][iRot] ) );
					H2_col_rate_out[iVib][iRot] += collup;

					{
						/*@-redef@*/
						enum {DEBUG_LOC=FALSE};
						/*@+redef@*/
						if( DEBUG_LOC && iRot==1 && iVib==0 )
						{
							fprintf(ioQQQ,"DEBUG  Cf vibhi\t%li\trothi\t%li\trate_in\t%.2e\tpop\t%.2e\tsum\t%.2e\n",
								iVibHi,iRotHi, colldown ,
								H2_old_populations[iElecHi][iVibHi][iRotHi] ,
								H2_col_rate_in[iVib][iRot]);
						}
					}
				}
			}
			/* =======================INSIDE POPULATIONS X CONVERGE LOOP =====================*/
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				/* this debug gives total rates after interactions with 
				 * higher levels within X have been included */
				if( DEBUG_LOC&&iRot==1 && iVib==0 )
				{
					fprintf(ioQQQ,"DEBUG aboveee vib\t%li\trot\t%li\tin\t%.2e\tout\t%.2e\trel pop\t%.2e\n",
						iVib,iRot, H2_col_rate_in[iVib][iRot] ,
						H2_col_rate_out[iVib][iRot],
						H2_old_populations[iElec][iVib][iRot]/hmi.H2_total );
				}
			}

			/* we now have total rate this state is populated from above, now get rate
			 * this state intereacts with levels that are below */
			iElecLo = 0;
			for( ipLo = 0; ipLo<nEner; ++ipLo )
			{
				ip = H2_ipX_ener_sort[ipLo];
				iVibLo = ipVib_H2_energy_sort[ip];
				iRotLo = ipRot_H2_energy_sort[ip];
				/* radiative interactions between this level and lower levels */
				/* the test on vib is needed - the energies are ok but the space does not exist */
				if( ((fabs(iRotLo-iRot) == 2)||(fabs(iRotLo-iRot) == 0))  && (iVibLo<=iVib)/**/)
				{
					H2_rad_rate_out[iVib][iRot] +=
						H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Aul*
						(H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Pesc + 
						H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Pelec_esc + 
						H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Pdest);
				}
				/* =======================INSIDE X POPULATIONS CONVERGE LOOP =====================*/
				/* collisions to level below us - 
				* VIB_COLLID is the highest vib level that has coll data */
				if( iVib <= VIB_COLLID && iVibLo <= VIB_COLLID)
				{
					colldown = 0.;
					for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
					{
						colldown +=
							H2_CollRate[nColl][iVib][iRot][iVibLo][iRotLo]*collider_density[nColl];
					}
					/* s-1 */
					colldown *= mole.lgColl_deexec_Calc ;
					H2_col_rate_out[iVib][iRot] += colldown;
					/* rate up into this level from below */
					collup = (float)(colldown* H2_old_populations[iElecLo][iVibLo][iRotLo] *	
						H2_stat[iElec][iVib][iRot] / H2_stat[iElecLo][iVibLo][iRotLo] *
						H2_Boltzmann[iElec][iVib][iRot] /
						SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] ));
					/* units cm-3 s-1 */
					H2_col_rate_in[iVib][iRot] += collup ;

					{
						/*@-redef@*/
						enum {DEBUG_LOC=FALSE};
						/*@+redef@*/
						if( DEBUG_LOC && iRot==3 && iVib==0 && iRotLo==1 && iVibLo==0 )
						{
							fprintf(ioQQQ,"DEBUG  C2 vibhi\t%li\trothi\t%li\trate_in\t%.2e\tpop\t%.2e\tsum\t%.2e\n",
								iVib,iRot, colldown , collup ,
								H2_col_rate_in[iVib][iRot] );
						}
						/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
						if( DEBUG_LOC && iRot==1 && iVib==0 && iRotLo==0 && iVibLo==0 )
						{
							fprintf(ioQQQ,"DEBUG  C2 vibhi\t%li\trothi\t%li\trate_in\t%.2e\tpop\t%.2e\tsum\t%.2e\n",
								iVib,iRot, 
								colldown , 
								collup ,
								H2_col_rate_in[iVib][iRot] );
						}
					}
				}
				/* >>chng 02 jul 29, add wild guess if lower level is close but delta v is large */
				/* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
				/*TODO	2	- did all gbar coll include test on spin changing?  Also, this
				 * does not have the back collision term */
				else if( mole.lgColl_gbar && 
					(H2_lgOrtho[0][iVib][iRot]-H2_lgOrtho[0][iVibLo][iRotLo]==0) &&
					((energy_wn[iElec][iVib][iRot]-energy_wn[0][iVibLo][iRotLo])<2.*phycon.te)
					)
				{
					/* "typical" rates for H, He, H2ortho, H2para, H+ */
					float csguess[N_X_COLLIDER]={1e-15f,1e-17f,1e-18f,1e-18f,3e-10f};

					colldown = 0.;
					for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
					{
						colldown +=
							csguess[nColl]*collider_density[nColl];
					}
					/* cm-3 s-1 */
					colldown *= mole.lgColl_deexec_Calc ;
					H2_col_rate_out[iVib][iRot] += colldown;
					/* rate up into this level from below, units s-1 */
					collup = (float)(colldown* H2_old_populations[iElecLo][iVibLo][iRotLo] *	
						H2_stat[iElec][iVib][iRot] / H2_stat[iElecLo][iVibLo][iRotLo] *
						H2_Boltzmann[iElec][iVib][iRot] /
						SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] ));
					/* units cm-3 s-1 */
					H2_col_rate_in[iVib][iRot] += collup ;
				}
			}
			/* =======================INSIDE X POPULATIONS CONVERGE LOOP =====================*/

			/* we now have the total rates into and out of this level, get its population */
			if( (H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[iVib][iRot])> SMALLFLOAT )
			{
				/* H2_populations has units cm-3 */
				H2_populations[iElec][iVib][iRot] = 
					(H2_col_rate_in[iVib][iRot]+ H2_rad_rate_in[iVib][iRot]) / 
					(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[iVib][iRot]) ;

			}
			else
			{
				H2_populations[iElec][iVib][iRot] = 0.;
			}
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
				if( DEBUG_LOC && iRot==1 && iVib==0 )
				{
					fprintf(ioQQQ,"DEBUG  vib\t%li\trot\t%li\tin\t%.2e\tout\t%.2e\trel pop\t%.2e\n",
						iVib,iRot, H2_col_rate_in[iVib][iRot] ,
						H2_col_rate_out[iVib][iRot],
						H2_populations[iElec][iVib][iRot]/hmi.H2_total );
				}
			}

			ASSERT( H2_populations[iElec][iVib][iRot] >= 0.  );
		}
		/* find ortho and para densites, sum of pops in each vib */
		iElecHi = 0;
		if(mole.lgH2_TRACE) 
		{
			fprintf(ioQQQ," Pop(e=%li):",iElecHi);
		}

		/* this will become total pop is X, which will be renormed to equal hmi.H2_total */
		pops_per_elec[0] = 0.;
		for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
		{
			double sumv;
			sumv = 0.;
			pops_per_vib[0][iVibHi] = 0.;

			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{
				pops_per_elec[0] += H2_populations[iElecHi][iVibHi][iRotHi];
				sumv += H2_populations[iElecHi][iVibHi][iRotHi];
				pops_per_vib[0][iVibHi] += H2_populations[iElecHi][iVibHi][iRotHi];
			}
			/* print sum of H2_populations in each vib if trace on */
			if(mole.lgH2_TRACE) fprintf(ioQQQ,"\t%.2e",sumv/hmi.H2_total);
		}
		/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
		if(mole.lgH2_TRACE) 
		{
			fprintf(ioQQQ,"\n");
			/* print the ground vib state */
			fprintf(ioQQQ," Pop(0,J):");
			iElecHi = 0;
			iVibHi = 0;
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{
				fprintf(ioQQQ,"\t%.2e",H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total);
			}
			fprintf(ioQQQ,"\n");
		}

		/* now find population in states done with matrix - this is only used to pass
		 * to matrix solver */
		iElec = 0;
		sum_pops_matrix = 0.;
		ip =0;
		for( i=0; i<nXLevelsMatrix; ++i )
		{   
			ip = H2_ipX_ener_sort[i];
			iVib = ipVib_H2_energy_sort[ip];
			iRot = ipRot_H2_energy_sort[ip];
			sum_pops_matrix += H2_populations[iElec][iVib][iRot];
		}
		/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
		/* this is self consistent since pops_per_elec[0] came from current soln,
		* as did the matrix.  pops will be renomalized by results from the chemistry
		* a few lines down */
		frac_matrix = sum_pops_matrix / pops_per_elec[0];

		/* assuming that all H2 population is in X, this is the
		 * ratio of H2 that came out of the chemistry network to what we just obtained -
		 * we need to multiply the pops by renorm to agree with the chemistry,
		 * this routine does not alter hmi.H2_total, but does change pops_per_elec */
		H2_renorm_conserve = hmi.H2_total/pops_per_elec[0];
		pops_per_elec[0] = hmi.H2_total;

		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				iElec = 0;
				iVib = 0;
				fprintf(ioQQQ,"DEBUG loop pops");
				for(iRot=0; iRot<3; ++iRot )
				{
					fprintf(ioQQQ,"\t%.3e\t%.3e\t%.3e\t%.3e", 
						H2_populations[iElec][iVib][iRot]/hmi.H2_total,

						(H2_col_rate_in[iVib][iRot]+H2_rad_rate_in[iVib][iRot])/
						SDIV(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[iVib][iRot])/hmi.H2_total,

						(H2_col_rate_in[iVib][iRot]+H2_rad_rate_in[iVib][iRot]),
						SDIV(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[iVib][iRot]) );
				}
				fprintf(ioQQQ,"\n");
			}
		}

		/* renormalize H2_populations  - the H2_populations were updated by renorm when routine entered, 
		 * before pops determined - but population determinations above do not have a sum rule on total
		 * population - this renorm is to preserve toal population */
		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 )
				{
					H2_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve;
					/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
				}
			}
		}

		/* now save this solution to use for next old pop,
		 * find largest change in old and new populations */
		sumold = 0.;
		for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
		{
			for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
			{
				for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
				{

					/* keep track of largest relative change in H2_populations to
					 * determince convergence */
					if( fabs(H2_populations[iElec][iVib][iRot] - 
						H2_old_populations[iElec][iVib][iRot])/
						/* on first call some very high J states can have zero pop */
						SDIV(H2_populations[iElec][iVib][iRot]) > fabs(PopChgMax) &&
						/* >>chng 03 jul 19, this had simply been H2_populations > SMALLFLOAT,
						* change to relative pops > 1e-15, spent too much time converging
						* levels at pops = 1e-37 */
						/* >>chng 03 dec 27, from rel pop 1e-15 to 1e-6 since converging heating will
						* be main convergence criteria */
						/*H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(hmi.H2_total)>1e-15 )*/
						H2_populations[iElec][iVib][iRot]/SDIV(hmi.H2_total)>1e-6 )
					{
						PopChgMax = 
							(H2_populations[iElec][iVib][iRot] - 
							H2_old_populations[iElec][iVib][iRot])/
							SDIV(H2_populations[iElec][iVib][iRot]);
						iRotMaxChng = iRot;
						iVibMaxChng = iVib;
					}


					kase = -1;
					/* update the populations - we used the old populations to update the
					 * current new populations - will do another iteration if they changed
					 * by much.  here old populations are updated for next sweep through molecule */
					/* pop oscillatinos have occurred - use small changes */
					/* >>chng 04 may 10, turn this back on - now with min on how small frac new
					 * can become */
					if( lgH2_pops_ever_oscil )
					{
						/* loscillations - take weighted average */
						/* frac_new_oscil was set to 0.25 upon entry to this routine,
						 * was set smaller each time an oscillation occurs */
						H2_old_populations[iElec][iVib][iRot] = 
							(1.f-frac_new_oscil)*H2_old_populations[iElec][iVib][iRot] +
							frac_new_oscil*H2_populations[iElec][iVib][iRot];
						kase = 1;
					}

					/* this branch very large changes, use mean of logs */
					else if( fabs( H2_old_populations[iElec][iVib][iRot] - 
						H2_populations[iElec][iVib][iRot] )/
						SDIV( H2_populations[iElec][iVib][iRot] ) > 3. &&
						H2_old_populations[iElec][iVib][iRot]*H2_populations[iElec][iVib][iRot]!=0  )
					{
						/* large changes or oscillations - take average in the log */
						H2_old_populations[iElec][iVib][iRot] = (float)pow( 10. , 
							log10(H2_old_populations[iElec][iVib][iRot])/2. +
							log10(H2_populations[iElec][iVib][iRot])/2. );
						kase = 2;
					}

					/* modest change, use means of old and new */
					else if( fabs( H2_old_populations[iElec][iVib][iRot] - 
						H2_populations[iElec][iVib][iRot] )/
						SDIV( H2_populations[iElec][iVib][iRot] ) > 0.5  )
					{
						/* large changes or oscillations - take average */
						H2_old_populations[iElec][iVib][iRot] = 
							FRAC_OLD*H2_old_populations[iElec][iVib][iRot] +
							(1.f-FRAC_OLD)*H2_populations[iElec][iVib][iRot];
						kase = 3;
					}
					else
					{
						/* small changes, use new value */
						H2_old_populations[iElec][iVib][iRot] = 
							H2_populations[iElec][iVib][iRot];
						kase = 4;
					}
					sumold += H2_old_populations[iElec][iVib][iRot];
				}
			}
		}
		/* will renormalize so that total population is correct */
		H2_renorm_conserve_init = hmi.H2_total/sumold;

		/* renormalize H2_populations  - the H2_populations were updated by renorm when routine entered, 
		 * before pops determined - but population determinations above do not have a sum rule on total
		 * population - this renorm is to preserve toal population */
		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 )
				{
					H2_old_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve_init;
					/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
				}
			}
		}
		/* get current ortho-para ratio, will be used as test on convergence */
		iElecHi = 0;
		h2.ortho_density = 0.;
		h2.para_density = 0.;
		for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
		{
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{
				if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] )
				{
					h2.ortho_density += H2_populations[iElecHi][iVibHi][iRotHi];
				}
				else
				{
					h2.para_density += H2_populations[iElecHi][iVibHi][iRotHi];
				}
				/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
			}
		}
		/* these will be used to determine whether soln has converged */
		ortho_para_older = ortho_para_old;
		ortho_para_old = ortho_para_current;
		ortho_para_current = h2.ortho_density / SDIV( h2.para_density );

		/* >>chng 04 mar 12, move creation of solomon process up to here from
		 * after convergence loop since will check on value for convergence */
		/* this is total rate (s-1) h2 dissoc into X continuum by Solomon process,
		* sum over all excited elec states */
		oldrate = hmi.H2_Solomon_dissoc_rate_BigH2_total;
		hmi.H2_Solomon_dissoc_rate_BigH2_total = 0.;

		for( iElec=1; iElec<mole.n_h2_elec_states; ++iElec )
		{
			for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
			{
				for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
				{
					/* this is the total rate of dissociation of excited elec states into 
					* the X continuum.  At this point the units are cm-3 s-1 */
					/* >>chng 04 mar 13, use old pops sinc ehave been renorm */
					hmi.H2_Solomon_dissoc_rate_BigH2_total += 
						/*H2_populations[iElec][iVib][iRot]*H2_dissprob[iElec][iVib][iRot];*/
						H2_old_populations[iElec][iVib][iRot]*H2_dissprob[iElec][iVib][iRot];
				}
			}
		}
		/* at this point units of hmi.H2_Solomon_dissoc_rate_BigH2_total are cm-3 s-1 
		 * since H2_populations are included -
		 * div by pops to get actual dissocation rate, s-1 */
		if( hmi.H2_total > SMALLFLOAT )
		{
			hmi.H2_Solomon_dissoc_rate_BigH2_total /= hmi.H2_total;
			/* will be used for H2s-> H + H */
			hmi.H2_Solomon_dissoc_rate_BigH2_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_total* ratio_H2s_pump;
			/* will be used for H2g-> H + H */
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_total* (1.- ratio_H2s_pump );
		}
		else
		{
			hmi.H2_Solomon_dissoc_rate_BigH2_total = 0.;
			hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0; 
			hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0;		
		}

		/* are changes too large? must decide whether population shave converged,
		 * will check whether H2_populations themselves have changed by much,
		 * but also change in heating by collisional deexcitation is stable */
		HeatChangeOld = HeatChange;
		HeatChange = old_HeatH2Dexc_BigH2 - hmi.HeatH2Dexc_BigH2;
		/* check whether pops are oscillating, as evidenced by change in
		 * heating changing sign */
		if( loop_h2_pops>2 && (
			(HeatChangeOld*HeatChange<0. ) ||
			(PopChgMax*PopChgMaxOld<0. ) ) )
		{
			lgH2_pops_oscil = TRUE;
			if( loop_h2_pops > 6 )
			{
				lgH2_pops_ever_oscil = TRUE;
				/* make this smaller in attempt to damp out oscillations,
				 * but don't let get too small*/
				frac_new_oscil *= 0.8f;
				frac_new_oscil = MAX2( frac_new_oscil , 0.1f);
				++n_pop_oscil;
			}
		}
		else
		{
			lgH2_pops_oscil = FALSE;
		}

		/* reevaluate heating - cooling if H2 molecule is significant source or either,
		 * since must have stable heating cooling rate */
		old_HeatH2Dexc_BigH2 = hmi.HeatH2Dexc_BigH2;
		if(fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > conv.HeatCoolRelErrorAllowed/10. ||
			hmi.HeatH2Dexc_BigH2==0. )
			H2_Cooling();

		/* begin check on whether soln is converged */
		lgConv_h2_soln = TRUE;
		lgPopsConv = TRUE;;
		lgHeatConv = TRUE;
		lgSolomonConv = TRUE;
		lgOrthoParaRatioConv = TRUE;

		/*fprintf(ioQQQ,"DEBUG orthopara\t%.4e\t%.4e\t%.4e\t%.4e\t%c\n",
			ortho_para_current,
			ortho_para_old,ortho_para_older,
			fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current),
			TorF((ortho_para_current-ortho_para_old)*(ortho_para_old-ortho_para_older)>0. ));*/
		/* these are all the convergence tests */
		if( fabs(PopChgMax)>0.1)
		{
			/*lgPopsConv = (fabs(PopChgMax)<=0.1);*/
			lgConv_h2_soln = FALSE;
			lgPopsConv = FALSE;
			/* >>chng 04 sep 08, set quant_new to new chng max gs */
            /*quant_old = PopChgMax;*/
			quant_old = PopChgMaxOld;
			/*quant_new = 0.;*/
            quant_new = PopChgMax;

			strcpy( chReason , "populations changed" );
		}

		/* >>chng 04 apr 30, look at change in ortho-para ratio, also that is not
		 * oscillating */
		else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> 1e-3 &&
			(ortho_para_current-ortho_para_old)*(ortho_para_old-ortho_para_older)>0. )
		{
			lgConv_h2_soln = FALSE;
			lgOrthoParaRatioConv = FALSE;
			quant_old = ortho_para_old;
			quant_new = ortho_para_current;
			strcpy( chReason , "ortho/para ratio changed" );
		}

		else if( fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot > 
			conv.HeatCoolRelErrorAllowed/5.
			/* >>chng 04 may 09, do not check on error in heating if constant temperature */
			/*&& !(thermal.lgTSetOn || phycon.te <= StopCalc.TeLowest  )*/ )
		{
			/* default on HeatCoolRelErrorAllowed is 0.02 */
			/*lgHeatConv = (fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot <=
			 * conv.HeatCoolRelErrorAllowed/5.);*/
			lgConv_h2_soln = FALSE;
			lgHeatConv = FALSE;
			quant_old = old_HeatH2Dexc_BigH2/thermal.ctot;
			quant_new = hmi.HeatH2Dexc_BigH2/thermal.ctot;
			strcpy( chReason , "heating changed" );
			/*fprintf(ioQQQ,"DEBUG old new trip \t%.4e \t %.4e\n",
				old_HeatH2Dexc_BigH2,
				hmi.HeatH2Dexc_BigH2);*/
		}

		/* check on Solomon rate,
		 * >>chng 04 aug 28, do not do this check if induced processes are disabled,
		 * since Solomon process is then irrelevant */
		/*else if( rfield.lgInducProcess &&
			fabs(oldrate - hmi.H2_Solomon_dissoc_rate_BigH2_total)/SDIV(hmi.H2_rate_destroy) > 
			conv.EdenErrorAllowed/5.)*/
		/* >>chng 04 sep 21, GS*/
         else if( rfield.lgInducProcess &&
			fabs( hmi.H2_Solomon_dissoc_rate_BigH2_total - oldrate)/SDIV(hmi.H2_Solomon_dissoc_rate_BigH2_total) > 
			conv.EdenErrorAllowed/5.)
		{
			/*lgSolomonConv = (fabs(oldrate - hmi.H2_Solomon_dissoc_rate_BigH2_total)/SDIV(hmi.H2_rate_destroy) <= 
			conv.EdenErrorAllowed/5.);*/
			lgConv_h2_soln = FALSE;
			lgSolomonConv = FALSE;
			quant_old = oldrate;
			quant_new = hmi.H2_Solomon_dissoc_rate_BigH2_total;
			strcpy( chReason , "Solomon rate changed" );
		}

		/* did we pass all the convergence test */
		if( !lgConv_h2_soln )
			/*(fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > 0.1 && fabs(PopChgMax)>0.03) )*/
			/* >>chng 03 dec 27, use actual deexcitation heating to determine whether pops are OK */
			/* >>chng 03 dec 28, only converge down to half the allowed error - will take the mean of the
			 * old and new heating rates - the numerical scheme cannot get high precision, would be better
			 * to go over to matrix inversion */
			 /* >>chng 04 jan 15, from heat error /2 to 100% of allowed, soln does start
			  * to oscillate when error is ~1%, but does get soln */
			 /* >>chng 04 feb 24, x-ray pdr had trouble converging - bring criteria
			  * back down to error / 5 */
			/*!lgHeatConv ||
			!lgSolomonConv)*/
		{
			/* this branch H2 H2_populations within X are not converged,
			 * print diagnostic */
			if( PRT_POPS || mole.lgH2_TRACE )
			{
				/*fprintf(ioQQQ,"temppp\tnew\t%.4e\tnew\t%.4e\t%.4e\n",
					hmi.HeatH2Dexc_BigH2,
					old_HeatH2Dexc_BigH2,
					fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot );*/
				fprintf(ioQQQ,"DEBUG H2 pops no conv, reason:%s ",chReason );
				if( !lgPopsConv )
					fprintf(ioQQQ," pops \t%2e",PopChgMax);
				else if( !lgHeatConv )
					fprintf(ioQQQ," heat \t%2e",(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot);
				else if( !lgSolomonConv )
					fprintf(ioQQQ," dest \t%2e",(oldrate - hmi.H2_Solomon_dissoc_rate_BigH2_total)/SDIV(hmi.H2_rate_destroy));
				else if( !lgOrthoParaRatioConv )
					fprintf(ioQQQ," current, old, older ratios are %.4e %.4e %.4e",
					ortho_para_current , ortho_para_old, ortho_para_older );
				else
					TotalInsanity();
				fprintf(ioQQQ,
					"\t%.2f\t%li\t%i\t%.2e\t PopChgMax\t%.4e\t renorm chem\t%.4e\t renorm cons\t%.4e\tJ\t%li\tv\t%li\tHoscil?\t%c\tsol\t%.3e\t%.3e\n",
					fnzone , 
					loop_h2_pops , 
					kase,
					fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot ,
					PopChgMax ,
					H2_renorm_chemistry , /*renorm chem */
					H2_renorm_conserve , /* renorm cons */
					iRotMaxChng ,
					iVibMaxChng,
					TorF(lgH2_pops_oscil) ,
					hmi.H2_Solomon_dissoc_rate_BigH2_total,
					hmi.H2_total);
			}
		}
		/* end convergence criteria */

		/*fprintf(ioQQQ,"DEBUG h2 heat\t%3li\t%.2f\t%.4e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
			loop_h2_pops,
			fnzone,
			phycon.te,
			dense.eden,
			hmi.HeatH2Dexc_BigH2,
			hmi.HeatH2Dexc_BigH2/thermal.ctot ,
			hmi.H2_total,
			H2_renorm_chemistry ,
			H2_renorm_conserve,
			hmi.H2_H2g_to_H2s_rate_BigH2);*/
		if( trace.lgTrConvg >= 5 )
		{
			fprintf( ioQQQ, 
				"     H2 5lev %li Conv?%c",
				loop_h2_pops ,
				TorF(lgConv_h2_soln) );

			if( fabs(PopChgMax)>0.1 )
				fprintf(ioQQQ," pops, rel chng %.3e",PopChgMax);
			else
				fprintf(ioQQQ," rel heat %.3e rel chng %.3e H2 heat/cool %.2e",
					hmi.HeatH2Dexc_BigH2/thermal.ctot ,
					fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot ,
					hmi.HeatH2Dexc_BigH2/thermal.ctot);

			fprintf( ioQQQ, 
				" Oscil?%c Ever Oscil?%c",
				TorF(lgH2_pops_oscil) ,
				TorF(lgH2_pops_ever_oscil) );
			if( lgH2_pops_ever_oscil )
				fprintf(ioQQQ," frac_new_oscil %.4f",frac_new_oscil);
			fprintf(ioQQQ,"\n");
		}

		if( mole.lgH2_TRACE) 
		{
			fprintf(ioQQQ,
			"H2 loop\t%li\tkase pop chng\t%i\tchem renorm fac\t%.4e\tortho/para ratio:\t%.3e\tfrac of pop in matrix: %.3f\n",
			loop_h2_pops,
			kase,
			H2_renorm_chemistry,
			h2.ortho_density / h2.para_density ,
			frac_matrix);

			/* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
			if( iVibMaxChng>=0 && iRotMaxChng>=0 )
				fprintf(ioQQQ,
					"end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n",
					loop_h2_pops,
					PopChgMax , 
					H2_old_populations[0][iVibMaxChng][iRotMaxChng],
					H2_populations[0][iVibMaxChng][iRotMaxChng],
					iVibMaxChng , iRotMaxChng
					);
		}
	}
	/* =======================END POPULATIONS CONVERGE LOOP =====================*/

	if( !lgConv_h2_soln )
	{
		conv.lgConvPops = FALSE;
		strcpy( conv.chConvIoniz, "H2 pop cnv" );
		fprintf(ioQQQ,
			"  H2_LevelPops:  H2_populations not converged in %li tries; due to %s, old, new are %.4e %.4e, iteration %li zone %.2f.\n",
			loop_h2_pops, 
			chReason,
			quant_old,
			quant_new ,
			iteration , 
			fnzone );
		ConvFail("pops","H2");
	}

	/* loop over all possible lines and set H2_populations, 
	 * and quantities that depend on them */
	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 )
			{
				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 )
						{
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi = 
								H2_populations[iElecHi][iVibHi][iRotHi]; 
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo = 
								H2_populations[iElecLo][iVibLo][iRotLo]; 
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc = 
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo - 
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi*
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo / 
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi;
							ASSERT(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi >= 0. &&
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo >=  0.);

							/* >>chng 03 feb 14, from !=0 to >0 */
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								/* following two heat exchange excitation, deexcitation */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cool = 0.;
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].heat = 0.;

								/* number of photons in the line */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots = 
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul * 
									(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc + 
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc) * 
									 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi; 

								/* intensity of line */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity = 
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots *
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg;

								if( iElecHi==0 )
								{
									/*TODO	2	- put H2Lines in outward beams in RT_diffuse */
									/* the ground electronic state, most excitations are not direct pumping 
									 * (rather indirect, which does not count for ColOvTot) */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot = 1.;
								}
								else
								{
									/* these are excited electronic states, mostly pumped, except for supras */
									/*TODO	2	put supra thermal excitation into excitation of electronic bands */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot = 0.;
								}
							}
						}
					}
				}
			}
		}
	}	

	/* solomon rate will be updated below, then will take mean of old and new - save
	 * current rate */
	oldrate = hmi.H2_Solomon_dissoc_rate_BigH2_total;
#	if 0
	/* more to within convergence loop since will test on stability of this rate */
	/* this is total rate (s-1) h2 dissoc into X continuum by Solomon process,
	 * sum over all excited elec states */
	hmi.H2_Solomon_dissoc_rate_BigH2_total = 0.;
	for( iElec=1; iElec<mole.n_h2_elec_states; ++iElec )
	{
		for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
		{
			for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
			{
				/* this is the total rate of dissociation of excited elec states into 
				 * the X continuum.  At this point the units are cm-3 s-1 */
				hmi.H2_Solomon_dissoc_rate_BigH2_total += 
					H2_populations[iElec][iVib][iRot]*H2_dissprob[iElec][iVib][iRot];
			}
		}
	}
	if( hmi.H2_total > SMALLFLOAT )
	{
		hmi.H2_Solomon_dissoc_rate_BigH2_total /= hmi.H2_total;
	}
#	endif

	/* add up H2 + hnu => 2H, continuum photodissociation,
	 * this is not the Solomon process, true continuum */
	hmi.H2_photodissoc_BigH2 = 0.;
	iElec = 0;
	hmi.H2_star_BigH2 = 0.;
	for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
	{
		for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
		{
			
			/* this is the total rate of direct photo-dissociation of excited elec states into 
			 * the X continuum - this is continuum photodissociation, not the Solomon process */
			hmi.H2_photodissoc_BigH2 += 
				H2_populations[iElec][iVib][iRot]/SDIV(hmi.H2_total)*rfield.flux_accum[H2_ipPhoto[iVib][iRot]-1]*0.25e-18f;

			/* >>chng 03 sep 03, make sum of pops of excited states */
			if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
			{
				/* sum of pops in levels with iVib > 0 for use in chemistry network */
				hmi.H2_star_BigH2 += H2_populations[iElec][iVib][iRot];
			}
		}
	}
	if( mole.lgH2_TRACE || (trace.lgTrace && trace.lgTr_H2_Mole) )
	{
		fprintf(ioQQQ," H2_LevelPops exits finding Sol dissoc %.2e (TH85 %.2e)",
			hmi.H2_Solomon_dissoc_rate_BigH2_total , 
			hmi.H2_Solomon_dissoc_rate_TH85);

		/* Solomon process rate from X into the X continuum with units s-1
		 * rates are total rate, and rates from H2g and H2s */ 
		fprintf(ioQQQ," H2g Sol %.2e H2s Sol %.2e",
			hmi.H2_Solomon_dissoc_rate_used_H2g , 
			hmi.H2_Solomon_dissoc_rate_BigH2_H2s );

		/* photoexcitation from H2g to H2s */
		fprintf(ioQQQ," H2g->H2s %.2e (TH85 %.2e)",
			hmi.H2_H2g_to_H2s_rate_BigH2 , 
			hmi.H2_H2g_to_H2s_rate_TH85);

		/* add up H2s + hnu => 2H, continuum photodissociation,
		 * this is not the Solomon process, true continuum, units s-1 */
		fprintf(ioQQQ," H2 con diss %.2e (TH85 %.2e)\n",
			hmi.H2_photodissoc_BigH2 , 
			hmi.H2_photodissoc_TH85);

	}
	/* >>chng 03 sep 01, add this population - before had just used H2star from chem network */
	/* if big H2 molecule is turned on and used for this zone, use its
	 * value of H2* (pops of all states with v > 0 ) rather than simple network */
	/*fprintf(ioQQQ," debuggg h2star\t%.2f\t%.3e\t%.3e\n",
		fnzone , hmi.Hmolec[ipMH2s] , hmi.H2_star_BigH2 );*/

	/* update number of times we have been called */
	++nCallH2_this_iteration;

	/* this will say how many times the large H2 molecule has been called in this zone -
	 * if not called (due to low H2 abundance) then not need to update its line arrays */
	++h2.nCallH2_this_zone;

	return;
}

/* 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 */
#	if 0
#	define N_H2_LINES_OUTPUT	1
	long int ip_H2_lines_output[N_H2_LINES_OUTPUT*6] = {
	/* R branch - J decreases by 1 in emission line */
	/* P branch - J increases by 1 in emission line */
	/* Q branch - J changes   by 0 in emission line */
	/* S branch - J decreases by 2 in emission line */
	/*0,0,2, 0,0,0 ,  S(0) */
	/* 0,0,3, 0,0,1 , S(1) */
	/*0,0,4, 0,0,2 ,  S(2) */
	/*0,0,5, 0,0,3 ,  S(3) */
	/*0,0,6, 0,0,4 ,  S(4) */
	0,1,3, 0,0,1 } /* 1-0 S(1) */;

	/* the most commonly observed H2 lines:
	 * 0-0 R(5),
	 * 0-0 S(0), 
	 * 0-0 S(1), 
	 * 0-0 S(2), 
	 * 0-0 S(3), 
	 * 0-0 S(5), 
	 * 0-0 S(19), 
	 * 0-0 S(25), 
	 * 1-0 Q(3), 
	 * 1-0 Q(5), 
	 * 1-0 R(5) 
	 * 1-0 S(1), 
	 * 1-0 S(7), 
	 * 1-0 Q(3), 
	 * 2-1 S(1), 
	 * 6-4 O(3),  
	 */
#	endif
	int iRotHi, iVibHi, iElecHi ,iRotLo, iVibLo, iElecLo,
		ipHi , ipLo;

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

	iElecHi = 0;
	iElecLo = 0;
#	if 0
	/* print all the v = 0, 1, ro-vib lines */
	for( iVibHi=0; iVibHi<2; ++iVibHi )
	{
		for( iRotHi=2; iRotHi<nRot_hi[0][0]; ++iRotHi )
		{
			/* all ground vib state rotation lines - first is J to J-2 */
			PutLine(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotHi-2]);
			/* next is J to J - may not have been defined for some rot,
			 * and certainly not defined for iVibHi = 0 - since J=J and v=v */
			if( iVibHi > 0 && H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotHi].Aul > 0. )
				PutLine(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotHi]);
		}
	}
#	endif
	/* print all lines from lowest n levels within X */
	for( ipHi=2; ipHi<nLevels_per_elec[iElecHi]; ++ipHi )
	{
		/* obtain the proper indices for the upper level */
		long int ip = H2_ipX_ener_sort[ipHi];
		iVibHi = ipVib_H2_energy_sort[ip];
		iRotHi = ipRot_H2_energy_sort[ip];
		for( ipLo=0; ipLo<ipHi; ++ipLo )
		{
			ip = H2_ipX_ener_sort[ipLo];
			iVibLo = ipVib_H2_energy_sort[ip];
			iRotLo = ipRot_H2_energy_sort[ip];
			if( iVibHi >= iVibLo && (abs(iRotLo-iRotHi)==2 || (iRotLo-iRotHi)==0) )
			{
				if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul>0. )
				{
					/* 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[iVibHi][iRotHi][iVibLo][iRotLo] =  0.;
					}
					else if( LineSave.ipass == 1 )
					{
						H2_SaveLine[iVibHi][iRotHi][iVibLo][iRotLo] += (float)(
							radius.dVeff*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity);
					}
				}
			}
		}
	}
	return;
#	if 0
	for( i=0; i<N_H2_LINES_OUTPUT; ++i )
	{
		/* the lower state */
		iElecLo = ip_H2_lines_output[i*6+3];
		iVibLo = ip_H2_lines_output[i*6+4];
		iRotLo = ip_H2_lines_output[i*6+5];
		/* the upper state */
		iElecHi = ip_H2_lines_output[i*6];
		iRotHi = ip_H2_lines_output[i*6+2];
		iVibHi = ip_H2_lines_output[i*6+1];

		/* dump the line */
		PutLine(&H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]);
	}

	/* only do this if molecule has been updated in this zone */
	/* this saves all lines for possible punch output */
	iElecLo = 0;
	iElecHi = 0;
	/* first branch is when called after space allocated, but before intensities
	 * known, or first zone computed - job is to zero out storage arrays */
	if( LineSave.ipass == 0 )
	{
		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 )
					{
						H2_SaveLine[iVibHi][iRotHi][iVibLo][iRotLo] =  0.;
					}
				}
			}
		}
	}
	else if( LineSave.ipass == 1 )
	{
		/* only add the flux into the save array if the H2 molecule was updated in this zone */
		if( h2.nCallH2_this_zone )
		{
			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( (iElecHi==iElecLo) && (iVibHi==iVibLo) && (iRotHi==iRotLo) )
								continue;
							ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul >= 0. );
							if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )
							{
								ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity>= 0.);
								ASSERT( H2_SaveLine[iVibHi][iRotHi][iVibLo][iRotLo] >= 0.);
								H2_SaveLine[iVibHi][iRotHi][iVibLo][iRotLo] += (float)(
									radius.dVeff*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity);
							}
						}
					}
				}
			}
		}
	}
#	endif
}
#if 0
/* this was attempt at lepp & shull coll diss rate, but formalism can't work */
/*H2_H_CollDissoc H2 collisional dissociation by H */
static double H2_H_CollDissoc( long iVib , long iRot )
{
	double rate;
	static float tused=-1;
	static double rate00=-1;
	/* collisional dissociation rates taken from 
	 * >>refer	H2	coll diss	Lepp, S., & Shull, J.M., 1983, ApJ, 270, 578-582 */

	if( tused == phycon.te )
	{
		tused = phycon.te;
		/* this is the low density rate, k_D(0,0) in Lepp Shull equation 4 */
		if( phycon.te > 7390. )
		{
			rate00 = 6.11e-14 * sexp( 2.93e4/phycon.te);
		}
		else
		{
			rate00 = 2.67e-15 * sexp( POW2(6750./phycon.te) );
		}
	}
	/* find rate for iVib and iRot from above, and eqn 4 of LS */
	rate = rate00 * exp( 13.45

	return rate;
}
#endif

/*H2_cooling evaluate cooling and heating due to H2 molecule, called by H2_LevelPops
 * in convergence loop when h2 heating is important, also called by coolr to get
 * final heating */
void H2_Cooling(void)
{
	long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
	double rate_dn_heat, 
		rate_up_cool;
	long int nColl,
		ipHi, ipLo;
	double Big1;
	long int ipVib_big_hi,ipVib_big_lo ,ipRot_big_hi ,ipRot_big_lo;
	static float old_HeatH2Dexc=0.;

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

	old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;

	/* nCallH2_this_iteration is not incremented until after the level
	 * populations have converged the first time.  so for the first n calls
	 * this will return zero, a good idea since populations will be wildly
	 * incorrect during search for first valid pops */
	if( !h2.lgH2ON || !nCallH2_this_iteration )
	{
		hmi.HeatH2Dexc_BigH2 = 0.;
		hmi.HeatH2Dish_BigH2 = 0.;
		hmi.deriv_HeatH2Dexc_BigH2 = 0.;
		return;
	}

	hmi.HeatH2Dish_BigH2 = 0.;
	/* heating due to dissociation of elec excited states */
	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 )
			{
				hmi.HeatH2Dish_BigH2 += 
					H2_populations[iElecHi][iVibHi][iRotHi] * 
					H2_dissprob[iElecHi][iVibHi][iRotHi] *
					H2_disske[iElecHi][iVibHi][iRotHi];
			}
		}
	}
	/* dissociation heating HeatH2Dish_BigH2 was in eV - 
	 * convert to ergs */
	hmi.HeatH2Dish_BigH2 *= (float)EN1EV;

	/* now work on collisional heating due to bound-bound
	 * collisional transitions within X */
	hmi.HeatH2Dexc_BigH2 = 0.;
	/* these are the colliders that will be considered as depopulating agents */
	/* the colliders are H, He, H2 ortho, H2 para, H+ */
	/* atomic hydrogen */
#	if 0
	collider_density[0] = dense.xIonDense[ipHYDROGEN][0];
	/* atomic helium */
	collider_density[1] = dense.xIonDense[ipHELIUM][0];
	/* all ortho h2 */
	collider_density[2] = (float)h2.ortho_density;
	/* all para H2 */
	collider_density[3] = (float)h2.para_density;
	/* protons - ionized hydrogen */
	collider_density[4] = dense.xIonDense[ipHYDROGEN][1];
#	endif

	/* now make sum of all collisions within X itself */
	iElecHi = 0;
	iElecLo = 0;
	Big1 = 0.;
	ipVib_big_hi = 0;
	ipVib_big_lo = 0;
	ipRot_big_hi = 0;
	ipRot_big_lo = 0;
	/* this will be derivative */
	hmi.deriv_HeatH2Dexc_BigH2 = 0.;
	for( ipHi=1; ipHi<nLevels_per_elec[iElecHi]; ++ipHi )
	{
		long int ip = H2_ipX_ener_sort[ipHi];
		iVibHi = ipVib_H2_energy_sort[ip];
		iRotHi = ipRot_H2_energy_sort[ip];
		if( iVibHi > VIB_COLLID )
			continue;

		for( ipLo=0; ipLo<ipHi; ++ipLo )
		{
			double oneline;
			ip = H2_ipX_ener_sort[ipLo];
			iVibLo = ipVib_H2_energy_sort[ip];
			iRotLo = ipRot_H2_energy_sort[ip];
			if( iVibLo > VIB_COLLID)
				continue;

			rate_dn_heat = 0.;
			rate_up_cool = 0.;

			/* this sum is total downward heating */
			for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
			{
				rate_dn_heat +=
					H2_populations[iElecHi][iVibHi][iRotHi] * 
					H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
					collider_density[nColl];

				/* now get upward collisional cooling by detailed balance */
				rate_up_cool += 
					H2_populations[iElecLo][iVibLo][iRotLo] *
					/* downward collision rate */
					H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
					collider_density[nColl]*
					/* rest converts into upward collision rate */
					H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElecLo][iVibLo][iRotLo] *
					H2_Boltzmann[iElecHi][iVibHi][iRotHi] /
					SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] );
			}/* end loop over all colliders */

			/* net heating due to collisions within X - this will usually be heating -
			 * in printout and in punch heating this is called "H2cX" */
			oneline = (rate_dn_heat - rate_up_cool)*
				(energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]) *
				ERG1CM;
			hmi.HeatH2Dexc_BigH2 += (float)oneline;

			/* deriv wrt temperature - assume exp wrt ground - this needs to be
			 * divided by square of temperature in wn - 
			 * done at end of loop */
			hmi.deriv_HeatH2Dexc_BigH2 +=  (float)(oneline *energy_wn[iElecHi][iVibHi][iRotHi]);

			/* option to keep track of strongest single heating agent */
#			define H2COOLDEBUG	FALSE
			if(H2COOLDEBUG)
			{
				if( fabs(oneline) > fabs(Big1 ) )
				{
					Big1 = oneline;
					ipVib_big_hi = iVibHi;
					ipVib_big_lo = iVibLo;
					ipRot_big_hi = iRotHi;
					ipRot_big_lo = iRotLo;
				}
			}

			/* this would be a major logical error */
			ASSERT( 
				(rate_up_cool==0 && rate_dn_heat==0) || 
				(energy_wn[iElecHi][iVibHi][iRotHi] > energy_wn[iElecLo][iVibLo][iRotLo]) );
		}/* end loop over lower levels, all collisions within X */
	}/* end loop over upper levels, all collisions within X */
	if(H2COOLDEBUG)
		fprintf(ioQQQ," H2_Cooling, total coll %.2e, frac 1 line %.2e %li %li %li %li \n",
			hmi.HeatH2Dexc_BigH2 , Big1/hmi.HeatH2Dexc_BigH2 ,
			ipVib_big_hi , ipRot_big_hi , ipVib_big_lo , ipRot_big_lo );

	/* this is inside h2 cooling, and is called extra times when H2 heating is important */
	if( PRT_POPS ) 
		fprintf(ioQQQ,
		"  DEBUG H2 heat fnzone\t%.2f\trenrom\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n",
		fnzone , 
		H2_renorm_chemistry , 
		phycon.te , 
		hmi.HeatH2Dexc_BigH2,
		hmi.HeatH2Dexc_BigH2/thermal.ctot);

	/* this is deriv of collisional heating wrt temperature - needs to be
	 * divided by square of temperature in wn */
	hmi.deriv_HeatH2Dexc_BigH2 /=  POW2(phycon.te_wn) ;
	/* no way we can get zero for this, unless collisions have been turned off,
	 * except for underflow */
	ASSERT( hmi.deriv_HeatH2Dexc_BigH2!=0. || !mole.lgColl_deexec_Calc );

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (fabs(hmi.HeatH2Dexc_BigH2) > SMALLFLOAT) )
		{
			int iVib = 0;

			/*fprintf(ioQQQ," H2_cooling pops\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
				H2_populations[0][iVib][0]/hmi.H2_total,
				H2_populations[0][iVib][1]/hmi.H2_total,
				H2_populations[0][iVib][2]/hmi.H2_total,
				H2_populations[0][iVib][3]/hmi.H2_total,
				H2_populations[0][iVib][4]/hmi.H2_total,
				H2_populations[0][iVib][5]/hmi.H2_total);*/

			iElecHi = iElecLo = 0;
			iVibHi = iVibLo = 0;
			iRotHi = 3;
			iRotLo = 1;
			rate_dn_heat = rate_up_cool = 0.;
			/* this sum is total downward heating */
			for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
			{
				rate_dn_heat +=
					H2_populations[iElecHi][iVibHi][iRotHi] * 
					H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
					collider_density[nColl];

				/* now get upward collisional cooling by detailed balance */
				rate_up_cool += 
					H2_populations[iElecLo][iVibLo][iRotLo] *
					/* downward collision rate */
					H2_CollRate[nColl][iVibHi][iRotHi][iVibLo][iRotLo]*
					collider_density[nColl]*
					/* rest converts into upward collision rate */
					H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElecLo][iVibLo][iRotLo] *
					H2_Boltzmann[iElecHi][iVibHi][iRotHi] /
					SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] );
			}

			fprintf(ioQQQ," H2_cooling pop31\t%.3e\tdn up 31\t%.3e\t%.3e\n",
				H2_populations[0][iVib][3]/H2_populations[0][iVib][1],
				rate_dn_heat,rate_up_cool
				);
		}
	}
#	undef	H2COOLDEBUG
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC  )
		{
			static long nzdone=-1 , nzincre;
			if( nzone!=nzdone )
			{
				nzdone = nzone;
				nzincre = -1;
			}
			++nzincre;
			fprintf(ioQQQ," H2 nz\t%.2f\tnzinc\t%li\tTe\t%.4e\tH2\t%.3e\tcXH\t%.2e\tdcXH/dt%.2e\tDish\t%.2e \n",
				fnzone, 
				nzincre,
				phycon.te,
				hmi.H2_total ,
				hmi.HeatH2Dexc_BigH2,
				hmi.deriv_HeatH2Dexc_BigH2 ,
				hmi.HeatH2Dish_BigH2);

		}
	}

	/* this can be noisy due to finite accuracy of soln, so take average with
	 * previous value */
	/*>>chng 04 mar 01, do not take average */
	if( 1 || nzone <1 || old_HeatH2Dexc==0. || nCallH2_this_iteration <2)
	{
		old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;
	}
	else
	{
		hmi.HeatH2Dexc_BigH2 = (hmi.HeatH2Dexc_BigH2+old_HeatH2Dexc)/2.f;
		old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;
	}
	if( mole.lgH2_TRACE) 
		fprintf(ioQQQ,
		" H2_Cooling Ctot\t%.4e\t HeatH2Dish_BigH2 \t%.4e\t HeatH2Dexc_BigH2 \t%.4e\n" ,
		thermal.ctot , 
		hmi.HeatH2Dish_BigH2 , 
		hmi.HeatH2Dexc_BigH2 );

	return;

}


/*cdH2_colden return column density in H2, negative -1 if cannot find state,
 * header is cddrive */
double cdH2_colden( long iVib , long iRot )
{			

	/*if iVib is negative, return
	 * total column density - iRot=0
	 * ortho column density - iRot 1
	 * para column density - iRot 2 
	 * else return column density in iVib, iRot */
	if( iVib < 0 )
	{
		if( iRot==0 )
		{
			/* return total H2 column density */
			return( ortho_colden + para_colden );
		}
		else if( iRot==1 )
		{
			/* return ortho H2 column density */
			return ortho_colden;
		}
		else if( iRot==2 )
		{
			/* return para H2 column density */
			return para_colden;
		}
		else
		{
			fprintf(ioQQQ," iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n");
			return -1.;
		}
	}
	else if( h2.lgH2ON )
	{
		/* this branch want state specific column density, which can only result from
		 * evaluation of big molecule */
		int iElec = 0;
		if( iRot <0 || iVib >nVib_hi[iElec] || iRot > nRot_hi[iElec][iVib])
		{
			fprintf(ioQQQ," iVib and iRot must lie within X, returning -2.\n");
			fprintf(ioQQQ," iVib must be <= %li and iRot must be <= %li.\n",
				nVib_hi[iElec],nRot_hi[iElec][iVib]);
			return -2.;
		}
		else
		{
			return H2_X_colden[iVib][iRot];
		}
	}
	/* error condition - no valid parameter */
	else
		return -1;
}

/*H2_Colden maintain H2 column densities within X */
void H2_Colden( char *chLabel )
{
	long int iVib , iRot;

	if( !h2.lgH2ON || !h2.nCallH2_this_zone )
		return;

	if( strcmp(chLabel,"ZERO") == 0 )
	{
		/* the column density (cm-2) of ortho and para H2 */
		ortho_colden = 0.;
		para_colden = 0.;
		/* zero out formation rates and column densites */
		for( iVib = 0; iVib <= nVib_hi[0] ; ++iVib )
		{
			for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
			{
				/* space for the rotation quantum number */
				H2_X_colden[iVib][iRot] = 0.;
				H2_X_colden_LTE[iVib][iRot] = 0.;
			}
		}
	}

	else if( strcmp(chLabel,"ADD ") == 0 )
	{
		ortho_colden += h2.ortho_density*radius.drad_x_fillfac;
		para_colden += h2.para_density*radius.drad_x_fillfac;
		/*  add together column densities */
		for( iVib = 0; iVib <= nVib_hi[0] ; ++iVib )
		{
			for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
			{
				/* state specific H2 column density */
				H2_X_colden[iVib][iRot] += (float)(H2_populations[0][iVib][iRot]*radius.drad_x_fillfac);
				/* LTE state specific H2 column density - H2_populations_LTE is normed to unity
				 * so must be multiplied by total H2 density */
				H2_X_colden_LTE[iVib][iRot] += (float)(H2_populations_LTE[0][iVib][iRot]*
					hmi.H2_total*radius.drad_x_fillfac);
			}
		}
	}

	/* we will not print column densities so skip that - if not print then we have a problem */
	else if( strcmp(chLabel,"PRIN") != 0 )
	{
		fprintf( ioQQQ, " H2_Colden does not understand the label %s\n", 
		  chLabel );
		puts( "[Stop in H2_Colden]" );
		cdEXIT(EXIT_FAILURE);
	}
}

/*H2_DR choose next zone thickness based on H2 big molecule */
double H2_DR(void)
{
	return BIGFLOAT;
#	if 0
#include "geometry.h" 
#include "doppvel.h" 
	double dtau , opacity  , dr_new;
	long int iElecHi ,iVibHi ,iRotHi , iVibLo, iRotLo, iElecLo;

	iElecHi = nEkcHi_SolomonRateMax;
	iVibHi = nVibHi_SolomonRateMax;
	iRotHi = nRotHi_SolomonRateMax;
	iVibLo = nVib_SolomonRateMax;
	iRotLo = nRot_SolomonRateMax;
	iElecLo = 0;

	/* want a dr that will have dtau less than this */
	/* >>chng 03 jun 17, from 0.1 to 1 */
	dtau = MAX2(1.,H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].opacity);

	opacity = H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].opacity *geometry.FillFac/ 
		DoppVel.doppler[ H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].nelem -1]*
		H2_populations[iElecLo][iVibLo][iRotLo];

	/* above is inverse opacity, convert to scale length */
	/* the above was tau = 1 - want much smaller than this since dr selected from change in H2 */
	dr_new = dtau/ SDIV( opacity);

	if( dr_BigMole > 0. && (TRUE || dr_BigMole==radius.drad) )
	{
#		define FRAC	0.5
		/* this branch, were using big H dr to set thickness, use mean of old and new */
		fprintf(ioQQQ,"bug mean used %li dr new %.3e dTau %.3e \n", nzone,dr_new, dtau );
		dr_BigMole = dr_new*FRAC + (1.-FRAC)*dr_BigMole;
#		undef FRAC
	}
	else
	{
		/* this branch, did not use big H2 dr, either because not important or because not known */
		dr_BigMole = dr_new;
		fprintf(ioQQQ,"bug value used %li\t%e\t%e\n", nzone, dr_BigMole,radius.drad);
	}

	return dr_BigMole;
#	endif
}

/*H2_RT_OTS - add H2 ots fields */
void H2_RT_OTS( void )
{

	long int iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo;

#	ifdef DEBUG_FUN
	fputs( "<+>H2_RT_OTS()\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;

	/* loop over all possible lines and set H2_populations, and quantities that depend on escape prob, dest, etc */
	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 )
			{
				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. )
							{
								if( iElecHi==0 )
								{
									/* ots destruction rate */
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ots = (float)(
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul *
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi *
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest);

									/* dump the ots rate into the stack - but only for ground elec state*/
									RT_OTS_AddLine(
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ots,
										H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont );
								}
							}
						}
					}
				}
			}
		}
	}

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

	return;
}


