/* 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 */
/*hmole_step do a step in chemical network */
/*hmirat compute radiative association rate for H- */
/* >> chng 02 nov 7 rjrw, Mole Moreliano:
 *   changes to linearized iterative form */
/* from Robin Williams:
The process for these kind of problems seems to be pretty uniform:
switch printsol on to check which terms in the chemical matrix change,
and next to switch on the prints in the matrix assembly which apply to
the species involved to find what reactions are involved.  It's a bit
of a pain grepping down to find the 47th reaction, so I guess some
kind of naming scheme for the reactions may come in handy (I'd thought
about generating the in and out vectors from a text string, e.g. "H +
H => H2" ;-), but you'd have to verify uniqueness).
*/
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "atmdat.h"
#include "grainvar.h"
#include "dense.h"
#include "secondaries.h"
#include "co.h"
#include "opacity.h"
#include "rfield.h"
#include "mole.h"
#include "thermal.h"
#include "punch.h"
#include "timesc.h"
#include "trace.h"
#include "phycon.h"
#include "doppvel.h"
#include "lapack.h"
#include "gammas.h"
#include "hmrate.h"
#include "h2.h"
#include "dynamics.h"
#include "radius.h"
#include "hextra.h" 
#include "hmi.h"

/* HP cc cannot compile following except in -O1 mode */
#if defined(__HP_cc)
  #pragma OPT_LEVEL 1
#endif

#define ABSLIM  1e-12

#define MAXREACTANTS 3
#define MAXRATES 3
#define MAXPRODUCTS  4
/* Calculate number of elements in an integer vector */
#define INTSZ(a)     (sizeof(a)/sizeof(int))

typedef struct reaction_s {
	int index;
	double ratek;
	int nreactants, nrates, nproducts;
	int reactants[MAXREACTANTS];
	int rate_species[MAXREACTANTS];
	int products[MAXPRODUCTS];
	struct reaction_s *next;
} reaction;

/* Generate new element for reaction list */
reaction *newreaction(
	/* reaction index, a number to reference current reaction */
	int rindex, 
	/* vector of ints that are incoming species, the reactants */
	int *in, 
	/* number of reactants */
	int nin, 
	/* same for products */
	int *out, 
	/* number of products */
	int nout, 
	/* rate determining species if they are not the same as the reactants,
	 * these are non-null in cases where a part of the chemical network
	 * acts as a catalyst in the reaction 
	 * as H2g + H2g -> H2g + H2* or where */
	int *rate, 
	/* number of these */
	int nrate)
{
	static reaction *list = NULL, *r;
	static int poolsize=1, index = 0;
	int i;

	/* this routine is called only to initialize structure with information
	 * on the reactants and products - this does not change, so no need
	 * to do this but once.  In later cases only the rate coefficients are
	 * updated in hmole_step 
	 * in hmole_stop on later call linked list is incremented upon each
	 * reaction, and rate coefficient for current temperature is stored 
	 * in r->ratek 
	 * r is pointer to structure of reaction information including
	 * r->next, which is the next reaction in the linked list */
	/* fprintf(ioQQQ,"New reaction %d %d %d\n",rindex,nin,nout); */

	/* default assumption for chemical kinetics,
	 * rate is determined by all incoming species that are in *in 
	 * when not null then the rate determining species are not the
	 * same as *in */
	if (rate == NULL) 
	{
		rate=in;
		nrate=nin;
	}

	/* space for the linked list "list" */
	/*lint -e701 shift left of signed quantity */
	if (list == NULL || index == poolsize) 
	{
		poolsize <<=1;
		if( (list = ((reaction *)MALLOC( (size_t)poolsize*sizeof(reaction) ))) == NULL )
			BadMalloc();
		index = 0;
	}
	/*lint +e701 */

	/* fprintf(ioQQQ,"Getting element %d+1 of %d\n",index,poolsize); */
	r = list+index;
	index++;
	r->next = NULL;
	r->index = rindex;
	assert (nin <= MAXREACTANTS && nout <= MAXPRODUCTS && nrate <= MAXRATES);

	r->nreactants =	nin;
	r->nrates =	nrate;
	r->nproducts  =	nout;

	/* incoming reactants */
	for (i=0; i<r->nreactants; i++)
		r->reactants[i] = in[i];

	/* rate determining species */
	for (i=0; i<r->nrates; i++)
		r->rate_species[i] = rate[i];

	/* outgoing products */
	for (i=0; i<r->nproducts; i++)
		r->products[i] = out[i];

	return r;
}

void hmole_step(int *nFixup, double *error)
{
	enum {PRINTSOL = FALSE};

	/* will punch debug output to this file */
	FILE*ioFile;

	int32 ipiv[N_H_MOLEC], merror;

	long int i, 
	  ipConserve, 
	  j, 
	  limit ,
	  nd,
	  mol;

	int nprot[N_H_MOLEC],
		printsol = PRINTSOL;

	int lgNegPop,
		iworst;
	float frac_H2star_grain,
		frac_H2star_hminus; 
	double sum_first_ions;

	double assoc_detach_backwards_grnd,
		assoc_detach_backwards_exct,
	  bh2dis,
	  bh2h22hh2,
	  bhneut, 
	  c3bod,
	  cionhm, 
	  corr, 
	  H2star_deexcit,
	  deexc_htwo,
	  deexc_hneut,
	  desh2p, 
	  eh2hh, 
	  eh3_h2h,
	  etmp,
	  eh3p_3h,
	  Boltz_fac_H2_H2star,
	  fac, 
	  fhneut, 
	  gamhd, 
	  gamheh, 
	  h1fnd, 
	  h1rat, 
	  h2pcin, 
	  h2ph3p, 
	  h2phhp, 
	  h2pion, 
	  H2star_excit ,
	  h32h2, 
	  h3petc, 
	  h3ph2p, 
	  radasc, 
	  radath, 
	  rh2h2p, 
	  h2hph3p,
	  hmihph2p,
	  h2phmh2h,
	  h2phmhhh,
	  h3phmh2hh,
	  h2h22hh2,
	  h2crphh,
	  h2crphphm,
	  h2crphpeh,
	  h3phm2h2,
	  h2pehh,
	  h3ph2ph,
	  h3ph2hp,
	  heph2hpheh,
	  heph2heh2p,
	  hphmhhpe,
	  h2hmhh2e,
	  hehph2h3phe,
	  hehmeheh,
	  hephmhhe,
	  sum,
		fracneg,
		fracnegtmp,
		fracnegfac,
		sum_H0_Hp,
		rtot,
		conserve,
		rate,
		ratek,
		rated,
		rate_deriv[MAXREACTANTS],
		sinkrate[MAXREACTANTS],
		TStancil;


	static double 
	  eh2hhm, 
	  eh2old, 
	  gamtwo, 
	  h2phet, 
	  protons;

	static double /*amat[N_H_MOLEC][N_H_MOLEC], */
	  b2pcin, 
	  *amat=NULL,
	  *bvec=NULL/*[N_H_MOLEC]*/, 
	  *Hmolec_old=NULL/*[N_H_MOLEC]*/, 
	  **c=NULL/*[N_H_MOLEC+1][N_H_MOLEC+1]*/, 
	  plte;

	/* if this is still true then must create space for arrays */
	static int lgMustMalloc = TRUE;

	static reaction *rlist = NULL;
	reaction *r;
	long int rindex, ratei, ratej;

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

	if( lgMustMalloc )
	{
		/* on very first call must create space */
		lgMustMalloc = FALSE;

		if( (bvec = ((double*)MALLOC( (size_t)N_H_MOLEC*sizeof(double) ))) == NULL )
			BadMalloc();
		if( (Hmolec_old = ((double*)MALLOC( (size_t)N_H_MOLEC*sizeof(double) ))) == NULL )
			BadMalloc();
		if( (amat = ((double*)MALLOC( (size_t)(N_H_MOLEC*N_H_MOLEC)*sizeof(double) ))) == NULL )
			BadMalloc();
		if( (c = ((double**)MALLOC( (size_t)N_H_MOLEC*sizeof(double *) ))) == NULL )
			BadMalloc();
		for( i=0; i<N_H_MOLEC; ++i )
		{
			/* this is the Jacobian array, MALLOC sets it to NaN, first reagents
			 * will be filled in, then Jacobian array set in loop below.  Search
			 * for Jacobian array */
			if( (c[i] = ((double*)MALLOC( (size_t)N_H_MOLEC*sizeof(double) ))) == NULL )
				BadMalloc();
		}
	}

	/* Assume no error for cases with abundances set */
	*error = 0;

	/* there are two "no molecules" options, the no co, which turns off everything,
	 * and the no n2, which only turns off the h2.  in order to not kill the co
	 * part we still need to compute the hydrogen network here, and then set h2 to
	 * small values */
	/* >> chng 03 jan 15 rjrw -- suddenly switching off molecules confuses the solvers... */
	if( hmi.lgNoH2Mole ) /*   || phycon.te > 1e5 ) */
	{
		/* thtmol = 0. */
		dense.xMolecules[ipHYDROGEN] = 0.;

		/* these are the molecular species */
		for (mol=0; mol<N_H_MOLEC; ++mol) 
		{
			hmi.Hmolec[mol] = 0.;
		}
		hmi.Hmolec[ipMH] = dense.xIonDense[ipHYDROGEN][0];
		hmi.Hmolec[ipMHp] = dense.xIonDense[ipHYDROGEN][1];
		hmi.H2_total = 0.;
		/* this is where the Emline struc expects to find the H2 abundance */
		dense.xIonDense[LIMELM+2][0] = 0.;
		hmi.hmihet = 0.;
		hmi.h2plus_heat = 0.;
		hmi.H2Opacity = 0.;
		hmi.hmicol = 0.;
		hmi.HeatH2Dish_TH85 = 0.;
		hmi.HeatH2Dexc_TH85 = 0.;
		hmi.deriv_HeatH2Dexc_TH85 = 0.;
		hmi.hmidep = 1.;
		hmi.rh2dis = 0.;
		hmi.HalphaHmin = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->hmole()\n", debug_fp );
#		endif
		return;
	}

	/* >>chng 02 jun 19, add option to force H2 abundance, for testing h2 molecules,
	 * hmi.frac_abund is fraction in molecules that is set by set h2 fraction command */
	if( hmi.frac_abund>0.)
	{
		for (mol=0;mol<N_H_MOLEC;mol++) 
		{
			hmi.Hmolec[mol] = 0.;
		}
		/* >>chng 03 jul 19, from 0 to SMALLFLOAT, to pass asserts in ConvBase,
		 * problem is that ion range has not been reset for hydrogen */
		dense.xIonDense[ipHYDROGEN][0] = dense.xIonDense[ipHYDROGEN][1] = 
			2.f*SMALLFLOAT*dense.gas_phase[ipHYDROGEN];
		/* put it all in the ground state */
		hmi.Hmolec[ipMH2g] = (float)(dense.gas_phase[ipHYDROGEN] * hmi.frac_abund);
		hmi.Hmolec[ipMH2s] = 0.;

		hmi.H2_total = hmi.Hmolec[ipMH2g] + hmi.Hmolec[ipMH2s];
		/* first guess at ortho and para densities */
		h2.ortho_density = 0.75*hmi.H2_total;
		h2.para_density = 0.25*hmi.H2_total;

		hmi.hmihet = 0.;
		hmi.h2plus_heat = 0.;
		hmi.H2Opacity = 0.;
		hmi.hmicol = 0.;
		hmi.HeatH2Dish_TH85 = 0.;
		hmi.HeatH2Dexc_TH85 = 0.;
		hmi.deriv_HeatH2Dexc_TH85 = 0.;
		hmi.hmidep = 1.;
		hmi.HalphaHmin = 0.;

		for( nd=0; nd < gv.nBin; nd++ )
		{
			gv.bin[nd]->rate_h2_form_grains_used = 0.;
		}

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

	/* update these two to current values of atomic and ionized hydrogen density */
	hmi.Hmolec[ipMH] = dense.xIonDense[ipHYDROGEN][0];
	hmi.Hmolec[ipMHp] = dense.xIonDense[ipHYDROGEN][1];
	/* now copy all of H moles into old array */
	for (mol=0;mol<N_H_MOLEC;mol++) 
	{
		Hmolec_old[mol] = hmi.Hmolec[mol];
	}
	for(i=0; i<MAXREACTANTS ; ++i )
	{
		rate_deriv[i] = 0.;
		sinkrate[i] = 0.;
	}

	/* collisional ionization of H-, rate from Janev, Langer et al. */
	if( phycon.te < 3074. )
	{
		cionhm = 1.46e-32*(powi(phycon.te,6))*phycon.sqrte*hmi.exphmi;
	}
	else if( phycon.te >= 3074. && phycon.te < 30000. )
	{

		/* >>chng 03 mar 07, from above to below */
		/*cionhm = 5.9e-19*phycon.te*phycon.te*phycon.sqrte*phycon.te03*
		  phycon.te01*phycon.te01;*/
		cionhm = 5.9e-19*phycon.tesqrd*phycon.sqrte*phycon.te05;
	}
	else
	{
		cionhm = 3e-7;
	}

	/* H2 formation on grains;
	 * rate from 
	 * >>refer	H2	grain formation	Hollenback, D., & McKee, C.F., 1979, ApJS, 41, 555 eq 3.4 3.8 */
	if( gv.lgDustOn )
	{

#ifndef IGNORE_QUANTUM_HEATING
		/* hmole is called before grains, so assure that all the grain stuff is properly initialized */
		GrainDrive();
#endif
		/* 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 */
		hmi.rate_h2_allX_2_J0_grains = 0.;
		/* rate (s-1) all v,J levels go to 0, regardless of nuclear spin */
		hmi.rate_h2_ortho_para_conserve = 0. ;

		/* loop over all grain species */
		for( nd=0; nd < gv.nBin; nd++ )
		{
#ifndef IGNORE_QUANTUM_HEATING
			long k, qnbin;
			double *qtemp, *qprob;
			int lgUseQHeat = gv.lgGrainPhysicsOn && gv.bin[nd]->lgQHeat;
#endif
			/* >>chng 02 feb 15, removed check tedust > 1.01, change in InitGrains
			 * guarantees that all relevant parameters are initialized, PvH */

			/* sticking probability, 2H + grain equation 3.7 of
			 * >>refer	grain	phys	Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555,
			 * fraction of H impacts on grain surface that stick */
			/* this sticking probability is used for both HM79 and CT02 */
			double sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) + 
				0.002*phycon.te + 8e-6*phycon.te*phycon.te);

#ifndef IGNORE_QUANTUM_HEATING
			/* >>chng 04 feb 21, included quantum heating in calculation of formation rate, PvH */
			if( lgUseQHeat )
			{
				if( ( qtemp = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
					BadMalloc();
				if( ( qprob = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
					BadMalloc();

				qheat(qtemp,qprob,&qnbin,nd);

				if( gv.bin[nd]->lgUseQHeat )
				{
					ASSERT( qnbin > 0 );
				}
				else
				{
					qnbin = 1;
					qprob[0] = 1.;
					qtemp[0] = gv.bin[nd]->tedust;
				}

				gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;

				for( k=0; k < qnbin; k++ )
				{
					/* fraction of impacts that produce H2 before evaporation from grain surface.
					 * this is equation 3.4 of
					 * >>refer	grain	phys	Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555
					 * 1e4 is ratio of total absorption sites to apropriate sites 
					 * 920 is D_H and chosen to get f_a = 0.5 at 100 K.
					 * factor of 0.6252 needed to obtain std ism rate to be 3e-17 at 100 K,
					 * the value deduced by
					 * >>refer	H2	grain physics	Jura, M., 1974, ApJ, 197, 581 */
					double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./qtemp[k]));
					sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) + 
								     0.002*phycon.te + 8e-6*phycon.te*phycon.te);

					gv.bin[nd]->rate_h2_form_grains_HM79 += qprob[k] * sticking_probability_H *
						conversion_efficiency_HM79 ;
				}

				/* NB IntArea is total, not projected, area, must div by 4 */
				/* gv.bin[nd]->rate_h2_form_grains_HM79 has units s^-1 since gv.bin[nd]->cnv_H_pCM3 has units cm-3 */
				/* cnv_H_pCM3 converts <unit>/H (default depl) -> <unit>/cm^3 (actual depl), units are cm-3 */
				gv.bin[nd]->rate_h2_form_grains_HM79 *= 0.5 * DoppVel.AveVel[ipHYDROGEN]*
					gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;

				ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
			}
			else
#endif
			{
				/* fraction of impacts that produce H2 before evaporation from grain surface.
				* this is equation 3.4 of
				* >>refer	grain	phys	Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555
				* 1e4 is ratio of total absorption sites to apropriate sites 
				* 920 is D_H and chosen to get f_a = 0.5 at 100 K.
				* factor of 0.6252 needed to obtain std ism rate to be 3e-17 at 100 K,
				* the value deduced by
				* >>refer	H2	grain physics	Jura, M., 1974, ApJ, 197, 581 */
				double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./gv.bin[nd]->tedust));

				/* NB IntArea is total, not projected, area, must div by 4 */
				/* units s^-1 since gv.bin[nd]->cnv_H_pCM3 has units cm-3 */
				gv.bin[nd]->rate_h2_form_grains_HM79 = 0.5 * DoppVel.AveVel[ipHYDROGEN]* gv.bin[nd]->IntArea/4. * 
					/* cnv_H_pCM3 converts <unit>/H (default depl) -> <unit>/cm^3 (actual depl), units are cm-3 */
					gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * conversion_efficiency_HM79 ;
				ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
			}

#ifndef IGNORE_QUANTUM_HEATING
			if( lgUseQHeat )
			{
				/* H2 formation on grains from 
				 * >>refer	H2	form	Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */
				/* number of monolayers per second - only affects efficienty at very low or high temperatures */
				double f = 1e-10;
				/* equation 17 
				double sqrt_term = POW2( 1. + sqrt( (10000.-200.)/(600.-200.) ) );*/
				double sqrt_term = 35.399494936611667;

				gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;

				for( k=0; k < qnbin; k++ )
				{
					double beta_alpha = 0.25 * sqrt_term *sexp(200./qtemp[k] );
					/* equation 16 */
					double xi =  1./ (1. + 1.3e13*sexp(1.5*1e4/qtemp[k])*sqrt_term/(2.*f) ) ;
					/* expression for beta comes from just after equation 5 */
					double beta = 3e12 * sexp( 320. / qtemp[k] );
					/* recombination efficiency given by their equation 15, they call
					 * this epsilon_H2 */
					double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );
					sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) + 
								     0.002*phycon.te + 8e-6*phycon.te*phycon.te);

/* 					printf( " k %ld Td %.6e re*sp %.6e\n", k, qtemp[k], recombination_efficiency_CT02* */
/* 						sticking_probability_H ); */

					gv.bin[nd]->rate_h2_form_grains_CT02 += qprob[k] * sticking_probability_H *
						recombination_efficiency_CT02;
				}

				/* gv.bin[nd]->IntArea integrated grain surface area Int(4pi*a^2), normalized per H, in cm^2/H,
				 * so x/4 is projected area of circle */
				/* gv.bin[nd]->cnv_H_pCM3 is H density [cm-3] times grain depletion factor */
				/* gv.bin[nd]->rate_h2_form_grains_CT02 units s-1 */
				gv.bin[nd]->rate_h2_form_grains_CT02 *= 0.5 * DoppVel.AveVel[ipHYDROGEN]*
					gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;

				ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );

				free(qtemp);
				free(qprob);
			}
			else
#endif
			{
				/* H2 formation on grains from 
				 * >>refer	H2	form	Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */
				/* number of monolayers per second - only affects efficienty at very low or high temperatures */
				double f = 1e-10;
				/* equation 17 
				double sqrt_term = POW2( 1. + sqrt( (10000.-200.)/(600.-200.) ) );*/
				double sqrt_term = 35.399494936611667;
				double beta_alpha = 0.25 * sqrt_term *sexp(200./gv.bin[nd]->tedust );
				/* equation 16 */
				double xi =  1./ (1. + 1.3e13*sexp(1.5*1e4/ gv.bin[nd]->tedust )*sqrt_term/(2.*f) ) ;
				/* expression for beta comes from just after equation 5 */
				double beta = 3e12 * sexp( 320. / gv.bin[nd]->tedust );
				/* recombination efficiency given by their equation 15, they call
				 * this epsilon_H2 */
				double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );

				/* gv.bin[nd]->IntArea integrated grain surface area Int(4pi*a^2), normalized per H, in cm^2/H,
				 * so x/4 is projected area of circle */
				/* gv.bin[nd]->cnv_H_pCM3 is H density [cm-3] times grain depletion factor */
				/* units s-1 */
				gv.bin[nd]->rate_h2_form_grains_CT02 = 0.5 * DoppVel.AveVel[ipHYDROGEN]* gv.bin[nd]->IntArea/4. * 
					gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_CT02 ;
				ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );
			}

#ifndef IGNORE_QUANTUM_HEATING
			/* reset sticking probability for code below */
			sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) + 
						     0.002*phycon.te + 8e-6*phycon.te*phycon.te);
#endif

			/*TODO	1	code below is inside loop over nd, but doesn't sum over nd, so hmi.rate_h2_allX_2_J0_grains
				  * and hmi.rate_h2_ortho_para_conserve are determined only by last grain bin, is this OK ???? */
			{
				/* ortho to para on grain surfaces, taken from 
				 * Le Bourlot, J., 2000, A&A, 360, 656-662 */
				/* For all grain temperatures, this process corresponds to high J going to
				 * either 0 or 1 preserving nuclear spin.  All ortho go to 1 and para go to 0.
				 * When the dust temperature is below Tcrit all 1 go to 0 and so all J go to 0.
				 */
				/* this temperature depends on grain composition, discussion left column of page 657,
				 * this is for a bare grain */
				/*TODO	2	- put in actual composition dependent Tad - this is only valid 
				 * for bare surfaces - not ice - for ice Tad is 555K */
				/* hmi.Tad is binding energy expressed as a temperature */
				/*double Tad = 800.;*/
				/* tau_nu the first equation in section 2.5
				double tnu = sqrt( BOLTZMANN / PROTON_MASS ) / PI / 2.6e-8; */
				/* equation one paragraph before equation 2 */
				/* ortho2para_final is ratio of ortho (J=1) to para (J=0) in final product,
				 * at low grain temperatures all end in para, J=1 */
				/* this is roughly 20 K */
				double T_ortho_para_crit = 2. * hmi.Tad / log( POW2(60. *1.1e11)*hmi.Tad); 

				/* AveVEl[LIMELM+2] is average speed of H2 molecules */
				/* for now assume that sticking probability for H2 on the grain is equal to
				 * that for H */
				if( gv.bin[nd]->tedust < T_ortho_para_crit )
				{
					/* rate (s-1) all v,J levels go to 0, regardless of nuclear spin */
					hmi.rate_h2_allX_2_J0_grains = DoppVel.AveVel[LIMELM+2]* gv.bin[nd]->IntArea/4. * 
						gv.bin[nd]->cnv_H_pCM3 * Hmolec_old[ipMH2g]/dense.gas_phase[ipHYDROGEN] * 
						sticking_probability_H  * mole.lgColl_deexec_Calc  ;

					/* rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
					hmi.rate_h2_ortho_para_conserve = 0. ;
				}
				else
				{
					/* rate (s-1) all v,J levels go to 0, regardless of nuclear spin */
					hmi.rate_h2_allX_2_J0_grains = 0.;

					/* rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
					hmi.rate_h2_ortho_para_conserve = DoppVel.AveVel[LIMELM+2]* gv.bin[nd]->IntArea/4. * 
						gv.bin[nd]->cnv_H_pCM3 * Hmolec_old[ipMH2g]/dense.gas_phase[ipHYDROGEN] * 
						sticking_probability_H * mole.lgColl_deexec_Calc ;
				}
			}
		}
		/*fprintf(ioQQQ," H2 grain form rate HM79 %.2e  %.2e CT02 %.2e  %.2e O-P grn %.2e %.2e\n", 
			gv.bin[nd]->rate_h2_form_grains_HM79 , 
			gv.bin[nd]->rate_h2_form_grains_HM79 ,
			gv.bin[nd]->rate_h2_form_grains_CT02 , 
			gv.bin[nd]->rate_h2_form_grains_CT02 , 
			hmi.rate_h2_allX_2_J0_grains,
			hmi.rate_h2_allX_2_J1_grains
			);*/
	}
	else
	{
		for( nd=0; nd < gv.nBin; nd++ )
		{
			gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;
			gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;
		}
		/* rate all H2 goes to either 0 or 1 depending on ortho/para */
		hmi.rate_h2_ortho_para_conserve = 0. ;
		/* at low temp, rate all H2 goes to J=0 */
		hmi.rate_h2_allX_2_J0_grains = 0.;
	}

	/* this is the one that is actually used in calculations */
	gv.rate_h2_form_grains_used_total = 0.;
	for( nd=0; nd < gv.nBin; nd++ )
	{
		if( hmi.chJura == 'T' )
		{
			/* activate first one to use Hollenbach & McKee 1979 */
			gv.bin[nd]->rate_h2_form_grains_used = hmi.grain_form_factor *
				gv.bin[nd]->rate_h2_form_grains_HM79;
			/* total rate H2 forms on grains, this is used in chem network */
			/* units are s-1*/
			gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
		}
		else if( hmi.chJura == 'C' )
		{
			/* use the new rate by Cazaux & Tielens */
			/* units are s-1*/
			gv.bin[nd]->rate_h2_form_grains_used = hmi.grain_form_factor *
				gv.bin[nd]->rate_h2_form_grains_CT02;
			/* total rate H2 forms on grains, this is used in chem network */
			/* units are s-1*/
			gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
		}
		else if( hmi.chJura == 'S' )
		{
			/* H2 formation rate from 
			 * >>refer	H2	form	Sternberg, A. & Neufeld, D.A. 1999, ApJ, 516, 371 */
			/* grain_form_factor is scale factor to change formation rate - normally unity */
			gv.bin[nd]->rate_h2_form_grains_used = hmi.grain_form_factor *
				3e-18 * phycon.sqrte / gv.nBin * dense.gas_phase[ipHYDROGEN];
			/* this is simple rate from Sternberg & Neufeld 99 */
			/* total rate H2 forms on grains, this is used in chem network */
			/* units are s-1*/
			gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
		}
	}
	ASSERT( gv.rate_h2_form_grains_used_total >= 0. );

#ifndef IGNORE_QUANTUM_HEATING
	printf( " fnzone %.2f H2 rate %.4e\n", fnzone, gv.rate_h2_form_grains_used_total );
#endif

	/* >>chng 03 sep 09, get ratio of excited to ground state H2 */
	if( h2.lgH2ON  && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
	{
		frac_H2star_grain = hmi.H2star_forms_grain / 
			SDIV(hmi.H2star_forms_grain+hmi.H2_forms_grain);
		/*printf( "H2s frac grain %.4e\n",frac_H2star_grain  );*/

		frac_H2star_hminus = hmi.H2star_forms_hminus / 
			SDIV(hmi.H2star_forms_hminus+hmi.H2_forms_hminus);
		/*printf( "H2s frac h- %.4e\n",frac_H2star_hminus  );*/
	}
	else
	{
		/* the large H2 molecule was not evaluated, so we can't use exact
		 * branching ratios.  These are the distribution fractions for around 500K */
		/*These depend on temperature and distribution function and the definition of ENERGY_H2_STAR.
		  So reset the values properly*/
		frac_H2star_grain = 0.51f;
		frac_H2star_hminus = 0.83f;
	}

	/* print rate coefficient */
	/*fprintf(ioQQQ," total grain h2 form rate %.3e\n",gv.rate_h2_form_grains_used_total);*/

	/* collisional dissociation, rate from 
	 * >>refer	H2	collisional dissociation	Dove, J.E., and Mandy, M. E., 1986, ApJ, 311, L93.
	 * corr is correction for approach to high density limit
	 * H2 + H => 3H - rate very uncertain */
	corr = MIN2(6.,14.44-phycon.alogte*3.08);

	if (corr > 0.)
		corr = pow(10.,corr*Hmolec_old[ipMH]/(Hmolec_old[ipMH]+1.6e4));
	else
		corr = 1.;
	hmi.rh2dis = (float)(1.55e-8/phycon.sqrte*sexp(65107./phycon.te)* corr);

	/* old hminus rate Hollenbach & McKee 1979
	 *>>chng 98 jan 02, from 2.12e4 to 2.123e4 */
	/*hmi.bh2h2p = 1.8e-12f*phycon.sqrte*phycon.te10/phycon.te01*2.f/16.f;
	rh2h2p = 1.8e-12*phycon.sqrte*phycon.te10/phycon.te01*sexp(2.123e4/
	  phycon.te);*/

	/* forward and back reactions for H2+ + H+ <=> H2+ + H */
	/*>>chng 02 oct 25, update rate from above (Hollenbach & McKee 1979) to
	 * >>refer	H2	form	Karpas, Z., Anicich, V., & Huntress, W.T. 1979, J Chem Phys, 70, 2877 
	 * following is from this paper.\:
	We note that the application of detailed balance is only strictly valid for 
	state-to-state reactions, i.e., when the , J level of the reactant and product 
	molecules are known. While typical laboratory conditions are such that  = 0 and J 
	is likely to be small for the reactant, the product , J is usually unknown. 
	For example, Krsti (2002) finds that for the reverse of reaction (1), the product 
	H2 is primarily formed into  = 4, not  = 0. Therefore, estimation of reaction (1) 
	by the application of detailed balance to the measured rate coefficient for the 
	reverse reaction gives the rate coefficient for H2( = 4), which can be as much as 
	an order of magnitude larger than for H2( = 0). Rate coefficients that are estimated 
	by detailed balance are therefore suspect. 
	*/

	/* H2+ + H => H2 + H+ */
	hmi.bh2h2p = 6.4e-10f;
#	if 0
	/* H2  +  H+  =>  H2+  +  H
	 * >>chng 04 Feb 24, get back reaction from above */
	 if (hmi.rel_pop_LTE_H2 != 0.)
		rh2h2p = hmi.bh2h2p*hmi.rel_pop_LTE_H2p/hmi.rel_pop_LTE_H2;
	 else
		rh2h2p = 0;
#	endif
	/*>> refer	H2	chem	Savin, D.W., Krstic, P.S., Haiman, Z., & Stancil, P.C., 2004,
	 *>>refercon	ApJL, 606L, 167, astro-ph/0404288 */
	if (phycon.te<=3.e4 )
	{
		/* this is lower bound to their temperature */
		/*TODO check that this is really natural log, not base 10 */
		double teused = MAX2(100., phycon.te );
		double telog = log(teused);
		rh2h2p = sexp(2.123715e4/teused)*(-3.3232183e-7 + 3.3735382e-7*log(teused) -
			1.4491368e-7*pow(telog,2) +3.4172805e-8*pow(telog,3) -
			4.781372e-9*pow(telog,4) + 3.9731542e-10*pow(telog,5) -
			1.8171411e-11*pow(telog,6) +3.5311932e-13*pow(telog,7));
	}
	else
		rh2h2p= 0;

	/* H2+  +  HNU  =>  H+  + H */
	gamtwo = GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.);
	
	/* this only occurs when set umits rates is entered */
	if (!co.lgUMISTrates)
		gamtwo = 5.7e-10*hmi.UV_Cont_rel2_Habing_TH85_face*(float)sexp((1.9*rfield.extin_mag_V_point))/1.66f;

	/*GammaPrt(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,ioQQQ,gamtwo,0.01);*/

	h2phet = thermal.HeatNet;

	/* >> chng 02 nov 15 rjrw: ionization fractions to multiply c[ipHo][*] terms
	 * as b[ipMHo] contains _both_ H0 and H+ */
	/* sum_H0_Hp = ((double)dense.xIonDense[ipHYDROGEN][0])+((double)dense.xIonDense[ipHYDROGEN][1]); */

	sum_H0_Hp = Hmolec_old[ipMH]+Hmolec_old[ipMHp];

	for (mol=0;mol<N_H_MOLEC;mol++)
		nprot[mol] = hmi.nProton[mol];

	/* rtot is total absolute value of reaction rates, for
	 *checking rate balance */

	rtot = 0.;
	rindex = 0;
	r = rlist;
	/* Special case, put null reaction at head of list */
	if (r == NULL) {
		int in[]={-1},out[]={-1};
		r = rlist = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	rindex++;

	/*-------------------------------------------------------------------- */

	/* H- H minus hminus balance equations
	 * (IHMI,IPMHO) == processes making H- from Ho =+sign
	 * radiative attachment: HI + NE => H-
	 * H + e -> H- + hnu */
	/* Use Newton-Raphson step to improve solution, so bvec[] contains reaction rates
	 * and c[][] components of the Jacobian of the rates */

	/* This block adds a reaction H => H- to the stack if it wasn't
	 * there already.
	 *
	 * >>>> ONLY CHANGE the elements of the in[] and out[] vectors and
	 * the rate constant, keep the rest fixed for all reactions 
	 * */
	if (r->next == NULL) {
		int in[]={ipMH},out[]={ipMHm};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = (hmi.hminus_rad_attach + hmi.HMinus_induc_rec_rate)*co.lgUMISTrates;

	/* >>chng 02 oct 29, add these two chemical processes */
	/* H- + H+ => H2+ + e
	 * equation (H6) from 
	 * >>refer	H2	chemistry	Galli,D., & Palla, F. 1998, A&A, 335, 403-420 
	 * hmihph2p = 6.9e-9f*(Tg)^(-0.35) for Tg<=8000
	 * hmihph2p = 6.9e-9f*(Tg)^(-0.9) for Tg>=8000  */
	/* >>chng 02 nov 07 rjrw, include H+ ion density in rate constant */
	if (phycon.te <= 7891.)
	{
	    /*hmihph2p = 6.9e-9*pow(phycon.te , -0.35);*/
	    hmihph2p = 6.9e-9 / (phycon.te30 * phycon.te05) ;
	}
	else 
	{
		/* >>chng 02 nov 18, had typo for leading coef here */
		/*hmihph2p = 9.6e-7*pow(phycon.te , -0.9);*/
		hmihph2p = 9.6e-7 / phycon.te90;
	}

	if (r->next == NULL) {
		int in[]={ipMHm,ipMHp},out[]={ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmihph2p*co.lgUMISTrates;

	/* >>chng 03 feb 6 */
    /* H2+ + H- => H2 + H
	 * equation (32) from
	 * >>refer  Stancil, P.C, & Lepp, S, & Dalgarno, A. 1998,ApJ, 509, 1-10
	 * h2phmh2h = 1.4e-7f*pow(phycon.te/300.0, -0.5) */
	/* >>chng 03 sep 01, rm the pow function */
	/*h2phmh2h = 1.4e-7f*pow(phycon.te/300 , -0.5);*/
	/* the fits in this paper cannot be used below 10K or above 1e4K.  Limit
	 * the range of evaluation */
	TStancil = MIN2(phycon.te, 1e4 );
	TStancil = MAX2( TStancil , 10. );

	h2phmh2h = 1.4e-7*17.305/phycon.sqrte;
    if (r->next == NULL) {
		int in[]={ipMH2p,ipMHm},out[]={ipMH2g,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2phmh2h*co.lgUMISTrates;



    /* >>chng 03 feb 7 */
    /* H2+ + H- => H + H + H
	 * equation (33) from
	 * >>refer  Stancil, P.C, & Lepp, S, & Dalgarno, A. 1998,ApJ, 509, 1-10
	 * h2phmhhh = 1.4e-7f*pow(phycon.te/300.0, -0.5) */
	/* >>chng 03 sep 01, rm the pow function */
	/*h2phmhhh = 1.4e-7f*pow(phycon.te/300 , -0.5);*/
	h2phmhhh = 1.4e-7f*17.3205/phycon.sqrte;
	if (r->next == NULL) {
		int in[]={ipMH2p,ipMHm},out[]={ipMH,ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2phmhhh;

	/* >>chng 03 sep 30 */
	/* H+ + H- => H + H+ + e
	 * >>refer Paolo Lenzuni, David F. Chernoff, Edwin E. Salpeter, 1991, ApJS, 76, 759L  (table 5)
	 * hphmhhpe = 4.75e-30*pow(phycon.te,3.1); */
	hphmhhpe = 4.75e-30*pow(phycon.te,3.1);
	if (r->next == NULL) {
		int in[]={ipMHp,ipMHm},out[]={ipMH,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =hphmhhpe ;

	/* >>chng 03 sep 30 */
	/* H2 + H- => H + H2 + e
	 * >>refer Paolo Lenzuni, David F. Chernoff, Edwin E. Salpeter, 1991, ApJS, 76, 759L (table 5)
	 * h2hmhh2e = 6.74e-17*pow(phycon.te,2)*sexp(19870/phycon.te); */
	h2hmhh2e = 6.74e-17*phycon.tesqrd*sexp(19870/phycon.te);
	if (r->next == NULL) {
		int in[]={ipMH2g,ipMHm},out[]={ipMH,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =h2hmhh2e ;

	/* (IHMI,IHMI) = processes destroying H- =-sign
	 * photodissociation, H- + H NU => H + NE */

	if (r->next == NULL) {
		int in[]={ipMHm},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.HMinus_photo_rate;

	/* mutual neutralization with heavies, rate from Dalgarno and Mcray
	 * all charged ions contribute equally,
	 * H- + A+ => H + A */
	/* >>chng 04 feb 19, do actual sum of first ions rather than following kludge */
	/* find the sum of all single ion densities for species heavier than helium */
	sum_first_ions = 0.;
	for( i=ipLITHIUM; i < LIMELM; i++ )
	{
		sum_first_ions += dense.xIonDense[i][1];
	}
#	if 0
	if( nzone > 140 )
	{
		fprintf(ioQQQ,"DEBUG sumfirstions\t%.2f\t%.4e\t%.4e\t%.4e",
			fnzone,phycon.te,
			sum_first_ions,
			sum_first_ions/dense.eden);
		for( i=ipLITHIUM; i < LIMELM; i++ )
		{
			if( dense.xIonDense[i][1]/sum_first_ions >0.1 )
				fprintf(ioQQQ,"\t%li\t%.3e",
					i,dense.xIonDense[i][1]/sum_first_ions);
		}
		fprintf(ioQQQ,"\n");
	}
#	endif
	/*hmi.hmin_ct_firstions = dense.eden-Hmolec_old[ipMHp]-dense.xIonDense[ipHELIUM][1]-
		2.*dense.xIonDense[ipHELIUM][2];*/
	hmi.hmin_ct_firstions = 4e-6f/phycon.sqrte;

	/* hmi.hmin_ct_firstions = 4e-6/phycon.sqrte*MAX2(0.,dense.eden-dense.xIonDense[ipHYDROGEN][1]-dense.xIonDense[ipHELIUM][1]-
		 2.*dense.xIonDense[ipHELIUM][2]); */

	if (r->next == NULL) {
		int in[]={ipMHm},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.hmin_ct_firstions*sum_first_ions;

	/* electron collisional ionization of H- */
	cionhm *= dense.eden;

	if (r->next == NULL) {
		int in[]={ipMHm},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = cionhm;

	/* inverse process; three body rec */
	c3bod = cionhm*(hmi.rel_pop_LTE_Hmin*dense.eden);

	if (r->next == NULL) {
		int in[]={ipMH},out[]={ipMHm};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = c3bod;

	/* form molecular hydrogen from H minus,
	 * associative detachment:  H- + H => H2 + E */
	/* make H2 from H- 
	 * associative detachment; H- + H => H2: 
	 * >>referold	H2	rates	Browne & Dalgarno J PHys B 2, 885 */
	/* rate coefficient from 
	 * >>refer	H2	form	Launay, J.R., Le Dourneuf, M., & Zeippen, C.J., 
	 * >>refercon	1991, A&A, 252, 842-852*/
	/* >>chng 02 oct 17, temp dependent fit to rate, updated reference,
	 * about 40% larger than before */
	{
		double y , x;
		x = MAX2(10., phycon.te );
		x = MIN2(1e4, x );
		y=545969508.1323510+x*71239.23653059864;
		hmi.assoc_detach = 1./y;
	}

	/* >>chng 02 nov 7 rjrw, example case of 2-body process */
	/* this one is into ground H2 */
	if (r->next == NULL) {
		int in[]={ipMH, ipMHm},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.assoc_detach*(1.-frac_H2star_hminus);
	
	/* >>chng 03 sep 10, multiply above by correction for excited state,
	 * add below reaction for population of excited state */
    /* this one is into excited H2 */
	if (r->next == NULL) {
		int in[]={ipMH, ipMHm},out[]={ipMH2s};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.assoc_detach*frac_H2star_hminus*co.lgUMISTrates;

	{
		/*@-redef@*/
		/* this debug print statement compares H2 formation through grn vs H- */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nzone>140 )
		{
			fprintf(ioQQQ," debuggggrn grn\t%.2f\t%.3e\t%.3e\tfrac\t%.3e\tH-\t%.3e\t%.3e\tfrac\t%.3e\t%.3e\t%.3e\t%.3e\n",
				fnzone ,
				gv.rate_h2_form_grains_used_total , 
				hmi.H2_forms_grain+hmi.H2star_forms_grain ,
				frac_H2star_grain,
				hmi.assoc_detach*dense.xIonDense[ipHYDROGEN][0]*hmi.Hmolec[ipMHm],
				hmi.H2star_forms_hminus+hmi.H2_forms_hminus,
				frac_H2star_hminus,
				hmi.assoc_detach,dense.xIonDense[ipHYDROGEN][0],hmi.Hmolec[ipMHm]
				);
		}
	}

	/* convert H2 into H- 
	 * the back reaction, H2(grnd) + e => H- + Ho */
	if( hmi.rel_pop_LTE_H2 > 0. )
	{ 
		assoc_detach_backwards_grnd = hmi.assoc_detach*hmi.rel_pop_LTE_Hmin/hmi.rel_pop_LTE_H2*dense.eden;
	}
	else
	{
		assoc_detach_backwards_grnd = 0.;
	}
	/* convert H2 into H- 
	 * the back reaction, H2(exct) + e => H- + Ho */
    if( hmi.rel_pop_LTE_H2s > 0. )
	{ 
		
	    assoc_detach_backwards_exct = hmi.assoc_detach*hmi.rel_pop_LTE_Hmin/hmi.rel_pop_LTE_H2s*dense.eden; 
	}
	else
	{
		assoc_detach_backwards_exct = 0.;
	}
	{
		/*@-redef@*/
		/* often the H- route is the most efficient formation mechanism for H2,
		 * will be through rate called Hmolec_old[ipMH]*hmi.assoc_detach
		 * this debug print statement is to trace h2 oscillations */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nzone>140/*&& iteration > 1*/)
		{
			/* rapid increase in H2 density caused by rapid increase in hmi.rel_pop_LTE_H2 */
			fprintf(ioQQQ,"assoc_detach_backwards_grnd\t%.2f\t%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
				/* total forward rate */
				fnzone,
				phycon.te, 
				dense.eden,
				/* rate H- + H => H2 + E */
				hmi.assoc_detach,
				assoc_detach_backwards_grnd, 
				assoc_detach_backwards_exct, 
				hmi.hminus_rad_attach,
				hmi.HMinus_induc_rec_rate,
				/* H0 */
				hmi.Hmolec[ipMH],
				/* H+ */
				hmi.Hmolec[ipMHp],
				/* H- */
				hmi.Hmolec[ipMHm],
				hmi.H2_total,
				hmi.rel_pop_LTE_Hmin,
				hmi.rel_pop_LTE_H2,
				hmi.rel_pop_LTE_H2s
				);
		}
	}

	/* >>chng 03 sep 11, resolve H2 and H2*, use fraction determined above */
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH,ipMHm};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = assoc_detach_backwards_grnd*(1.-frac_H2star_hminus)*co.lgUMISTrates;

	/* >>chng 03 sep 11, resolve H2 and H2*, add new destruction process for H2* */
	if (r->next == NULL) {
		int in[]={ipMH2s},out[]={ipMH,ipMHm};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	/* >>chng 04 jan 28, had wrong Boltzmann factor for this reaction,
	 * fixed by Gargi Shaw */
	r->ratek = assoc_detach_backwards_exct*frac_H2star_hminus*co.lgUMISTrates;

#define HNEUT 7e-8
	/* mutual neut, mostly into n=3; rates from Janev et al
	 * H- + H+ => H + H(n=3) */
	/*TODO	2	process is net source term for H(n=3) states, must be added in */
	fhneut = Hmolec_old[ipMHp]*HNEUT; /* dense.xIonDense[ipHYDROGEN][1]*7e-8; */

	if (r->next == NULL) {
		int in[]={ipMHm,ipMHp},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = HNEUT*co.lgUMISTrates;

	/* back reaction from excited state H */
	if( phycon.te > 1000. )
	{
		/* HBN(3,1) is defined; when <HydTempLimit then set to 1 */
		bhneut = (HNEUT*hmi.rel_pop_LTE_Hmin*dense.eden)*iso.DepartCoef[ipH_LIKE][ipHYDROGEN][3];
	}
	else
	{
		bhneut = 0.;
	}

	/* mutual neut, mostly into n=3; rates from Janev et al
	 * H + H(n=3) => H- + H+ */
	/*TODO	2	process is net ionization term for H(n=3) states */
	/* this is the back reaction, forming H- from Ho */

	if (r->next == NULL) {
		int in[]={ipMH,ipMH},out[]={ipMHm,ipMHp}, ratesp[]={ipMH,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = bhneut; 
	bhneut *= Hmolec_old[ipMHp];


	/* the processes H2(v>=4) + e- => H + H-,
	 * >>refer	H2	rates	Lenzuni et al. apj sup 76, 759
	 * >>refer	H2	rates	Janev et al.
	 * density dep is for non-lte, guess from dalgarno and roberge apl 233, 25
	 * extra expo factor added for low temps */
	if( nzone <= 1 )
	{
		/* this is initial setup of code, so set rate coef to actual val */
		eh2hhm = 2.7e-8*pow(10.,-0.7*POW2(phycon.alogte - 3.615))*
		  dense.eden*(dense.gas_phase[ipHYDROGEN]/(1e7 + dense.gas_phase[ipHYDROGEN]))*sexp(52000./
		  phycon.te);
		eh2old = eh2hhm;
	}
	else
	{
		/* this is deeper into the cloud, and there is danger of oscillation */
		eh2old = eh2hhm;
		eh2hhm = 2.7e-8*pow(10.,-0.7*POW2(phycon.alogte - 3.615))*
		  dense.eden*(dense.gas_phase[ipHYDROGEN]/(1e7 + dense.gas_phase[ipHYDROGEN]))*sexp(52000./
		  phycon.te);
		fac = 0.5;
		eh2hhm = eh2hhm*fac + eh2old*(1. - fac);
	}

	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH,ipMHm};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = eh2hhm;

	/*--------------------------------------------------------------------
	 *
	 * molecular hydrogen H2 htwo balance equation
	 * (IPMH2,IPMHO)==create H2 from Ho =+ */

	/* H2 formation on grains */
	/* >>chng 01 jan 05, remove from matrix part and add hden to hmi.rate_h2_form_grains_used, */
	/* the large molecule keeps explicit track of the fraction that goes into 
	 * excited vs ground H2.  Use that ratio if H2 turned on, else use an
	 * estimate of it */
	/* The reaction rate is only proportional to one of the ipMH, due to
	 * surface saturation (?) */
	
#define	CATALYST	TRUE
	if( CATALYST )
	{
		/* this is the method used by the code for most of its history.  The grain
		 * is only a catalyttic agent, and so the rate goes as the square of the
		 * incoming H0 density */
		/* This goes to excited H2 */
		if (r->next == NULL) {
			int in[]={ipMH,ipMH},out[]={ipMH2s},ratesp[]={ipMH};
			r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
		}
		r = r->next;
		rindex++;
		r->ratek = gv.rate_h2_form_grains_used_total*frac_H2star_grain;
		
		/* >>chng 03 sep 10, multiply above by correction for excited state,
		* add below reaction for population of excited state */

		/* This goes to ground H2 */
		if (r->next == NULL) {
			int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH};
			r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
		}
		r = r->next;
		rindex++;
		r->ratek = gv.rate_h2_form_grains_used_total*(1. - frac_H2star_grain);

	}
	else
	{
		/* >>chng 03 nov 25, go to this formalism */
		/* the grain is not a true catalyst, but rather a target loaded with H atoms
		 * ready to react.  So the rate is the number of these grains, times their
		 * cross section, times the number of incident H atoms.  The number of grains
		 * is replaced with the total hydrogen density, which is not in the network
		 * but is a constant */
		/* This goes to excited H2 */
		if (r->next == NULL) {
			int in[]={ipMH,ipMH},out[]={ipMH2s},ratesp[]={ipMH};
			r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
		}
		r = r->next;
		rindex++;
		r->ratek = gv.rate_h2_form_grains_used_total*frac_H2star_grain*
			dense.gas_phase[ipHYDROGEN]/SDIV(dense.xIonDense[ipHYDROGEN][0]);
		
		/* >>chng 03 sep 10, multiply above by correction for excited state,
		* add below reaction for population of excited state */

		/* This goes to ground H2 */
		if (r->next == NULL) {
			int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH};
			r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
		}
		r = r->next;
		rindex++;
		r->ratek = gv.rate_h2_form_grains_used_total*(1. - frac_H2star_grain)*
			dense.gas_phase[ipHYDROGEN]/SDIV(dense.xIonDense[ipHYDROGEN][0]);

	}

	/* excited atom radiative association,
	 * H(n=2) + H(n=1) => H2 + hnu
	 * written as H(n=1)*pop ratio + H(n=1) -> H2 + hnu but ratio of pops is in
	 * terms of pop to ion, so initial term is
	 * n(H+) * pop2ion, 
	 * >>refer	H2	rates	Latter, W.B., & Black, J.H., 1991, Ap.J. 372, 161 */
	/* radasc = ((iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p] + iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2s])*dense.xIonDense[ipHYDROGEN][1])*3e-14; */
	
	radasc = ((iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p] + iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2s]))*3e-14;
	/* >>chng 02 nov 7 rjrw: correct for n^2 behaviour w.r.t. H 
	   >>chng 02 nov 7 rjrw, correct stoichiometry */

	/* Possible that changing to a rate proportional to ipMHp would be more consistent */
	if (r->next == NULL) {
		int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = radasc;  
	radasc *= Hmolec_old[ipMHp];

	/* photo-destroy H2 */
	/* >>chng 00 nov 25 factor of 0.1, assume pump is total, and 10% distroy H2 is 21*/
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	/* >>chng 03 mar 07, had factor of 0.1 for branching ratio from H2** to H+H, 
	 * but branching is now already included */
	/*r->ratek = hmi.H2_Solomon_dissoc_rate_used*0.1;*/
	r->ratek = hmi.H2_Solomon_dissoc_rate_used_H2g;

	/* >>chng 03 sep 11, add this process */
    /* photo-destroy H2* by Solomon process at same rate as H2ground dissociation,
		see above eqn A12 in TH85 */
	/* >>chng 00 nov 25 factor of 0.1, assume pump is total, and 10% distroy H2 */
	if (r->next == NULL) {
		int in[]={ipMH2s},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	/* >>chng 03 mar 07, had factor of 0.1 for branching ratio from H2** to H+H, 
	 * but branching is now already included */
	/*r->ratek = hmi.H2_Solomon_dissoc_rate_used*0.1; is #22*/
	r->ratek = hmi.H2_Solomon_dissoc_rate_used_H2s;

	/* the process H2 + e- => H + H + e
	 * >>refer	H2	rates	Lenzuni et al. apj sup 76, 759, quoted from Janev et al. 
	 * >>chng 02 nov 7 rjrw, correct stoichiometry */
	eh2hh = 1.3e-18*phycon.te*phycon.te*sexp(52000./phycon.te)*dense.eden;
	
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = eh2hh;

	/* H2 + H+ => H3+ HNU
	 * equation H21 from 
	 * >>refer	H2	chemistry	Galli,D., & Palla, F. 1998, A&A, 335, 403-420 */
	/* >>chng 02 nov 07 rjrw, include H+ ion density in rate constant */
	h2hph3p = 1.0e-16f;

	if (r->next == NULL) {
		int in[]={ipMH2g,ipMHp},out[]={ipMH3p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2hph3p*co.lgUMISTrates;

	/* collisional dissociation, rate from 
	 * >>refer	H2	collisional dissociation	Dove, J.E., and Mandy, M. E., 1986, ApJ, 311, L93.
	 * H_2 + H => 2H + H 	 
	 * >>chng 02 nov 7 rjrw, correct stoichiometry */

	/* Rate is catalyzed by an additional H */
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH,ipMH},ratesp[]={ipMH,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.rh2dis;
  
	/* >>chng 04 apr 21 */
	/* 2H + H2 => H2 + H2
	 * equation (5) from 
	 * >>refer	H2	chemistry Palla, F., Salpeter, E.E., & Stahler, S.W., 1983, ApJ,271, 632-641
	 * bh2h22hh2= 5.5e-29/(8*phycon.te) */
	bh2h22hh2 = 5.5e-29/(8.*phycon.te);

	if (r->next == NULL) {
		int in[]={ipMH,ipMH,ipMH2g},out[]={ipMH2g,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = bh2h22hh2;

	/* >>chng 04 apr 21 */
	/* H2 + H2 => 2H + H2
	 * equation (5) from 
	 * >>refer	H2	chemistry Palla, F., Salpeter, E.E., & Stahler, S.W., 1983, ApJ,271, 632-641
	 * h2h22hh2 = bh2h22hh2/hmi.rel_pop_LTE_H2 */
	if( hmi.rel_pop_LTE_H2 > 0. )
	{
		h2h22hh2 = bh2h22hh2/hmi.rel_pop_LTE_H2;
	}
	else
	{
		h2h22hh2 =0.;
	}

	if (r->next == NULL) {
		int in[]={ipMH2g,ipMH2g},out[]={ipMH,ipMH,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2h22hh2;


	/* back rate, three body recombination, 2H + S => H_2 + S 	 
	 * >>chng 02 nov 7 rjrw: correct for n^2 behaviour w.r.t. H 
	 * >>chng 02 nov 7 rjrw, correct stoichiometry 
	 * >>chng 02 nov 7 rjrw, correct for n^3 behaviour w.r.t. H !! */
	/* bh2dis = hmi.rh2dis*hmi.rel_pop_LTE_H2*dense.xIonDense[ipHYDROGEN][0]*dense.xIonDense[ipHYDROGEN][0]; */

	bh2dis = hmi.rh2dis*hmi.rel_pop_LTE_H2*Hmolec_old[ipMH]*Hmolec_old[ipMH];	

	if (r->next == NULL) {
		int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH,ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.rh2dis*hmi.rel_pop_LTE_H2;

    

	/* H2 + HNU=>  H2+ + E
	 * photoionization by hard photons, crossection=3*HI */
	gamhd = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s];

	/* cosmic rays predominantly H2 + cr -> H2+ + e, as per table 10 of TH85 */
	/* >>chng 00 nov 28, factor of 0.93 from
	 >>refer	cosmic ray	ionization rate	Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722
	 * also cosmic rays producing secondary ionization csupra */
	/* >>chng 00 nov 28, factor of 0.93 from
	 >>refer	cosmic ray	ionization rate	Maloney, P.R., Hollenbach, D., & Tielens, A. G. G. M., 1998, ApJ, 466, 561
	 */
	/* >>chng 04jan 26, assume ion(H2) = 2x ion(H), H ion rate of
	 * 2.5e-17 s-1, as per
	 * >>refer	cosmicray	ionization	Williams, J.P., Bergin, E.A., Caselli, P., 
	 * >>refercon	Myers, P.C., & Plume, R. 1998, ApJ, 503, 689 
	 * so H2 secondary ionzation rate is 5e-17 s-1 */
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	/* ratio of H2 to H cr rates from table 10 of tielens & hollenbach 1985 */
	
	if (co.lgUMISTrates)
	{
		r->ratek = gamhd + secondaries.csupra[ipHYDROGEN][0]*2.02;
	}

	else
	{
		r->ratek = 4.4e-17;
	}



	/* >>chng 04 apr 22, add H2 + cr -> H+ H + e, TH85 table 10 */
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	/* ratio of H2 to H cr rates from table 10 of tielens & hollenbach 1985 */

	if (co.lgUMISTrates)
	{
		r->ratek = secondaries.csupra[ipHYDROGEN][0]*0.0478;
	}
	else
	{
		r->ratek = 1e-19;
	}
	/* >>chng 03 apr 11 */
	/* H2 + CRP => H + H ; CRP=Cosmic Ray proton 
	 * equation (3643) from
	 * >>refer  Millar, T.J. et.al, 1997,A&AS, 121, 139
	 * h2crphh = 1.3e-18f */  
	h2crphh = 1.3e-18 * hextra.cryden_ov_background;

	 if (!co.lgUMISTrates)
		 h2crphh = 5e-18;

	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2crphh;

	/* >>chng 03 apr 11 */
    /* H2 + CRP => H+ + H_ ; CRP=Cosmic Ray Proton
	 * equation (3644) from
	 * >>refer  T.J.Millar et.al, 1997,A&AS, 121, 139
	 * h2crphphm = 3.9e-21 */
	 h2crphphm = 3.9e-21 * hextra.cryden_ov_background;

	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMHp,ipMHm};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2crphphm*co.lgUMISTrates;

    /* >>chng 03 apr 11 */
    /* H2 + CRP => H+ + H + e ; CRP=Cosmic Ray Proton
	 * equation (3641) from
	 * >>refer  T.J.Millar et.al, 1997,A&AS, 121, 139
	 * h2crphpeh = 2.2e-19f */ 
	 h2crphpeh = 2.2e-19 * hextra.cryden_ov_background;

	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMHp,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2crphpeh*co.lgUMISTrates;

    /* >>chng 03 apr 11 */
    /* H3+ + H_=> H2 + H2 ; 
	 * equation (5,table 9) from
	 * >>refer  Maloney et.al, 1996,ApJ, 466, 561
	 * h3phm2h2 = 1.3e-7f*pow(phycon.te/1000., -0.5) */ 
	
    h3phm2h2 = 1.3e-7 / (phycon.sqrte/31.62278);/*pow(phycon.te/1000., -0.5);*/
	if (r->next == NULL) {
		int in[]={ipMH3p,ipMHm},out[]={ipMH2g,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h3phm2h2*co.lgUMISTrates;

    /* >>chng 03 sep 30 */
    /* H3+ + HNU=> H2+ + H ; 
	 * equation (table 5) from
	 * >>refer	H2	dissoc	Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722  
	 * h3ph2ph = 7.9e-9*hmi.UV_Cont_rel2_Habing_TH85_depth;*/
	/* >>chng 04 jun 13 --  update this rate to match that in the UMIST database */

	if (co.lgUMISTrates)
	{
		h3ph2ph = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_depth/1.66f;
	}

	else
	{
		h3ph2ph = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_face*(float)sexp((2.3*rfield.extin_mag_V_point))/1.66f;
	}

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2p,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =h3ph2ph ;


    /* >>chng 03 sep 30 */
    /* H3+ + HNU=> H2 + H+; 
	 * equation (table 5) from
	 * >>refer	H2	dissoc	Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722  
	 * h3ph2hp = 2.0e-8*hmi.UV_Cont_rel2_Habing_TH85_depth;*/
	/* >>chng 04 jun 13 --  update this rate to match that in the UMIST database */

	if (co.lgUMISTrates)
	{
		h3ph2hp = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_depth/1.66f;
	}

	else
	{
		h3ph2hp = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_face*(float)sexp((1.8*rfield.extin_mag_V_point))/1.66f;
	}

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =h3ph2hp ;


	/* >> chng 02 nov 15 rjrw: multiply c[ipMHo][*] terms by ionization fraction
	 * as b[ipMHo] contains _both_ H0 and H+ */
	/* H2  +  H+  =>  H2+  +  H */

	if (r->next == NULL) {
		int in[]={ipMH2g,ipMHp},out[]={ipMH,ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = rh2h2p*co.lgUMISTrates;

	/* (3,IPMH2P) == destroy H2+ = -sign
	 * H + H2+ => H+ + H2 */
	/* >>chng 02 nov 7 rjrw, remove the destruction rate 
	 *	 c[ipMHo][ipMHo] += -hmi.bh2h2p*hmi.Hmolec[ipMH2p]; 
	 * twice -- reaction changes state of H within single [H0,H+] `species' */

	if (r->next == NULL) {
		int in[]={ipMH,ipMH2p},out[]={ipMHp,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.bh2h2p;

	/*TODO	2	this rate drives numerical instability in such models as secondary1 and 2.in */
	/* this rate couples H2+ and H3+, and tends to destabalize the matrix in both highly
	 * ionized and fully molecular conditions.  Setting this to zero had no effect - the th85
	 * predictions were identical.  
	 *
	 */
	/* H + H3+ => H2 + H2+ */

	h3ph2p = hmrate(2.08e-9,0.,1.88e4);

	if (r->next == NULL) {
		int in[]={ipMH,ipMH3p},out[]={ipMH2g,ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h3ph2p*co.lgUMISTrates;


	/* >>chng 03 feb 7 */
		/* H3+ + H- => H2 + H + H
	 * equation (50) from
	 * >>refer  Stancil, P.C, & Lepp, S, & Dalgarno, A. 1998,ApJ, 509, 1-10
	 * h3phmh2hh = 2.3e-7f*pow(phycon.te/300.0, -0.5) */
	h3phmh2hh = 2.3e-7f*pow(phycon.te/300 , -0.5);
	if (r->next == NULL) {
		int in[]={ipMH3p,ipMHm},out[]={ipMH2g,ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h3phmh2hh*co.lgUMISTrates;

	/* H2 + H3+ => H2 + H2+ + H */
	h3petc = hmrate(3.41e-11,0.5,7.16e4);

	if (r->next == NULL) {
		int in[]={ipMH2g,ipMH3p},out[]={ipMH2g,ipMH2p,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h3petc*co.lgUMISTrates;

	/* H2 + H3+ => H2 + H+ + H2 */
	h32h2 = hmrate(3.41e-11,0.5,5.04e4);

	if (r->next == NULL) {
		int in[]={ipMH2g,ipMH3p},out[]={ipMHp,ipMH2g,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h32h2*co.lgUMISTrates;

	/* e + H3+ => H2 + H */
	/* e + H3+ => 3H was supposed to be included in this rate, 
	 * and stoichiometric factor 2* on sink rate seemed wrong */
	/* >>chng 03 feb 10, increase rate by factor of 13.6 to agree with
	 * >>refer	H3+	DR	McCall, B.J., et al. 2003, Nature, in press (astro-ph 0302106)*/
	/* >>chng 03 feb 13, extra 0.2 since 20% of these go to H2 + H, Stancil private comm */
	/* >>chng 04 apr 22 , update the next two rates to match that of:
	   >>refer H3+ Stancil, P. C., Lepp, S., and Dalgarno, A 509, 1-10; Table 1, reactions #48 and #49 */
	
	eh3_h2h = hmrate(4.00e-8,-0.5,0.)*dense.eden;

	if (!co.lgUMISTrates)
		eh3_h2h = hmrate(2.5e-8,-0.3,0.)*dense.eden;


	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = eh3_h2h;

    /* e + H3+ => 3H */
	eh3p_3h = hmrate(1.6e-7,-0.5,0.)*dense.eden;

	if (!co.lgUMISTrates)
		eh3p_3h = hmrate(7.5e-8,-0.3,0.)*dense.eden;

	/* >>chng 03 feb 10, increase rate by factor of 13.6 to agree with
	 * >>refer	H3+	DR	McCall, B.J., et al. 2003, Nature, in press (astro-ph 0302106)*/
	/* >>chng 03 feb 13, extra 0.8 since 80% of these go to 3H, Stancil private comm */

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH,ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = eh3p_3h;

	if( (trace.lgTrace && trace.lgTr_H2_Mole) || punch.lgPunH2 )
	{
		if( punch.lgPunH2 )
		{
			ioFile = punch.ipPunH2;
		}
		else
		{
			ioFile = ioQQQ;
		}

		if( hmi.H2_rate_destroy > SMALLFLOAT )
		{
			fprintf( ioFile, 
			  " H2 destroy rate=%.2e DIS;%.3f bat;%.3f h2dis;%.3f gamhd;%.3f h2h2p;%.3f E-h;%.3f E-h-;%.2f h2hph3p;%.3f sec;%.3f\n", 
			  hmi.H2_rate_destroy, 
			  hmi.H2_Solomon_dissoc_rate_used / hmi.H2_rate_destroy, 
			  assoc_detach_backwards_grnd / hmi.H2_rate_destroy, 
			  hmi.rh2dis*dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_destroy, 
			  gamhd / hmi.H2_rate_destroy, 
			  rh2h2p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy, 
			  eh2hh /hmi.H2_rate_destroy, 
			  eh2hhm / hmi.H2_rate_destroy ,
			  h2hph3p / hmi.H2_rate_destroy ,
			  secondaries.csupra[ipHYDROGEN][0]*2.02 / hmi.H2_rate_destroy
			  );
		}
		else
		{
			fprintf( ioFile, " Destroy H2: rate=0\n" );
		}
	}

	/*------------------------------------------------------------------- */

	/* h2plus H2+ balance equations */

	/*TODO	2	must add process H2+ + H- => H2 + H, Dalgarno&Lepp 87 */
	/* >>refer	H2+	chemistry	Dalgarno, A., & Lepp, S., 1987, in Astrochemistry, eds. 
	 * >>refercon	M.S. Vardya & S.P. Tarafar, Reidel, Dordrecht, p 109 */
	/* rate = 5e-7 * sqrt(100. / phycon.te); */

	/*TODO	2	put in H2+ + gamma => H + H+ */
	/* >>refer	H2+	chemistry	Stancil, P.C., 1994, ApJ, 430, 360 */
	/* cross section is log10( cs_25) = -1.6547717e6 + 1.8660333e5 ln(nu) - 7.8986431e3*ln(nu)^2
	 * 148.73693 * ln(nu)^3 - 1.0513032*ln(nu)^4 */

	/* make H2+ from Ho
	 * H+  + H  =>  H2+ + HNU
	 * approximation was from Kurucz thesis, not meant for hot gas 
	 * >>chng 02 nov 7 rjrw, stoichiometric factor */
	radath = MAX2(0.,2.325*MIN2(5000.,phycon.te)-1375.)*1e-20;

	if( !co.lgUMISTrates)
		radath = hmrate(5.3e-19,1.85,0);

	if( r->next == NULL) {
		int in[]={ipMH,ipMHp},out[]={ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	/* Printf("O: %ld %g %g\n",rindex,radath*dense.xIonDense[ipHYDROGEN][1]*dense.xIonDense[ipHYDROGEN][0],
		 radath*dense.xIonDense[ipHYDROGEN][1]); */
	rindex++;
	r->ratek = radath;

	/* H2+  +  H+  => H + H+ + H+; Janev et al. 3.2.6 */
	/* >>chng 02 nov 7 rjrw, stoichiometric factor */
	h2pion = 2.4e-27*POW3(phycon.te);

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH,ipMHp},ratesp[]={ipMHp,ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = h2pion*co.lgUMISTrates;

	/* H2+  +  E  => H + H+ + e-; Janev et al. */
	/* >>chng 02 nov 7 rjrw, stoichiometric factor */
	h2pcin = 2e-7*sexp(30720./phycon.te);

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2pcin*dense.eden*co.lgUMISTrates;

	/* >>chng 04 jul 06 -- NPA, include H ionization/recombination due
	to O-H charge transfer in the molecular solver.  */
	/* O + H+ => H + O+ */
	{
		static double oatomic = -1., oion = -1.;
#		if 0
		static long int iter_eval = -2, nzone_eval = -2;
		if( iteration != iter_eval || conv.lgSearch )
		{
			/* new zone or iteration */
			oatomic = dense.xIonDense[ipOXYGEN][0];
			oion = dense.xIonDense[ipOXYGEN][1]; 
			iter_eval = iteration;
		}
		else
		{
			/* same zone and iteration, take mean of old and new abund */
#			define OLD_FRAC	0.5
			oatomic = oatomic*OLD_FRAC + dense.xIonDense[ipOXYGEN][0]*(1.-OLD_FRAC);
			oion = oion*OLD_FRAC + dense.xIonDense[ipOXYGEN][1]*(1.-OLD_FRAC); 
		}
		nzone_eval = nzone;
#		endif
		{
			oatomic = dense.xIonDense[ipOXYGEN][0];
			oion = dense.xIonDense[ipOXYGEN][1]; 
		}

		if (r->next == NULL) {
			int in[]={ipMHp},out[]={ipMH};
			r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
		}

		r = r->next;
		rindex++;
		/*r->ratek = atmdat.HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][0]; */
		r->ratek = atmdat.HCharExcIonOf[ipOXYGEN][0]*oatomic; 

		/* O+ + H => H+ + O */
		if (r->next == NULL) {
			int in[]={ipMH},out[]={ipMHp};
			r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
		}
		r = r->next;
		rindex++;
		/*r->ratek = atmdat.HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1]; */
		r->ratek = atmdat.HCharExcRecTo[ipOXYGEN][0]*oion; 
	}

	/* >>chng 03 aug 22 */
	/* H2+ + e=> H + H ; 
	 * equation (6,table 4) from
	 * >>refer  Maloney et.al, 1996,ApJ, 466, 561
	 * h2pehh = 2.8e-8f*pow(phycon.te/1000., -0.37) */
	/* >>chng 03 sep 01, rm the pow function */
	/*h2pehh = 2.8e-8f*pow(phycon.te/1000., -0.37);*/
	h2pehh = 2.8e-8*12.882/(phycon.te30*phycon.te07);

	if (!co.lgUMISTrates)
		h2pehh = hmrate(1.6e-8,-0.43,0);

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =h2pehh*dense.eden ;

	/* back reaction, H + H+ + e => h2+ + e */
	b2pcin = h2pcin*hmi.rel_pop_LTE_H2p;
	/* this is the hot reaction at high densities */

	if (r->next == NULL) {
		int in[]={ipMH,ipMHp},out[]={ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = b2pcin*dense.eden*co.lgUMISTrates;

	/* H2+  +  HNU  =>  H+  + H */

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = gamtwo;

	/* photoionization by hard photons, crossection =2*HI (wild guess) */

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = 2.*iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*co.lgUMISTrates;

	/* H2 + H2+ => H + H3+ */
	h2ph3p = 1.40e-9*(1. - sexp(9940./phycon.te));

	if (!co.lgUMISTrates)
		h2ph3p = 2.08e-9;


	if (r->next == NULL) {
		int in[]={ipMH2g,ipMH2p},out[]={ipMH,ipMH3p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2ph3p;

	/* destroy H2+ via H2+ + H2 => H + H+ + H2 */
	h2phhp = 2.41e-12*phycon.sqrte*sexp(30720./phycon.te);

	if (r->next == NULL) {
		int in[]={ipMH2g,ipMH2p},out[]={ipMH,ipMHp,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h2phhp*co.lgUMISTrates;

	/*------------------------------------------------------------------ */

	/* H3+ balance equations*/

	/* this reaction does not belong, it is actually the UMIST 
	   reaction CO + H3+ => HCO+ + H2.  Since this is already included below
	   it is not necessary here */

#	if 0

	/* CO + H3+ => H + H2 + CO+ */
	/* >>chng 02 nov 03, add this, which is important for CO dest */
	/* this also appears in hmole.c */
	h3pcop = 1.70e-9*co.hevmol[ipCO];

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = h3pcop;

#	endif
	/* photoionization by hard photons, crossection =2*HI (wild guess)
	 * -- rjrw: where do they go??? 
	 * -- H3+ + hv => H2+ + H+ + e, best guess (P. Stancil, priv comm) */

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2p,ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = 2.*iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*co.lgUMISTrates;

	/*------------------------------------------------------------------ */

	/* vib excited H2, called H2* balance equations, these closely follow
	 * >>refer	mh2	fits	Tielens, A.G.G.M., & Hollenbach, D., 1985a, ApJ 291, 722 */
	/* population of vib-excited H2, from discussion on pp 736-737 of TH85 */

	/* deexcitation rate from upper level, H2* => H2 */
	/* >>chng 03 sep 12, next three expressions did not have factor of 6 -
	 * caught by Gargi Shaw */
	/*H2star_deexcit = Hmolec_old[ipMH2g]*(1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6. +
		 Hmolec_old[ipMH] * (1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.;*/
	deexc_htwo = (1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6.;
	deexc_hneut =  (1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.;

	/* depopulate H2_star, 2e-7 is spontaneous deexcitation rate,
	 * which also appears in lines where intensity of vib lines is entered into line stack */
	/* H2* + H2g -> H2g + H2g */

	if (r->next == NULL) {
		int in[]={ipMH2s},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = 2e-7;

	if (r->next == NULL) {
		int in[]={ipMH2s},out[]={ipMH2g},ratesp[]={ipMH2s,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = deexc_htwo;

	if (r->next == NULL) {
		int in[]={ipMH2s},out[]={ipMH2g},ratesp[]={ipMH2s,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = deexc_hneut;

	/* collisional excitation of vib from ground, 
	 * H2g + H2g -> H2* + H2g, so H2g acts as catalyst
	 * stat weight of ground 1, excit 6, as per TH discussion
	 * this must normally be zero */
	/* H2 producing H2_star */
	/* >>chng 03 sep 11, had been 6, changed to 1 */
	Boltz_fac_H2_H2star = 1.*sexp( 30172./phycon.te);
	/* total excitation rate to H2*, s-1, NB - this is also used in the cooling - heating
	 * rate below */
	/*H2star_excit = Boltz_fac_H2_H2star * H2star_deexcit;*/

	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH2s},ratesp[]={ipMH2g,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = deexc_htwo*Boltz_fac_H2_H2star;

	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH2s},ratesp[]={ipMH,ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
	}
	r = r->next;
	rindex++;
	r->ratek = deexc_hneut*Boltz_fac_H2_H2star;

    /* >>chng 03 aug 28 */
    /* H2* + H => H + H + H
	 * equation  from table 9
	 * >>refer  Tielens and Hollenbach 1985,ApJ, 291, 722-746*/
	 
	if (r->next == NULL) {
		int in[]={ipMH2s,ipMH},out[]={ipMH,ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmrate(9.8e-12,0.5,2.7e4);


    /* >>chng 03 aug 28 */
    /* H2* + H2 => H2 + H + H
	 * equation from table 9
	 * >>refer  Tielens and Hollenbach 1985,ApJ, 291, 722-746*/
	 
	if (r->next == NULL) {
		int in[]={ipMH2s,ipMH2g},out[]={ipMH2g,ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmrate(9.8e-12,0.5,2.7e4);

    /* >>chng 03 aug 28 */
    /* H2* + H2* => H2 + H + H
	 * equation from table 9
	 * >>refer  Tielens and Hollenbach 1985,ApJ, 291, 722-746*/
	if (r->next == NULL) {
		int in[]={ipMH2s,ipMH2s},out[]={ipMH2g,ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmrate(9.8e-12,0.5,0.);

	/* >>03 may 22, hmi.H2_Solomon_dissoc_rate_used is now rate electronic excited states
	 * decay into X continuum, not the rate that elec excited states are excited.
	 * assuming that 10% of excitations lead to ionization, this rate, for relax
	 * into ro-vib excited states, is 10x the Solomon rate */
	/* assume that 0.9 of H2 dissociations lead to H2_star,
	 * H2 + 0.9*hmi.H2_Solomon_dissoc_rate_used => h2_star */
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH2s};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	/* >>chng 03 may 22, see comment just above */
	/*r->ratek = 0.9*hmi.H2_Solomon_dissoc_rate_used;
	r->ratek = 10.*hmi.H2_Solomon_dissoc_rate_used;*/
	/* >>chng 03 sep 11, use real evaluated rate */
	r->ratek = hmi.H2_H2g_to_H2s_rate_used;

	/* rate of photodissoc of vib-excit H2, A12 of TH85 */
	if (r->next == NULL) {
		int in[]={ipMH2s},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hmi.H2_photodissoc_used;

	/*---------------------------------------------------------------- */

	/* He H+ formation rates taken from Flower+Roueff, Black */

	/* He+ + H => HeH+
	 * radiative association from 
	 * >>refer	mheh+	rate	Zygelman, B., and Dalgarno, A. 1990, ApJ 365, 239 */

	if (r->next == NULL) {
		int in[]={ipMH},out[]={ipMHeHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = 1e-15*dense.xIonDense[ipHELIUM][1];


#	if 0
/* >>chng 04 jun 09, is also in iso_recom, then in ion solver, do not double count */
    /* >>chng 03 sep 30 */
    /* He+ + H=> He + H+; 
	 * equation (table 6) from
	 * >>refer	H2	dissoc	Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722  
	 * hephhphe = 1.9e-15*dense.xIonDense[ipHELIUM][1];*/
	/* >>chng 04 jun 10 -- put in the UMIST rate for this reaction */
	hmi.hephhehp = (float)hmrate(4.85e-15, 0.18, 0.);
    hephhphe = hmi.hephhehp*dense.xIonDense[ipHELIUM][1];
	if (r->next == NULL) {
		int in[]={ipMH},out[]={ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = hephhphe;

#	endif

    /* >>chng 03 sep 30 */
	/* back reaction of He+ + H=> He + H+	*/
    /* He + H+=> He+ + H;   
	 * bhephhphe = hephhphe*dense.xIonDense[ipHELIUM][1]*dense.xIonDense[ipHYDROGEN][0]/dense.xIonDense[ipHYDROGEN][1];*/
    /*bhephhphe = hephhphe*dense.xIonDense[ipHELIUM][1]*dense.xIonDense[ipHYDROGEN][0]/SDIV(dense.xIonDense[ipHYDROGEN][1]);
	if (r->next == NULL) {
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = bhephhphe;*/

    /* >>chng 03 sep 30 */
    /* He + H-=> He + H + e; 
	 * equation (table 5) from
	 * >>refer Paolo Lenzuni, David F. Chernoff, Edwin E. Salpeter, 1991, ApJS, 76, 759L  
	 * hehmeheh = 4.1e-17*pow(phycon.te,2)*sexp(19870/phycon.te)*dense.xIonDense[ipHELIUM][0];*/
	hehmeheh = 4.1e-17*phycon.tesqrd*sexp(19870/phycon.te)*dense.xIonDense[ipHELIUM][0];
	if (r->next == NULL) {
		int in[]={ipMHm},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =hehmeheh  ;

    /* >>chng 03 sep 30 */
    /* He+ + H2=> He + H+ + H; 
	 * equation (table 6) from
	 * >>refer	H2	dissoc	Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722  
	 * heph2hpheh = 1.5e-13*dense.xIonDense[ipHELIUM][1];*/
	/*hmi.rheph2hpheh = 1.5e-13f;*/
    /* >>chng 04 jun 10 */
	/* use the rate for this reaction as defined in UMIST */
	hmi.rheph2hpheh = (float)hmrate(3.7e-14, 0., 35);
    heph2hpheh = hmi.rheph2hpheh*dense.xIonDense[ipHELIUM][1];
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMHp,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = heph2hpheh;

	/* >>chng 04 jun 10 -- This reaction was not orignally included.  It is in
	    the UMIST database and seems to be an important He+ destruction mechanism */
    /* He+ + H2=> He + H2+  
	 * use the rate for this reaction as defined in UMIST */
	/*hmi.heph2heh2p = (float)hmrate(3.7e-14, 0., 35);*/
	hmi.heph2heh2p = (float)hmrate(7.2e-15, 0., 0);
    heph2heh2p = hmi.heph2heh2p*dense.xIonDense[ipHELIUM][1];
	if (r->next == NULL) {
		int in[]={ipMH2g},out[]={ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = heph2heh2p;


	/* He + H+ => HeH+ */	
	if (r->next == NULL) {
		int in[]={ipMHp},out[]={ipMHeHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = 1e-20*dense.xIonDense[ipHELIUM][0];

	/* H2+ + HE => HEH+ + H0 */
	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH,ipMHeHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = 3e-10*exp(-6717./phycon.te)*dense.xIonDense[ipHELIUM][0]*co.lgUMISTrates;

	/* photodissociation through 1.6->2.3 continuum */

	/* why is this in a look instead of GammaK?
	 * to fix must set opacities into stack */
	gamheh = 0.;
	limit = MIN2(hmi.iheh2-1 , rfield.nflux );
	for( i=hmi.iheh1-1; i < limit; i++ )
	{
		gamheh += rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i];
	}
	gamheh *= 4e-18;

	/* hard radiation */
	gamheh += 3.*iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s];

	/* recombination, HeH+  +  e => He + H */
	gamheh += dense.eden*1e-9;

	if (r->next == NULL) {
		int in[]={ipMHeHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = gamheh;

	/* HeH+  +  H => H2+  + He */
	if (r->next == NULL) {
		int in[]={ipMH,ipMHeHp},out[]={ipMH2p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = 1e-10*co.lgUMISTrates;

    /* >>chng 03 sep 30 */
    /* HeH+  +  H2 => H3+  + He 
	 * equation (He13) from
	 * >>refer D.Galli & F.Palla, 1998, A&A,335,403-420 
	 * hehph2h3phe = 1.3e-9;*/
	hehph2h3phe = 1.3e-9;
	if (r->next == NULL) {
		int in[]={ipMH2g,ipMHeHp},out[]={ipMH3p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =hehph2h3phe*co.lgUMISTrates;

    /* >>chng 03 sep 30 */
    /* He+  +  H- => H  + He 
	 * equation (20) from
	 * >>refer P.C.Stancil, S.Lepp, A.Dalgarno,1998, ApJ,509,1-10
	 * hephmhhe = 2.32e-7*pow(phycon.te/300,-.52)*sexp(phycon.te/-22400.)*dense.xIonDense[ipHELIUM][1];*/
	/* >>chng 03 oct 22, only add for Te < 1e5 - caused exception as T -> inf */
	hephmhhe = 2.32e-7*pow(phycon.te/300,-.52)*exp(TStancil/22400.)*dense.xIonDense[ipHELIUM][1];
	if (r->next == NULL) {
		int in[]={ipMHm},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =hephmhhe;


	/* >>chng 04 may 26, NA add this code */
	/* Many reactions in co.c also will create or destroy molecules that are predicted 
	   in hmole.  These reactions need to be inserted into this solver to obtain a more
	   accurate solution. These are the rates for reactions in co.c that create or destroy
	   species that are predicted in hmole_step.c.*/
		
	co.H_CH_C_H_H				= hmrate(0.000000006,0,40200)*co.hevmol[ipCH];
	co.H_OH_O_H_H				= hmrate(0.000000006,0,50900)*co.hevmol[ipOH];
	co.H_H2O_OH_H_H				= hmrate(0.0000000058,0,52900)*co.hevmol[ipH2O];
	co.H_COP_CO_HP				= hmrate(0.00000000075,0,0)*co.hevmol[ipCOP];
	co.H_CH_C_H2				= hmrate(0.000000000027,0.38,0)*co.hevmol[ipCH];
	co.H_CHP_CP_H2				= hmrate(0.00000000075,0,0)*co.hevmol[ipCHP];
	co.H_CH2_CH_H2				= hmrate(0.0000000000664,0,0)*co.hevmol[ipCH2];
	co.H_CH3P_CH2P_H2			= hmrate(0.0000000007,0,10560)*co.hevmol[ipCH3P];
	co.H_OH_O_H2				= hmrate(0.0000000000000699,2.8,1950)*co.hevmol[ipOH];
	co.H_H2O_OH_H2				= hmrate(0.0000000000159,1.2,9610)*co.hevmol[ipH2O];
	co.Hminus_HCOP_CO_H2	   = hmrate(0.00000023,-0.5,0)*co.hevmol[ipHCOP];
	co.Hminus_H3OP_H2O_H2	  = hmrate(0.00000023,-0.5,0)*co.hevmol[ipH3OP];
	co.Hminus_H3OP_OH_H2_H	 = hmrate(0.00000023,-0.5,0)*co.hevmol[ipH3OP];
	co.HP_CH_CHP_H			    = hmrate(0.0000000019,0,0)*co.hevmol[ipCH];
	co.HP_CH2_CH2P_H		    = hmrate(0.0000000014,0,0)*co.hevmol[ipCH2];
	co.HP_H2O_H2OP_H		    = hmrate(0.0000000069,0,0)*co.hevmol[ipH2O];
	co.HP_O2_O2P_H			    = hmrate(0.000000002,0,0) *co.hevmol[ipO2];
	co.HP_OH_OHP_H			    = hmrate(0.0000000021,0,0)*co.hevmol[ipOH];
	co.HP_SiO_SiOP_H		    = hmrate(0.0000000033,0,0)*co.hevmol[ipSIO];
	co.HP_CH2_CHP_H2		    = hmrate(0.0000000014,0,0)*co.hevmol[ipCH2];
	co.HP_SiH_SiP_H2		    = hmrate(0.0000000017,0,0)*co.hevmol[ipSIH];
	co.H2_CHP_CH2P_H		    = hmrate(0.0000000012,0,0)*co.hevmol[ipCHP];
	co.H2_CH2P_CH3P_H		    = hmrate(0.0000000016,0,0)*co.hevmol[ipCH2P];
	co.H2_OHP_H2OP_H		    = hmrate(0.00000000101,0,0)*co.hevmol[ipOHP];
	co.H2_H2OP_H3OP_H		    = hmrate(0.00000000064,0,0)*co.hevmol[ipH2OP];
	co.H2_COP_HCOP_H		    = hmrate(0.0000000018,0,0)*co.hevmol[ipCOP];
	co.H2_OP_OHP_H			    = hmrate(0.0000000017,0,0)*co.hevmol[ipOP];
	co.H2_SiOP_SiOHP_H		    = hmrate(0.00000000032,0,0)*co.hevmol[ipSIOP];
	co.H2_C_CH_H				= hmrate(0.000000000664,0,11700)*co.hevmol[ipATC];
	co.H2_CP_CHP_H				= hmrate(0.0000000001,0,4640)*co.hevmol[ipCP];
	co.H2_CH_CH2_H				= hmrate(0.000000000546,0,1943)*co.hevmol[ipCH];
	co.H2_OH_H2O_H				= hmrate(0.00000000000205,1.52,1736)*co.hevmol[ipOH];
	co.H2_O_OH_H				= hmrate(0.000000000000314,2.7,3150)*co.hevmol[ipATO];
	co.H2_CH_C_H2_H				= hmrate(0.000000006,0,40200)*co.hevmol[ipCH];
	co.H2_OH_O_H2_H				= hmrate(0.000000006,0,50900)*co.hevmol[ipOH];
	co.H2_H2O_OH_H2_H			= hmrate(0.0000000058,0,52900)*co.hevmol[ipH2O];
	co.H2_O2_O_O_H2				= hmrate(0.000000006,0,52300)*co.hevmol[ipO2];
	co.H2_C_CH_H				= hmrate(0.000000000664,0,11700)*co.hevmol[ipATC];
	co.H2s_CH_C_H2_H			= hmrate(0.000000006,0,0)*co.hevmol[ipCH]*hmi.lgLeiden_Keep_ipMH2s ;
	co.H2s_OH_O_H2_H			= hmrate(0.000000006,0,0)*co.hevmol[ipOH]*hmi.lgLeiden_Keep_ipMH2s ;
	co.H2s_H2O_OH_H2_H			= hmrate(0.0000000058,0,0)*co.hevmol[ipH2O]*hmi.lgLeiden_Keep_ipMH2s ;
	co.H2s_O2_O_O_H2			= hmrate(0.000000006,0,0)*co.hevmol[ipO2]*hmi.lgLeiden_Keep_ipMH2s ;
	co.H2P_C_CHP_H				 = hmrate(0.0000000024,0,0)*co.hevmol[ipATC];
	co.H2P_CH_CH2P_H			    = hmrate(0.00000000071,0,0)*co.hevmol[ipCH];
	co.H2P_CH2_CH3P_H			= hmrate(0.000000001,0,0)*co.hevmol[ipCH2];
	co.H2P_OH_H2OP_H			= hmrate(0.00000000076,0,0)*co.hevmol[ipOH];
	co.H2P_H2O_H3OP_H			= hmrate(0.0000000034,0,0)*co.hevmol[ipH2O];
	co.H2P_CO_HCOP_H			= hmrate(0.00000000216,0,0)*co.hevmol[ipCO];
	co.H2P_O_OHP_H				= hmrate(0.0000000015,0,0)*co.hevmol[ipATO];
	co.H2P_CH_CHP_H2			= hmrate(0.00000000071,0,0)*co.hevmol[ipCH];
	co.H2P_CH2_CH2P_H2			= hmrate(0.000000001,0,0)*co.hevmol[ipCH2];
	co.H2P_CO_COP_H2			= hmrate(0.00000000064,0,0)*co.hevmol[ipCO];
	co.H2P_H2O_H2OP_H2			= hmrate(0.0000000039,0,0)*co.hevmol[ipH2O];
	co.H2P_O2_O2P_H2			= hmrate(0.0000000008,0,0)*co.hevmol[ipO2];
	co.H2P_OH_OHP_H2			= hmrate(0.00000000076,0,0)*co.hevmol[ipOH];
	co.H3P_C_CHP_H2				= hmrate(0.000000002,0,0)*co.hevmol[ipATC];
	co.H3P_CH_CH2P_H2			= hmrate(0.0000000012,0,0)*co.hevmol[ipCH];
	co.H3P_CH2_CH3P_H2			= hmrate(0.0000000017,0,0)*co.hevmol[ipCH2];
	co.H3P_OH_H2OP_H2			= hmrate(0.0000000013,0,0)*co.hevmol[ipOH];
	co.H3P_H2O_H3OP_H2			= hmrate(0.0000000059,0,0)*co.hevmol[ipH2O];
	co.H3P_CO_HCOP_H2			= hmrate(0.0000000017,0,0)*co.hevmol[ipCO];
	co.H3P_O_OHP_H2				= hmrate(0.0000000008,0,0)*co.hevmol[ipATO];
	co.H3P_SiH_SiH2P_H2			= hmrate(0.000000002,0,0)*co.hevmol[ipSIH];
	co.H3P_SiO_SiOHP_H2			= hmrate(0.000000002,0,0)*co.hevmol[ipSIO];
	co.h2s_ch_ch2_h				= hmrate(0.000000000546,0, 0. )*co.hevmol[ipCH]*hmi.lgLeiden_Keep_ipMH2s ;
	co.h2s_o_oh_h				= hmrate(0.000000000000314,2.7, 0. )*co.hevmol[ipATO]*hmi.lgLeiden_Keep_ipMH2s ;
	co.h2s_oh_h2o_h				= hmrate(0.00000000000205,1.52, 0. )*co.hevmol[ipOH]*hmi.lgLeiden_Keep_ipMH2s ;
	co.h2s_c_ch_h				= hmrate(0.000000000664,0, 0. )*co.hevmol[ipATC]*hmi.lgLeiden_Keep_ipMH2s ;
	co.h2s_cp_chp_h				= hmrate(0.0000000001,0, 0. )*co.hevmol[ipCP]*hmi.lgLeiden_Keep_ipMH2s ;
	co.H_CH3_CH2_H2				= hmrate(0.0000000001,0,7600)*co.hevmol[ipCH3];
	co.H_CH4P_CH3P_H2			= hmrate(0.00000000001,0,0)*co.hevmol[ipCH4P];
	co.H_CH5P_CH4P_H2			= hmrate(0.00000000002,0,0)*co.hevmol[ipCH5P];
	co.H2_CH2_CH3_H				= hmrate(0.0000000000518,0.17,6400)*co.hevmol[ipCH2];
	co.H2_CH3_CH4_H				= hmrate(0.0000000000000686,2.74,4740)*co.hevmol[ipCH3];
	co.H2_CH4P_CH5P_H			= hmrate(0.000000000033,0,0)*co.hevmol[ipCH4P];
	co.H2s_CH2_CH3_H			= hmrate(0.0000000000518,0.17,0)*co.hevmol[ipCH2]*hmi.lgLeiden_Keep_ipMH2s ;
	co.H2s_CH3_CH4_H			= hmrate(0.0000000000000686,2.74,0)*co.hevmol[ipCH3]*hmi.lgLeiden_Keep_ipMH2s ;
	co.H2P_CH4_CH3P_H2			= hmrate(0.0000000023,0,0)*co.hevmol[ipCH4];
	co.H2P_CH4_CH4P_H2			= hmrate(0.0000000014,0,0)*co.hevmol[ipCH4];
	co.H2P_CH4_CH5P_H			= hmrate(0.000000000114,0,0)*co.hevmol[ipCH4];
	co.H3P_CH3_CH4P_H2			= hmrate(0.0000000021,0,0)*co.hevmol[ipCH];
	co.H3P_CH4_CH5P_H2			= hmrate(0.0000000024,0,0)*co.hevmol[ipCH4];
	co.HP_CH3_CH3P_H			= hmrate(0.0000000034,0,0)*co.hevmol[ipCH3];
	co.HP_CH4_CH3P_H2			= hmrate(0.0000000023,0,0)*co.hevmol[ipCH4];
	co.HP_CH4_CH4P_H			= hmrate(0.0000000015,0,0)*co.hevmol[ipCH4];
	

    if (r->next == NULL) 
	{
		int in[]={ipMH},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek = co.H_CH_C_H_H;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_OH_O_H_H;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMH,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_H2O_OH_H_H;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_COP_CO_HP;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_CH_C_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_CHP_CP_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_CH2_CH_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_CH3P_CH2P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_OH_O_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_H2O_OH_H2;

	
	if (r->next == NULL) 
	{
		int in[]={ipMHm},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.Hminus_HCOP_CO_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMHm},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.Hminus_H3OP_H2O_H2;

	if (r->next == NULL)
	{
		int in[]={ipMHm},out[]={ipMH2g, ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.Hminus_H3OP_OH_H2_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_CH_CHP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_CH2_CH2P_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_H2O_H2OP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_O2_O2P_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_OH_OHP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_SiO_SiOP_H;

	if (r->next == NULL) 
	{
		int in[]={ipMHp},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_CH2_CHP_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMHp},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.HP_SiH_SiP_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_CHP_CH2P_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_CH2P_CH3P_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_OHP_H2OP_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_H2OP_H3OP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_COP_HCOP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_OP_OHP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_SiOP_SiOHP_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_C_CH_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_CP_CHP_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_CH_CH2_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_OH_H2O_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_O_OH_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH2g,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_CH_C_H2_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2g},out[]={ipMH2g, ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_OH_O_H2_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH2g, ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_H2O_OH_H2_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2g},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_O2_O_O_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_C_CH_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2s},out[]={ipMH2g,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2s_CH_C_H2_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2s},out[]={ipMH2g,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2s_OH_O_H2_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2s},out[]={ipMH2g,ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2s_H2O_OH_H2_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2s},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2s_O2_O_O_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_C_CHP_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_CH_CH2P_H;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_CH2_CH3P_H;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_OH_H2OP_H;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_H2O_H3OP_H;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_CO_HCOP_H;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_O_OHP_H;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_CH_CHP_H2;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_CH2_CH2P_H2;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_CO_COP_H2;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_H2O_H2OP_H2;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_O2_O2P_H2;

	if (r->next == NULL) {
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2P_OH_OHP_H2;

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_C_CHP_H2;

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_CH_CH2P_H2;

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_CH2_CH3P_H2;

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_OH_H2OP_H2;

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_H2O_H3OP_H2;

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_CO_HCOP_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_O_OHP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_SiH_SiH2P_H2;

	if (r->next == NULL) {
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H3P_SiO_SiOHP_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH2s},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.h2s_ch_ch2_h;

	if (r->next == NULL)
	{
		int in[]={ipMH2s},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.h2s_o_oh_h;

	if (r->next == NULL) 
	{
		int in[]={ipMH2s},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.h2s_oh_h2o_h;


	if (r->next == NULL) 
	{
		int in[]={ipMH2s},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.h2s_c_ch_h;

	if (r->next == NULL)
	{
		int in[]={ipMH2s},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.h2s_cp_chp_h;

	if (r->next == NULL) 
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_CH3_CH2_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_CH4P_CH3P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H_CH5P_CH4P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2_CH2_CH3_H;

	
	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_CH3_CH4_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_CH4P_CH5P_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2s},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =co.H2s_CH2_CH3_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2s},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2s_CH3_CH4_H;

	if (r->next == NULL) 
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_CH4_CH3P_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_CH4_CH4P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_CH4_CH5P_H;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_CH3_CH4P_H2;

	if (r->next == NULL) 
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_CH4_CH5P_H2;

	if (r->next == NULL) {
		int in[]={ipMHp},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_CH4_CH3P_H2;

	if (r->next == NULL) {
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_CH4_CH4P_H;

	/*  >>chng 04 jul 14 -- NPA. This is the contribution that the heavy element 
	 *  network makes to the hydrogen chemistry due to reactions with N and S. */

		co.H2_N_NH_H           = hmrate(0.00000000169,0,18095)*co.hevmol[ipATN];
		co.H2_NH_NH2_H         = hmrate(0.0000000000596,0,7782)*co.hevmol[ipNH];
		co.H2_NH2_NH3_H        = hmrate(0.00000000000000205,3.89,1400)*co.hevmol[ipNH2];
		co.H2_CN_HCN_H         = hmrate(0.000000000000404,2.87,820)*co.hevmol[ipCN];
		co.HP_HNO_NOP_H2       = hmrate(0.000000004,0,0)*co.hevmol[ipHNO];
		co.HP_HS_SP_H2         = hmrate(0.0000000016,0,0)*co.hevmol[ipHS];
		co.H_HSP_SP_H2         = hmrate(0.00000000011,0,0)*co.hevmol[ipHSP];
		co.H2P_N_NHP_H         = hmrate(0.0000000019,0,0)*co.hevmol[ipATN];
		co.H2_NP_NHP_H         = hmrate(0.000000001,0,85)*co.hevmol[ipNP];
		co.H2_NHP_N_H3P        = hmrate(0.000000000225,0,0)*co.hevmol[ipNHP];
		co.H2P_NH_NH2P_H       = hmrate(0.00000000076,0,0)*co.hevmol[ipNH];
		co.H2_NHP_NH2P_H       = hmrate(0.00000000128,0,0)*co.hevmol[ipNHP];
		co.H2_NH2P_NH3P_H      = hmrate(0.00000000027,0,0)*co.hevmol[ipNH2P];
		co.H2_NH3P_NH4P_H      = hmrate(0.0000000000002,0,0)*co.hevmol[ipNH3P];
		co.H2P_CN_HCNP_H       = hmrate(0.0000000012,0,0)*co.hevmol[ipCN];
		co.H2_CNP_HCNP_H       = hmrate(0.000000001,0,0)*co.hevmol[ipCNP];
		co.H2P_NO_HNOP_H       = hmrate(0.0000000011,0,0)*co.hevmol[ipNO];
		co.H2_SP_HSP_H         = hmrate(0.00000000011,0,9860)*co.hevmol[ipSP];
		co.H2_CSP_HCSP_H       = hmrate(0.00000000045,0,0)*co.hevmol[ipCSP];
		co.H3P_NH_NH2P_H2      = hmrate(0.0000000013,0,0)*co.hevmol[ipNH];
		co.H3P_NH2_NH3P_H2     = hmrate(0.0000000018,0,0)*co.hevmol[ipNH2];
		co.H3P_NH3_NH4P_H2     = hmrate(0.0000000027,0,0)*co.hevmol[ipNH3];
		co.H3P_CN_HCNP_H2      = hmrate(0.000000002,0,0)*co.hevmol[ipCN];
		co.H3P_NO_HNOP_H2      = hmrate(0.0000000011,0,0)*co.hevmol[ipNO];
		co.H3P_S_HSP_H2        = hmrate(0.0000000026,0,0)*co.hevmol[ipATS];
		co.H3P_CS_HCSP_H2      = hmrate(0.0000000029,0,0)*co.hevmol[ipCS];
		co.H3P_NO2_NOP_OH_H2   = hmrate(0.0000000007,0,0)*co.hevmol[ipNO2];
		co.HP_NH_NHP_H         = hmrate(0.0000000021,0,0)*co.hevmol[ipNH];
		co.HP_NH2_NH2P_H       = hmrate(0.0000000029,0,0)*co.hevmol[ipNH2];
		co.HP_NH3_NH3P_H       = hmrate(0.0000000011,0,0)*co.hevmol[ipNH3];
		co.H_CNP_CN_HP         = hmrate(0.00000000019,0,0)*co.hevmol[ipCNP];
		co.HP_HCN_HCNP_H       = hmrate(0.0000000105,-0.13,0)*co.hevmol[ipHCN];
		co.H_HCNP_HCN_HP       = hmrate(0.000000000037,0,0)*co.hevmol[ipHCNP];
		co.H_N2P_N2_HP         = hmrate(0.00000000012,0,0)*co.hevmol[ipN2P];
		co.HP_NO_NOP_H         = hmrate(0.0000000029,0,0)*co.hevmol[ipNO];
		co.HP_HS_HSP_H         = hmrate(0.0000000016,0,0)*co.hevmol[ipHS];
		co.HP_SIN_SINP_H       = hmrate(0.000000003,0,0)*co.hevmol[ipSIN];
		co.HP_CS_CSP_H         = hmrate(0.0000000049,0,0)*co.hevmol[ipCS];
		co.HP_NS_NSP_H         = hmrate(0.0000000047,0,0)*co.hevmol[ipNS];
		co.HP_SO_SOP_H         = hmrate(0.0000000032,0,0)*co.hevmol[ipSO];
		co.HP_OCS_OCSP_H       = hmrate(0.0000000021,0,0)*co.hevmol[ipOCS];
		co.HP_S2_S2P_H         = hmrate(0.000000003,0,0)*co.hevmol[ipS2];
		co.H2P_NH_NHP_H2       = hmrate(0.00000000076,0,0)*co.hevmol[ipNH];
		co.H2P_NH2_NH2P_H2     = hmrate(0.0000000021,0,0)*co.hevmol[ipNH2];
		co.H2P_NH3_NH3P_H2     = hmrate(0.0000000057,0,0)*co.hevmol[ipNH3];
		co.H2P_CN_CNP_H2       = hmrate(0.0000000012,0,0)*co.hevmol[ipCN];
		co.H2P_HCN_HCNP_H2     = hmrate(0.0000000027,0,0)*co.hevmol[ipHCN];
		co.H2P_NO_NOP_H2       = hmrate(0.0000000011,0,0)*co.hevmol[ipNO];
	
	/* The following terms are for the reactions with the rates defined above.  The 
	 * actual reaction is self-contained in the rate label. */

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_N_NH_H ;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_NH_NH2_H ;

	
	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_NH2_NH3_H ;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_CN_HCN_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_HNO_NOP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_HS_SP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H_HSP_SP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_N_NHP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_NP_NHP_H ;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH3p};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_NHP_N_H3P ;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_NH_NH2P_H;

	
	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_NHP_NH2P_H ;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_NH2P_NH3P_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_NH3P_NH4P_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_CN_HCNP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_CNP_HCNP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_NO_HNOP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_SP_HSP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2g},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2_CSP_HCSP_H;


	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_NH_NH2P_H2;

	
	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_NH2_NH3P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_NH3_NH4P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_CN_HCNP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_NO_HNOP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_S_HSP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_CS_HCSP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH3p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H3P_NO2_NOP_OH_H2;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_NH_NHP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_NH2_NH2P_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_NH3_NH3P_H;

	
	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H_CNP_CN_HP;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_HCN_HCNP_H;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H_HCNP_HCN_HP;

	if (r->next == NULL)
	{
		int in[]={ipMH},out[]={ipMHp};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H_N2P_N2_HP;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_NO_NOP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_HS_HSP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_SIN_SINP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_CS_CSP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_NS_NSP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_SO_SOP_H;

	
	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_OCS_OCSP_H;

	if (r->next == NULL)
	{
		int in[]={ipMHp},out[]={ipMH};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.HP_S2_S2P_H;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_NH_NHP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_NH2_NH2P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_NH3_NH3P_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_CN_CNP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_HCN_HCNP_H2;

	if (r->next == NULL)
	{
		int in[]={ipMH2p},out[]={ipMH2g};
		r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
	}
	r = r->next;
	rindex++;
	r->ratek =	co.H2P_NO_NOP_H2;


	/*  The following reactions in co.c include a formation or destruction
	    process, but not both, for one of the molecules in hmole_step.c.
		Because of this, the following rates have to be put into the bvec[i]
		part of the matrix. */

	/* Zero out the CO contribution to the sources and sinks in hmole */

	
	for (i=0; i<N_H_MOLEC; ++i)
	{
		co.hydro_source[i] = 0;
		co.hydro_sink[i] = 0;
	}

	/* C + H3OP =  HCOP + H2 */   
	
	co.C_H3OP_HCOP_H2_1 =  hmrate(1.0e-11,0,0)*co.hevmol[ipH3OP]*co.hevmol[ipATC];
	co.hydro_source[ipMH2g] += (float)co.C_H3OP_HCOP_H2_1 ;


	/* C + OH =  CO + H  */  

	co.C_OH_CO_H_1 = hmrate(1.1e-10,0.5,0)*co.hevmol[ipOH]*co.hevmol[ipATC];
	co.hydro_source[ipMH] += (float)co.C_OH_CO_H_1 ;
	
	/* CP + OH =>  CO + HP */   

	co.CP_OH_CO_HP_1 = hmrate(7.7e-10,0,0)*co.hevmol[ipOH]*co.hevmol[ipCP];
	co.hydro_source[ipMHp] += (float)co.CP_OH_CO_HP_1 ;

	/* CP + H2O =>  HCOP + H  */  
	
	co.CP_H2O_HCOP_H_1 = hmrate(9.0e-10,0,0)*co.hevmol[ipH2O]*co.hevmol[ipCP];
	co.hydro_source[ipMH] += (float)co.CP_H2O_HCOP_H_1 ;
	
	/* CP + OH =>  COP + H  */  
	
	co.CP_OH_COP_H_1 = hmrate(7.7e-10,0,0)*co.hevmol[ipOH]*co.hevmol[ipCP];
	co.hydro_source[ipMH] += (float)co.CP_OH_COP_H_1 ;
	
	/* O + CH =>  CO + H  */  
	
	co.O_CH_CO_H_1 = hmrate(6.6e-11,0,0)*co.hevmol[ipCH]*co.hevmol[ipATO];
	co.hydro_source[ipMH] += (float)co.O_CH_CO_H_1 ;
	
	/* O + CHP =>  COP + H  */  
	
	co.O_CHP_COP_H_1 = hmrate(3.5e-10,0,0)*co.hevmol[ipCHP]*co.hevmol[ipATO];
	co.hydro_source[ipMH] += (float)co.O_CHP_COP_H_1 ;
	
	/* O + CH2 =>  CO + H + H  */  
	
	co.O_CH2_CO_H_H_1 = hmrate(1.33e-10,0,0)*co.hevmol[ipCH2]*co.hevmol[ipATO];
	co.hydro_source[ipMH] += (float)co.O_CH2_CO_H_H_1 ;
	
	/* O + CH2 =>  CO + H2 */  
	
	co.O_CH2_CO_H2_1 = hmrate(8.0e-11,0,0)*co.hevmol[ipCH2]*co.hevmol[ipATO];
	co.hydro_source[ipMH2g] += (float)co.O_CH2_CO_H2_1 ;
	
	/* O + CH2P =>  HCOP + H  */  
	
	co.O_CH2P_HCOP_H_1 = hmrate(7.5e-10,0,0)*co.hevmol[ipCH2P]*co.hevmol[ipATO];
	co.hydro_source[ipMH] += (float)co.O_CH2P_HCOP_H_1 ;
	
	/* O + CH3P =>  HCOP + H2  */   
	
	co.O_CH3P_HCOP_H2_1	= hmrate(4.0e-10,0,0)*co.hevmol[ipCH3P]*co.hevmol[ipATO];
	co.hydro_source[ipMH2g] += (float)co.O_CH3P_HCOP_H2_1 ;
	
	/* O + H2OP =>  O2P + H2  */   
	
	co.O_H2OP_O2P_H2_1 = hmrate(4.0e-11,0,0)*co.hevmol[ipH2OP]*co.hevmol[ipATO];
	co.hydro_source[ipMH2g] += (float)co.O_H2OP_O2P_H2_1 ;
	
	/* O + OH =>  O2 + H  */  
	
	co.O_OH_O2_H_1 = hmrate(4.34e-11,-0.5,30)*co.hevmol[ipOH]*co.hevmol[ipATO];
	co.hydro_source[ipMH] += (float)co.O_OH_O2_H_1 ;
	
	/* O + OHP =>  O2P + H  */  
	
	co.O_OHP_O2P_H_1 = hmrate(7.1e-10,0,0)*co.hevmol[ipOHP]*co.hevmol[ipATO];
	co.hydro_source[ipMH] += (float)co.O_OHP_O2P_H_1 ;
	
	/* O + SiH =>  SiO + H  */  
	
	co.O_SiH_SiO_H_1 = hmrate(4.0e-11,0.5,0)*co.hevmol[ipSIH]*co.hevmol[ipATO];
	co.hydro_source[ipMH] += (float)co.O_SiH_SiO_H_1 ;
	
	/* O + SiH2P =>  SiOHP + H  */  
	
	co.O_SiH2P_SiOHP_H_1 = hmrate(6.3e-10,0,0)*co.hevmol[ipSIH2P]*co.hevmol[ipATO];	
	co.hydro_source[ipMH] += (float)co.O_SiH2P_SiOHP_H_1 ;
	
	/* OP + CH =>  COP + H  */  
	
	co.OP_CH_COP_H_1 = hmrate(3.5e-10,0,0)*co.hevmol[ipCH]*co.hevmol[ipOP];
	co.hydro_source[ipMH] += (float)co.OP_CH_COP_H_1 ;
	
	/* OP + OH =>  O2P + H  */  
	
	co.OP_OH_O2P_H_1 = hmrate(3.6e-10,0,0)*co.hevmol[ipOH]*co.hevmol[ipOP];
	co.hydro_source[ipMH] += (float)co.OP_OH_O2P_H_1 ;
	
	/* Si + OH =>  SiO + H  */  
	
	co.Si_OH_SiO_H_1 = hmrate(2.0e-10,0.5,0)*co.hevmol[ipOH]*co.hevmol[ipATSI];
	co.hydro_source[ipMH] += (float)co.Si_OH_SiO_H_1 ;
	
	/* SiP + H2O =>  SiOHP + H  */  
	
	co.SiP_H2O_SiOHP_H_1 = hmrate(2.3e-10,0,0)*co.hevmol[ipH2O]*co.hevmol[ipSIP];			
	co.hydro_source[ipMH] += (float)co.SiP_H2O_SiOHP_H_1 ;
	
	/* SiP + OH =>  SiOP + H  */  
		
	co.SiP_OH_SiOP_H_1 = hmrate(6.3e-10,0,0)*co.hevmol[ipOH]*co.hevmol[ipSIP];		
	co.hydro_source[ipMH] += (float)co.SiP_OH_SiOP_H_1 ;
	
	/* CHP + H2O =>  HCOP +  H2  */   
		
	co.CHP_H2O_HCOP_H2_1 = hmrate(2.9e-9,0,0)*co.hevmol[ipH2O]*co.hevmol[ipCHP];			
	co.hydro_source[ipMH2g] += (float)co.CHP_H2O_HCOP_H2_1 ;
	
	/* CHP + OH =>  COP + H2  */   
		
	co.CHP_OH_COP_H2_1 = hmrate(7.5e-10,0,0)*co.hevmol[ipOH]*co.hevmol[ipCHP];		
	co.hydro_source[ipMH2g] += (float)co.CHP_OH_COP_H2_1 ;
	
	/* H + C =>  CH + nu  */
		
	co.H_C_CH_nu = hmrate(0.00000000000000001,0,0)*dense.xIonDense[ipHYDROGEN][0]*co.hevmol[ipATC];
	co.hydro_sink[ipMH] += (float)co.H_C_CH_nu ;
	
	/* H + CP =>  CHP + nu  */
		
	co.H_CP_CHP_nu = hmrate(1.7e-17,0,0)*co.hevmol[ipCP]*dense.xIonDense[ipHYDROGEN][0];	
	co.hydro_sink[ipMH] += (float)co.H_CP_CHP_nu ;
	
	/* H + OH =>  H2O + nu  */
		
	co.H_OH_H2O_nu = hmrate(5.26E-18,-5.22,90)*co.hevmol[ipOH]*dense.xIonDense[ipHYDROGEN][0];	
	co.hydro_sink[ipMH] += (float)co.H_OH_H2O_nu ;
	
	/* Hminus + CH =>  CH2 + e  */
		
	co.Hminus_CH_CH2_e = hmrate(1.0e-10,0,0)*co.hevmol[ipCH]*hmi.Hmolec[ipMHm];		
	co.hydro_sink[ipMHm] += (float)co.Hminus_CH_CH2_e ;
	
	/* Hminus + C =>  CH + e  */
		
	co.Hminus_C_CH_e = hmrate(1.0e-9,0,0)*co.hevmol[ipATC]*hmi.Hmolec[ipMHm];	
	co.hydro_sink[ipMHm] += (float)co.Hminus_C_CH_e ;
	
	/* Hminus + OH =>  H2O + e  */
		
	co.Hminus_OH_H2O_e = hmrate(1.0e-10,0,0)*co.hevmol[ipOH]*hmi.Hmolec[ipMHm];			
	co.hydro_sink[ipMHm] += (float)co.Hminus_OH_H2O_e ;
	
	/* Hminus + O =>  OH + e  */
		
	co.Hminus_O_OH_e = hmrate(1.0e-9,0,0)*co.hevmol[ipATO]*hmi.Hmolec[ipMHm];				
	co.hydro_sink[ipMHm] += (float)co.Hminus_O_OH_e ;
	
	/* H2 + C =>  CH2 + nu  */
		
	co.H2_C_CH2_nu = hmrate(1.0e-17,0,0)*co.hevmol[ipATC]*hmi.Hmolec[ipMH2g];		
	co.hydro_sink[ipMH2g] += (float)co.H2_C_CH2_nu ;
	
	/* H2 + CP =>  CH2P + nu  */
		
	co.H2_CP_CH2P_nu = hmrate(4.0e-16,-0.2,0)*co.hevmol[ipCP]*hmi.Hmolec[ipMH2g];			
	co.hydro_sink[ipMH2g] += (float)co.H2_CP_CH2P_nu ;
	
	/* H2 + SiP =>  SiH2P + nu  */
		
	co.H2_SiP_SiH2P_nu = hmrate(3.0e-18,0,0)*co.hevmol[ipSIP]*hmi.Hmolec[ipMH2g];	
	co.hydro_sink[ipMH2g] += (float)co.H2_SiP_SiH2P_nu ;
	
	/* H2 + O2 =>  OH + OH */
		
	co.H2_O2_OH_OH = hmrate(3.16e-10,0,21890)*co.hevmol[ipO2]*hmi.Hmolec[ipMH2g];
	co.hydro_sink[ipMH2g] += (float)co.H2_O2_OH_OH ;
	
	/* HeP + CH =>  CP + He + H  */
		
	co.HeP_CH_CP_He_H = hmrate(1.1e-9,0,0)*co.hevmol[ipCH]*dense.xIonDense[ipHELIUM][1];	
	co.hydro_source[ipMH] += (float)co.HeP_CH_CP_He_H ;
	
	/* HeP + CH2 =>  CHP + He + H  */
		
	co.HeP_CH2_CHP_He_H = hmrate(7.5e-10,0,0)*co.hevmol[ipCH2]*dense.xIonDense[ipHELIUM][1];	
	co.hydro_source[ipMH] += (float)co.HeP_CH2_CHP_He_H ;
	
	/* HeP + OH =>  OP + He + H  */
		
	co.HeP_OH_OP_He_H = hmrate(1.1e-9,0,0)*co.hevmol[ipOH]*dense.xIonDense[ipHELIUM][1];	
	co.hydro_source[ipMH] += (float)co.HeP_OH_OP_He_H ;
	
	/* HeP + H2O  OHP + He + H  */

	co.HeP_H2O_OHP_He_H = hmrate(2.86e-10,0,0)*co.hevmol[ipH2O]*dense.xIonDense[ipHELIUM][1];	
	co.hydro_source[ipMH] += (float)co.HeP_H2O_OHP_He_H ;
	
	/* HeP + SiH =>  SiP + He + H  */
	
	co.HeP_SiH_SiP_He_H = hmrate(1.8e-9,0,0)*co.hevmol[ipSIH]*dense.xIonDense[ipHELIUM][1];	
	co.hydro_source[ipMH] += (float)co.HeP_SiH_SiP_He_H ;
	
	/* HeP + H2O =>  OH + He + HP */
		
	co.HeP_H2O_OH_He_HP = hmrate(2.04e-10,0,0)*co.hevmol[ipH2O]*dense.xIonDense[ipHELIUM][1];	
	co.hydro_source[ipMHp] += (float)co.HeP_H2O_OH_He_HP ;
	
	/* HeP + CH2 =>  CP + He + H2  */
		
	co.HeP_CH2_CP_He_H2 = hmrate(7.5e-10,0,0)*co.hevmol[ipCH2]*dense.xIonDense[ipHELIUM][1];	
	co.hydro_source[ipMH2g] += (float)co.HeP_CH2_CP_He_H2 ;
	
	/* crnu + CH =>  C + H  */
		
	co.crnu_CH_C_H = co.hevmol[ipCH]*secondaries.csupra[ipHYDROGEN ][0] * 2. * 756;
	co.hydro_source[ipMH] += (float)co.crnu_CH_C_H ;
	
	/* crnu + CHP =>  CP + H  */
		
	co.crnu_CHP_CP_H = co.hevmol[ipCHP]*secondaries.csupra[ipHYDROGEN][0] * 2. * 183;	
	co.hydro_source[ipMH] += (float)co.crnu_CHP_CP_H ;
	
	/* crnu + H2O =>  OH + H  */
		
	co.crnu_H2O_OH_H = co.hevmol[ipH2O]*secondaries.csupra[ipHYDROGEN][0] * 2. * 979;	
	co.hydro_source[ipMH] += (float)co.crnu_H2O_OH_H ;
	
	/* crnu + OH =>  O + H  */
		
	co.crnu_OH_O_H = co.hevmol[ipOH]*secondaries.csupra[ipHYDROGEN][0] * 2. *522;
	co.hydro_source[ipMH] += (float)co.crnu_OH_O_H ;
	
	/* crnu + SiH =>  Si + H  */
		
	co.crnu_SiH_Si_H = co.hevmol[ipSIH]*secondaries.csupra[ipHYDROGEN][0] * 2. *500;	
	co.hydro_source[ipMH] += (float)co.crnu_SiH_Si_H ;
	
	/* nu + CH =>  C + H  */
		
	co.nu_CH_C_H = hmrate(8.6e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH];
	co.hydro_source[ipMH] += (float)co.nu_CH_C_H ;
	
	/* nu + CHP =>  CP + H  */
		
	co.nu_CHP_CP_H = hmrate(2.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCHP];
	co.hydro_source[ipMH] += (float)co.nu_CHP_CP_H ;
	
	/* nu + CH2 =>  CH + H  */
		
	co.nu_CH2_CH_H = hmrate(7.2e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH2];	
	co.hydro_source[ipMH] += (float)co.nu_CH2_CH_H ;
	
	/* nu + CH2P =>  CHP + H  */
		
	co.nu_CH2P_CHP_H = hmrate(1.7e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH2P];		
	co.hydro_source[ipMH] += (float)co.nu_CH2P_CHP_H ;
	
	/* nu + CH3P =>  CH2P + H  */
		
	co.nu_CH3P_CH2P_H = hmrate(1.0e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH3P];		
	co.hydro_source[ipMH] += (float)co.nu_CH3P_CH2P_H ;
	
	/* nu + CH3P =>  CHP + H2  */
		
	co.nu_CH3P_CHP_H2 = hmrate(1.0e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH3P];			
	co.hydro_source[ipMH2g] += (float)co.nu_CH3P_CHP_H2 ;
	
	/* nu + H2O =>  OH + H  */
		
	co.nu_H2O_OH_H = hmrate(5.9e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipH2O];		
	co.hydro_source[ipMH] += (float)co.nu_H2O_OH_H ;
	
	/* nu + OH =>  O + H  */
		
	co.nu_OH_O_H = hmrate(3.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipOH];			
	co.hydro_source[ipMH] += (float)co.nu_OH_O_H ;
	
	/* nu + OHP =>  O + HP */
		
	co.nu_OHP_O_HP = hmrate(1.0e-12,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipOHP];		
	co.hydro_source[ipMHp] += (float)co.nu_OHP_O_HP ;
	
	/* nu + SiH =>  Si + H  */
		
	co.nu_SiH_Si_H = hmrate(2.8e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipSIH];	
	co.hydro_source[ipMH] += (float)co.nu_SiH_Si_H ;
	
	/* e + CHP =>  C + H  */
		
	co.e_CHP_C_H = hmrate(1.5e-7,-0.42,0)*dense.eden*co.hevmol[ipCHP];	
	co.hydro_source[ipMH] += (float)co.e_CHP_C_H ;
	
	/* e + CH2P =>  CH + H  */
		
	co.e_CH2P_CH_H = hmrate(1.6e-7,-0.6,0)*dense.eden*co.hevmol[ipCH2P];	
	co.hydro_source[ipMH] += (float)co.e_CH2P_CH_H ;
	
	/* e + CH2P =>  C + H + H  */
		
	co.e_CH2P_C_H_H = hmrate(4.03e-7,-0.6,0)*dense.eden*co.hevmol[ipCH2P];	
	co.hydro_source[ipMH] += (float)co.e_CH2P_C_H_H ;
	
	/* e + CH2P =>  C + H2  */
		
	co.e_CH2P_C_H2 = hmrate(7.68e-8,-0.6,0)*dense.eden*co.hevmol[ipCH2P];	
	co.hydro_source[ipMH2g] += (float)co.e_CH2P_C_H2 ;
	
	/* e + CH3P =>  C + H2 + H  */
		
	co.e_CH3P_C_H2_H = hmrate(1.05e-7,-0.5,0)*dense.eden*co.hevmol[ipCH3P];		
	co.hydro_source[ipMH] += (float)co.e_CH3P_C_H2_H ;
	co.hydro_source[ipMH2g] += (float)co.e_CH3P_C_H2_H ;
	
	/* e + CH3P =>  CH2 + H  */
		
	co.e_CH3P_CH2_H	= hmrate(1.4e-7,-0.5,0)*dense.eden*co.hevmol[ipCH3P];
	co.hydro_source[ipMH] += (float)co.e_CH3P_CH2_H ;
	
	/* e + CH3P =>  CH + H + H  */
		
	co.e_CH3P_CH_H_H = hmrate(5.6e-8,-0.5,0)*dense.eden*co.hevmol[ipCH3P];	
	co.hydro_source[ipMH] += (float)co.e_CH3P_CH_H_H ;
	
	/* e + CH3P =>  CH + H2  */
		
	co.e_CH3P_CH_H2	= hmrate(4.9e-8,-0.5,0)*dense.eden*co.hevmol[ipCH3P];
	co.hydro_source[ipMH2g] += (float)co.e_CH3P_CH_H2 ;
	
	/* e + H2OP =>  OH + H  */
		
	co.e_H2OP_OH_H = hmrate(7.92e-8,-0.5,0)*dense.eden*co.hevmol[ipH2OP];
	co.hydro_source[ipMH] += (float)co.e_H2OP_OH_H ;
	
	/* e + H2OP =>  O + H + H  */
		
	co.e_H2OP_O_H_H = hmrate(2.45e-7,-0.5,0)*dense.eden*co.hevmol[ipH2OP];
	co.hydro_source[ipMH] += (float)co.e_H2OP_O_H_H ;
	
	/* e + H2OP =>  O + H2  */
		
	co.e_H2OP_O_H2 = hmrate(3.6e-8,-0.5,0)*dense.eden*co.hevmol[ipH2OP];
	co.hydro_source[ipMH2g] += (float)co.e_H2OP_O_H2 ;
	
	/* e + H3OP =>  H2O + H  */
		
	co.e_H3OP_H2O_H = hmrate(1.08e-7,-0.5,0)*dense.eden*co.hevmol[ipH3OP];
	co.hydro_source[ipMH] += (float)co.e_H3OP_H2O_H ;
	
	/* e + H3OP =>  OH + H + H  */
		
	co.e_H3OP_OH_H_H = hmrate(2.58e-7,-0.5,0)*dense.eden*co.hevmol[ipH3OP];	
	co.hydro_source[ipMH] += (float)co.e_H3OP_OH_H_H ;
	
	/* e + H3OP =>  OH + H2  */
		
	co.e_H3OP_OH_H2	= hmrate(6.45e-8,-0.5,0)*dense.eden*co.hevmol[ipH3OP];
	co.hydro_source[ipMH2g] += (float)co.e_H3OP_OH_H2 ;
	
	/* e + H3OP =>  O + H2 + H  */
		
	co.e_H3OP_O_H2_H = hmrate(5.59e-9,-0.5,0)*dense.eden*co.hevmol[ipH3OP];	
	co.hydro_source[ipMH] += (float)co.e_H3OP_O_H2_H ;
	
	/* e + HCOP =>  CO + H  */
		
	co.e_HCOP_CO_H = hmrate(1.1e-7,-1,0)*dense.eden*co.hevmol[ipHCOP];
	co.hydro_source[ipMH] += (float)co.e_HCOP_CO_H ;
	
	/* e + OHP =>  O + H  */
		
	co.e_OHP_O_H = hmrate(3.75e-8,-0.5,0)*dense.eden*co.hevmol[ipOHP];
	co.hydro_source[ipMH] += (float)co.e_OHP_O_H ;
	
	/* e + SiH2P =>  SiH + H  */
		
	co.e_SiH2P_SiH_H = hmrate(1.5e-7,-0.5,0)*dense.eden*co.hevmol[ipSIH2P];	
	co.hydro_source[ipMH] += (float)co.e_SiH2P_SiH_H ;
	
	/* e + SiH2P =>  Si + H + H  */
		
	co.e_SiH2P_Si_H_H = hmrate(2.0e-7,-0.5,0)*dense.eden*co.hevmol[ipSIH2P];	
	co.hydro_source[ipMH] += (float)co.e_SiH2P_Si_H_H ;
	
	/* e + SiH2P =>  Si + H2  */
		
	co.e_SiH2P_Si_H2 = hmrate(1.5e-7,-0.5,0)*dense.eden*co.hevmol[ipSIH2P];	
	co.hydro_source[ipMH2g] += (float)co.e_SiH2P_Si_H2 ;
	
	/* e + SiOHP =>  SiO + H  */
		
	co.e_SiOHP_SiO_H = hmrate(1.5e-7,-0.5,0)*dense.eden*co.hevmol[ipHOSIP];	
	co.hydro_source[ipMH] += (float)co.e_SiOHP_SiO_H ;
	
	/* H2 + CH =>  CH3 + nu  */
		
	co.H2_CH_CH3_nu	= hmrate(5.09E-18,-0.71,11.6)*co.hevmol[ipMH2g]*co.hevmol[ipCH];
	co.hydro_sink[ipMH2g] += (float)co.H2_CH_CH3_nu ;
	
	/* H2 + CH3P =>  CH5P + nu  */
		
	co.H2_CH3P_CH5P_nu = hmrate(1.3e-15,-1,0)*co.hevmol[ipCH3P]*hmi.Hmolec[ipMH2g];	
	co.hydro_sink[ipMH2g] += (float)co.H2_CH3P_CH5P_nu ;
	
	/* H2s + CH =>  CH3 + nu  */
		
	co.H2s_CH_CH3_nu = hmrate(5.09E-18,-0.71,0)*co.hevmol[ipCH]*hmi.Hmolec[ipMH2s]*hmi.lgLeiden_Keep_ipMH2s ;	
	co.hydro_sink[ipMH2s] += (float)co.H2s_CH_CH3_nu ;
	
	/* Hminus + CH2 =>  CH3 + e  */
		
	co.Hminus_CH2_CH3_e	= hmrate(1.0e-9,0,0)*co.hevmol[ipCH2]*hmi.Hmolec[ipMHm];	
	co.hydro_sink[ipMHm] += (float)co.Hminus_CH2_CH3_e ;
	
	/* Hminus + CH3 =>  CH4 + e  */
		
	co.Hminus_CH3_CH4_e	= hmrate(1.0e-9,0,0)*co.hevmol[ipCH3]*hmi.Hmolec[ipMHm];	
	co.hydro_sink[ipMHm] += (float)co.Hminus_CH3_CH4_e ;
	
	/* nu + CH3 =>  CH2 + H  */
		
	co.nu_CH3_CH2_H	= hmrate(2.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH3];
	co.hydro_source[ipMH] += (float)co.nu_CH3_CH2_H ;
	
	/* nu + CH3 =>  CH + H2  */
		
	co.nu_CH3_CH_H2	= hmrate(2.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH3];
	co.hydro_source[ipMH2g] += (float)co.nu_CH3_CH_H2 ;
	
	/* nu + CH4 =>  CH3 + H  */
		
	co.nu_CH4_CH3_H	= hmrate(2.2e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH4];
	co.hydro_source[ipMH] += (float)co.nu_CH4_CH3_H ;
	
	/* nu + CH4 =>  CH2 + H2  */
		
	co.nu_CH4_CH2_H2 = hmrate(9.8e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH4];	
	co.hydro_source[ipMH2g] += (float)co.nu_CH4_CH2_H2 ;
	
	/* nu + CH4 =>  CH + H2  */
		
	co.nu_CH4_CH_H2	= hmrate(2.2e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*co.hevmol[ipCH4];
	co.hydro_source[ipMH2g] += (float)co.nu_CH4_CH_H2 ;
	
	/* crnu + CH3 =>  CH2 + H  */
		
	co.crnu_CH3_CH2_H = co.hevmol[ipCH3]*secondaries.csupra[ipHYDROGEN][0] * 2.*500;	
	co.hydro_source[ipMH] += (float)co.crnu_CH3_CH2_H ;
	
	/* crnu + CH3 =>  CH + H2  */
		
	co.crnu_CH3_CH_H2 = co.hevmol[ipCH3]*secondaries.csupra[ipHYDROGEN][0] * 2.*500;	
	co.hydro_source[ipMH2g] += (float)co.crnu_CH3_CH_H2 ;
	
	/* crnu + CH4 =>  CH2 + H2  */
		
	co.crnu_CH4_CH2_H2 = co.hevmol[ipCH4]*secondaries.csupra[ipHYDROGEN][0] * 2.*2272;	
	co.hydro_source[ipMH2g] += (float)co.crnu_CH4_CH2_H2 ;
	
	/* e + CH5P =>  CH3 + H2  */
		
	co.e_CH5P_CH3_H2 = hmrate(5.5e-7,-0.3,0)*dense.eden*co.hevmol[ipCH5P];	
	co.hydro_source[ipMH2g] += (float)co.e_CH5P_CH3_H2 ;
	
	/* e + CH5P =>  CH4 + H  */
		
	co.e_CH5P_CH4_H	= hmrate(5.5e-7,-0.3,0)*dense.eden*co.hevmol[ipCH5P];
	co.hydro_source[ipMH] += (float)co.e_CH5P_CH4_H ;
	
	/* e + CH4P =>  CH3 + H  */
		
	co.e_CH4P_CH3_H	= hmrate(1.75e-7,-0.5,0)*dense.eden*co.hevmol[ipCH4P];
	co.hydro_source[ipMH] += (float)co.e_CH4P_CH3_H ;
	
	/* e + CH4P =>  CH2 + H + H  */
		
	co.e_CH4P_CH2_H_H = hmrate(1.75e-7,-0.5,0)*dense.eden*co.hevmol[ipCH4P];	
	co.hydro_source[ipMH] += (float)co.e_CH4P_CH2_H_H ;
	
	/******* END OF REACTIONS ********/
#	if 0
	if( iteration>1 )
	{
		r = rlist;
		i = 0;
		while (r->next != NULL)
		{
			++i;
			r = r->next;
			fprintf(ioQQQ,"DEBUG r\t%li\t%.3e\n", i, r->ratek);
		}
	}
#	endif

	/* Generate chemical error vector and Jacobian array from reaction list */
	for( i=0; i < N_H_MOLEC; i++ )
	{
		bvec[i] = 0.;
		for( j=0; j < N_H_MOLEC; j++ )
		{
			c[j][i] = 0.;
		}
	}
	/* Subtotal rates for H_0 and H^+ within ipMHo */
	for (i=0;i<2;i++) 
	{
		mole.source[ipHYDROGEN][i] = mole.sink[ipHYDROGEN][i] = 0.;
	}

	/* reinitialize linked list and move through it */
	/* set up Jacobian matrix */
	rtot = 0.;
	r = rlist;

	/*  Comments about this line of code, made by Nick Abel.  This section of the code
	    does the same thing as in co.c, where it decouples products of two densities
		that the code is trying to predict. Reactions of the form n(X)*n(Y) cannot be 
		treated properly in a linear solver otherwise.
		
		This section of the code starts off with a while
		statement that loops over all reactions in the network (hence the r->next variable).
		For each reaction, the first thing that is done is to state that for the reaction
		of interest, the rate coefficient k is equal to r->ratek, which is the stored value from
		when the reaction was originally set up above.  After this a loop over all reactants is 
		performed.  This is actually a double loop.  If the number of reactants is greater than 
		one, then this code will generate (via the i!=j if statement) two "decoupled" products of 
		the rate coefficient times the previous solution for the density.  This is explained in co.c 
		and in that code it is just the reaction rates that end in _1 or _2. 
		
		After this loop the	variable "rate" is multiplied by the the density that was not multiplied
		in the loop.  This sets up "rate" to be what co.c calls "bvec" for that reaction. The bvec's 
		are needed because in the process of "decoupling" a reaction, their is a leftover term of
		the form k*n(X)old*n(Y)old which has to go into the solution vector for the matrix equation Ax = b.
		
		The next two loops set up the bvec so that it goes with the proper species.  Also, in the
		case of H and H+, these rates are saved so that they can be fed into the main ionization 
		solver.  This allows Cloudy to account for formation and destruction processes
		for Hydrogen that are due to reaction with molecules.  

		The last for statement inside the while loop takes the rates that were calculated previously 
		and stores them in the appropriate part of the matrix A. */
		
	/* all reactions have been stored, now fill in the c[][] matrix */
	while (r->next != NULL)
	{
		r = r->next;
		/* >>chng 04 feb 05, this was option to cut chemistry short in testing
		if (r->index == rindex)
			break;*/
		ratek = r->ratek;
		/* There's an O(n) algorithm for this -- but it doesn't improve
		 * things unless nreactants is >= 4...!*/
		/* loop over all rate determinig species */
		for (i=0;i<r->nrates;i++)
		{
			rate_deriv[i] = ratek;
			for (j=0;j<r->nrates;j++)
			{
				/* Hmolec_old was previous abundance,
				 * rate_deriv[i] is deriv of rates coefficient by species r->rate_species[i] */
				if (i!=j)
				{
					rate_deriv[i] *= Hmolec_old[r->rate_species[j]];
				}
			}
		}

		/* this is total rate, rate_deriv times old population */
		rate = rate_deriv[0]*Hmolec_old[r->rate_species[0]];

		/* sum of all rates used for error scaling and wrning flags */
		rtot += rate;

		/* Get sink terms (i.e. rate/abundance) to include in ionization ladders */
		for (i=0;i < r->nreactants;i++)
		{
			int ok = 0;
			for (j=0;j < r->nrates && !ok;j++)
			{
				if (r->rate_species[j] == r->reactants[i]) 
				{
					
					sinkrate[i] = rate_deriv[j];
					ok = TRUE;
				}
			}
			if (!ok) 
			{
				/* Odd, the rate didn't depend on one of the species it used
				 * at all!  An alternative way of getting sinkrate is
				 *
				 * sinkrate[i] = rate/Hmolec_old[r->reactants[i]]; 
				 *
				 * but this uses the possibly zero Hmolec_old, and is prone to underflow of rate.
				 * */
				fprintf(ioQQQ,"A chemical rate in hmole was independent of the species it used\n");
				fprintf(ioQQQ,"This probably shouldn't happen (so you shouldn't see this message).\n");
				puts( "[Stop in hmole]" );
				cdEXIT(EXIT_FAILURE);
			}
		}


		/* if (nzone == 416)
			 fprintf(ioQQQ,"Adding reaction %d rate %g\n",r->index,rate); */

		/* loop over all reactions, find total consuption rate,
		 * also keep track of rates that use up H0 or H+ */
		for (i=0;i<r->nreactants;i++)
		{
			ratei = r->reactants[i];
			bvec[ratei] -= rate;
			/*if ((nzone == 421 || nzone == 422) && ratei == ipMHm)
					fprintf(ioQQQ,"snk %s %d %g\n",hmi.chLab[ratei],r->index,rate);*/
			/* mole.sink[ipHYDROGEN] is how chemical reaction network reacts with ionization
			 * network, so this keeps track of total rates */
			if (ratei == ipMH || ratei == ipMHp)
				mole.sink[ipHYDROGEN][ratei] += sinkrate[i];
		}

		/* loop over all reactions, find total production rate,
		 * also keep track of rates that produce H0 or H+ */
		for (i=0;i<r->nproducts;i++)
		{
			ratei = r->products[i];
			bvec[ratei] += rate;
			/*if ((nzone == 421 || nzone == 422) && ratei == ipMHm)
					fprintf(ioQQQ,"src %s %d %g\n",hmi.chLab[ratei],r->index,rate); */
			if (ratei == ipMH || ratei == ipMHp)
				mole.source[ipHYDROGEN][ratei] += rate;
		}

		/* fill Jacobian rate matrix */
		for (j=0;j<r->nrates;j++)
		{
			ratej = r->rate_species[j];
			rated = rate_deriv[j];
			for (i=0;i<r->nreactants;i++)
			{
				c[ratej][r->reactants[i]] -= rated;
			}
			for (i=0;i<r->nproducts;i++)
			{
				c[ratej][r->products[i]] += rated;
			}
		}
	}
	/* the c[][] matrix Jacobian array has now been filled with all reagents */

	/* save rate H2 is destroyed */
	hmi.H2_rate_destroy = -c[ipMH2g][ipMH2g];

#	if 0
	/* add in some more reactions from co.c that create or destroy a 
	   molecule predicted in hmole and belong in the b vector part of the
	   matrix equation Ax = b */

	bvec[ipMH2g] -= co.hydro_sink[ipMH2g];
	bvec[ipMH2g] += co.hydro_source[ipMH2g];
	
	bvec[ipMH2s] -= co.hydro_sink[ipMH2s];
	bvec[ipMH2s] += co.hydro_source[ipMH2s];
	
	bvec[ipMHm] -= co.hydro_sink[ipMHm];
	bvec[ipMHm] += co.hydro_source[ipMHm];
	
	bvec[ipMH] -= co.hydro_sink[ipMH];
	bvec[ipMH] += co.hydro_source[ipMH];
	
	bvec[ipMHp] -= co.hydro_sink[ipMHp];
	bvec[ipMHp] += co.hydro_source[ipMHp];
	
	mole.source[ipHYDROGEN][ipMH] +=  co.hydro_source[ipMH];
	mole.sink[ipHYDROGEN][ipMH]   +=  co.hydro_sink[ipMH];

	mole.source[ipHYDROGEN][ipMHp] +=  co.hydro_source[ipMHp];
	mole.sink[ipHYDROGEN][ipMHp]   +=  co.hydro_sink[ipMHp];
#	endif

	{
		/* following should be set true to print populations */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			if( DEBUG_LOC && (nzone > 570) ) 
			{
				printsol = 1;
				fprintf(ioQQQ,"Temperature %g\n",phycon.te);
				fprintf(ioQQQ," Net mol ion rate [%g %g] %g\n",mole.source[ipHYDROGEN][1],mole.sink[ipHYDROGEN][1],
								mole.source[ipHYDROGEN][1]-mole.sink[ipHYDROGEN][1]*Hmolec_old[ipMHp]);
			}
		}
	}

	/* save total H2P destruction rate for possible later printout:
	 * NB this must come last */
	desh2p = -c[ipMH2p][ipMH2p]; 

	/* Check that matrix and vector generated in above loops make sense */
	/*if (!defined(NDEBUG))	*/
	/* in std C NDEBUG is a macro set at compile time */
	/*TODO	0	had to comment following test out - NA change to hmole caused
	 * massive prints */
#	if 0
/*#	if !defined(NDEBUG)*/
#	ifndef NDEBUG
	{
		double total, mtotal;
		for (i=0;i<N_H_MOLEC;i++) 
		{
			total = 0.;
			for( j=0;j<N_H_MOLEC;j++) 
			{
				total += c[i][j]*nprot[j];
			}
			if( fabs(total) > 1e-5*fabs(c[i][i]*nprot[i])) 
			{
				fprintf(ioQQQ,"PROBLEM Subtotal1 %.2e\n",fabs(total)/fabs(c[i][i]*nprot[i]));
				fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*nprot[i]);
			}
			else if( fabs(total) > 1e-6*fabs(c[i][i]*nprot[i]) && phycon.te< 1e6 ) 
			{
				fprintf(ioQQQ,"NOTE Subtotal1 %.2e Te=%.4e\n",
					fabs(total)/fabs(c[i][i]*nprot[i]),phycon.te);
				fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*nprot[i]);
			}
		}
		total = mtotal = 0.;
		for (j=0;j<N_H_MOLEC;j++) 
		{ 
			total += bvec[j]*nprot[j]; 
			mtotal += fabs(bvec[j]*nprot[j]); 
		}
		if (fabs(total) > 1e-30 && fabs(total) > 1e-10*rtot) 
		{ 
			fprintf(ioQQQ,"PROBLEM Subtotal2 %.2e\n",fabs(total)/mtotal);
			fprintf(ioQQQ,"RHS Total %g cf %g\n",total,mtotal);
		} 
		else if (fabs(total) > 1e-7*mtotal)  
		{
			fprintf(ioQQQ,"WARNING zone %li Hmole RHS conservation error %.2e of %.2e\n",nzone,total,mtotal);
			fprintf(ioQQQ,"(may be due to high rate equilibrium reactions)\n");
		}
	}
#		endif
#	endif


#define MOLMIN  1
#define N_H_MAT (N_H_MOLEC-MOLMIN)
	/* Will collapse ipMH and ipMHp into single species, as don't include
	 * all ionizations and recombinations here */
	/* last test - do not include advection if we have overrun the radius scale 
	 * of previous iteration */
	if( iteration >= 2 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth ) 
	{
		/* Don't use conservation form in matrix solution */
		ipConserve = -1; 
		/* Add rate terms for dynamics to equilibrium, makes c[][] non-singular */
		for (i=0;i<N_H_MOLEC;i++) 
		{
			c[i][i] -= dynamics.Rate;
			bvec[i] -= (Hmolec_old[i]*dynamics.Rate-dynamics.Molec[i]);
		}

		/* Dynamics implies conservation of advected material */
		protons = 0.;
		for (i=0; i<N_H_MOLEC;i++)
		{
			protons += nprot[i]*dynamics.Molec[i]/dynamics.Rate;
		}
		/* bring H+ and H0 together since their ratio is set in H atom solver,
		 * we determine sum of two here */
		for (i=0;i<N_H_MOLEC;i++) 
		{
			c[ipMHp][i] = (Hmolec_old[ipMH]*c[ipMH][i]+Hmolec_old[ipMHp]*c[ipMHp][i])/sum_H0_Hp;
			c[ipMH][i] = 0.;
		}
		for (i=1;i<N_H_MOLEC;i++) 
		{
			c[i][ipMHp] += c[i][ipMH];
			c[i][ipMH] = 0.;
		}
		bvec[ipMHp] += bvec[ipMH];
		bvec[ipMH] = 0.;
		Hmolec_old[ipMHp] += Hmolec_old[ipMH];
		Hmolec_old[ipMH] = 0.;
	}
	else
	{
		/* usual branch, no advection */
		/* bring H+ and H0 together since their ratio is set in H atom solver,
		 * we determine sum of two here */
		for (i=0;i<N_H_MOLEC;i++) 
		{
			/* >>chng 04 feb 04, sum_H0_Hp goes to zero when no ionization,
			 * add test on SMALLFLOAT */
			if( sum_H0_Hp > SMALLFLOAT )
				c[ipMHp][i] = (Hmolec_old[ipMH]*c[ipMH][i]+Hmolec_old[ipMHp]*c[ipMHp][i])/sum_H0_Hp;
			c[ipMH][i] = 0.;
		}
		Hmolec_old[ipMHp] += Hmolec_old[ipMH];
		bvec[ipMH] = Hmolec_old[ipMH] = 0.;
		ipConserve = ipMHp;
		/* For Newton-Raphson method, want the change in populations to be zero,
		 * so the conserved component must also be zero */
		bvec[ipConserve] = 0.;  
		protons = 0.;
		for (i=MOLMIN;i<N_H_MOLEC;i++) 
		{
			c[i][ipConserve] = nprot[i];
			protons += nprot[i]*Hmolec_old[i];
		}
	}

	{
		/* following should be set true to print populations */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			/* these are the raw results */
			fprintf( ioQQQ, " HMOLE h2 %.2e h2* %.2e\n" , Hmolec_old[ipMH2g] ,Hmolec_old[ipMH2s] );
		}
	}
	
	/*------------------------------------------------------------------ */
	if(printsol || (trace.lgTrace && trace.lgTr_H2_Mole ))
	{
		/* print the full matrix */
		fprintf( ioQQQ, "                ");
		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
			fprintf( ioQQQ, "      %s", hmi.chLab[i] );
		}
		fprintf( ioQQQ, " \n" );

		/*

		[0][0]  [0][1]  [0][2]  [0][3]  [0][4]  [0][5]
		[1][0]  [1][1]  [1][2]  [1][3]  [1][4]  [1][5]
		[2][0]  [2][1]  [2][2]  [2][3]  [2][4]  [2][5]
		[3][0]  [3][1]  [3][2]  [3][3]  [3][4]  [3][5]
		[4][0]  [4][1]  [4][2]  [4][3]  [4][4]  [4][5]
		[5][0]  [5][1]  [5][2]  [5][3]  [5][4]  [5][5]

		[ipMHo][ipMHo]  [ipMHo][ipMHm]  [ipMHo][ipMH2g]  [ipMHo][ipMH2p]  [ipMHo][ipMH3p]  [ipMHo][ipMH2s]
		[ipMHm][ipMHo] [ipMHm][ipMHm] [ipMHm][ipMH2g] [ipMHm][ipMH2p] [ipMHm][ipMH3p] [ipMHm][ipMH2s]
		[ipMH2g][ipMHo]  [ipMH2g][ipMHm]  [ipMH2g][ipMH2g]  [ipMH2g][ipMH2p]  [ipMH2g][ipMH3p]  [ipMH2g][ipMH2s]
		[ipMH2p][ipMHo] [ipMH2p][ipMHm] [ipMH2p][ipMH2g] [ipMH2p][ipMH2p] [ipMH2p][ipMH3p] [ipMH2p][ipMH2s]
		[ipMH3p][ipMHo] [ipMH3p][ipMHm] [ipMH3p][ipMH2g] [ipMH3p][ipMH2p] [ipMH3p][ipMH3p] [ipMH3p][ipMH2s]
		[ipMH2s][ipMHo] [ipMH2s][ipMHm] [ipMH2s][ipMH2g] [ipMH2s][ipMH2p] [ipMH2s][ipMH3p]  [ipMH2s][ipMH2s]

		*/

		fprintf(ioQQQ,"       MOLE old abundances\t%.2f",fnzone);
		for( i=0; i<N_H_MOLEC; i++ )
			fprintf(ioQQQ,"\t%.2e", Hmolec_old[i] );
		fprintf(ioQQQ,"\n" );

		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
			fprintf( ioQQQ, "       MOLE%2ld %s", i-MOLMIN ,hmi.chLab[i] );
			for( j=MOLMIN; j < N_H_MOLEC; j++ )
			{
				fprintf( ioQQQ, "%10.2e", c[j][i] );
			}
			fprintf( ioQQQ, "%10.2e", bvec[i] );
			fprintf( ioQQQ, "\n" );
		}
	}

	/*fprintf(ioQQQ,"DEBUG %.2e %.2e %.2e %.2e\n", 
		c[ipMH][ipMH2g] , c[ipMH][ipMH2s],
		c[ipMH2g][ipMH] , c[ipMH2s][ipMH]);
	fprintf(ioQQQ,"DEBUG %.2e %.2e %.2e %.2e\n\n", 
		c[ipMHp][ipMH2g] , c[ipMHp][ipMH2s],
		c[ipMH2g][ipMHp] , c[ipMH2s][ipMHp]);*/
	/* establish local timescale for H2 molecule destruction */
	if( -c[ipMH2g][ipMH2g] > SMALLFLOAT )
	{
		/* units are now seconds */
		timesc.time_H2_Dest_here = -1./c[ipMH2g][ipMH2g];
	}
	else
	{
		timesc.time_H2_Dest_here = 0.;
	}

	/* local timescale for H2 formation 
	 * both grains and radiative attachment */
	timesc.time_H2_Form_here = gv.rate_h2_form_grains_used_total * 
		/* this corrects for fact that we the timescale for H2 to form from an atomic gas.
		 * The rate becomes very small when gas is fully molecular, and ratio of total hydrogen
		 * to atomic hydrogen corrections for this. */
		dense.gas_phase[ipHYDROGEN]/SDIV(dense.xIonDense[ipHYDROGEN][0]) + 
		hmi.hminus_rad_attach;
	/* timescale is inverse of this rate */
	if( timesc.time_H2_Form_here > SMALLFLOAT )
	{
		/* units are now seconds */
		timesc.time_H2_Form_here = 1./timesc.time_H2_Form_here;
	}
	else
	{
		timesc.time_H2_Form_here = 0.;
	}

#	ifdef MAT
#		undef MAT
#	endif
#	define MAT(a,I_,J_)	(*((a)+(I_)*(N_H_MAT)+(J_)))

	/* copy contents over to 1D array */
	for( j=0; j < N_H_MAT; j++ )
	{
		for( i=0; i < N_H_MAT; i++ )
		{
			MAT(amat,i,j) = c[i+MOLMIN][j+MOLMIN];
		}
	}
	
	if (printsol)
	{
		double total=0;
		fprintf(ioQQQ,"Zone %.2f input:\n",fnzone);
		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
				fprintf(ioQQQ,"%15.7g\t",Hmolec_old[i]);
				total += nprot[i]*Hmolec_old[i];
		}
		fprintf(ioQQQ,"sum = %15.7g\n",total);
	}

	merror = 0;

	/* now invert the matrix */
	getrf_wrapper(N_H_MAT,N_H_MAT,(double*)amat,N_H_MAT,ipiv,&merror);
	/*DGETRF(N_H_MAT,N_H_MAT,(double*)amat,N_H_MAT,ipiv, &merror);*/

	getrs_wrapper('N',N_H_MAT,1,(double*)amat,N_H_MAT,ipiv,bvec+MOLMIN,N_H_MAT,&merror);
	/*DGETRS('N',N_H_MAT,1,(double*)amat,N_H_MAT,ipiv,bvec+MOLMIN,N_H_MAT,&merror);*/

	if( merror != 0 )
	{
		fprintf( ioQQQ, " hmole_step: dgetrs finds singular or ill-conditioned matrix\n" );
		puts( "[Stop in hmole]" );
		cdEXIT(EXIT_FAILURE);
	}

	if (printsol)
	{
		double total=0;
		fprintf(ioQQQ,"solution:\n");
		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
			fprintf(ioQQQ,"%15.7g\t",bvec[i]);
			total += nprot[i]*bvec[i];
		}
		fprintf(ioQQQ,"sum = %15.7g\n",total);
#if 0
		fprintf(ioQQQ,"Verify:\n");
		for (i=MOLMIN; i < N_H_MOLEC; i++ ) 
		{
			total = 0;
			for (j=MOLMIN; j<N_H_MOLEC; j++ ) 
			{
				total += c[j][i]*bvec[j];
				fprintf(ioQQQ,"%14.7g ",c[j][i]*bvec[j]);
			}
			fprintf(ioQQQ,"\t");
		  fprintf(ioQQQ,"%15.7g\n",total);
		}
		fprintf(ioQQQ,"\n");
#endif
	}

#	if 0
	/* this was fixup to limit changes in solution, but in the end was
	 * not used and is candidate for removal */
	/* For N-R iteration, output from matrix solution is now _change_ in solution */
	if(0 && !( iteration >= 2 && dynamics.lgAdvection )) 
	{
		/* Limit the size of the step to up by 100%, down by 50% */
		double nrfac = 1.;
		static int zone = -1;

		/* Be careful on first call */
		if (nzone == 0 && zone != nzone) 
		{
			nrfac = 1e-3;
			zone = nzone;
		}

		iworst = -1;
		for (i=MOLMIN; i< N_H_MOLEC; i++) 
		{
			if (bvec[i] > 0.) 
			{
				/* Only test `significant' species -- let others go -ve to be caught below */
				if (nrfac*bvec[i] > 0.5*Hmolec_old[i] && Hmolec_old[i] > 1e-6*protons) 
				{
					nrfac = 0.5*Hmolec_old[i]/bvec[i];
					iworst = i;
				}
			} 
			else 
			{
				if (-nrfac*bvec[i] > Hmolec_old[i] && Hmolec_old[i] > 1e-6*protons) 
				{
					nrfac = -Hmolec_old[i]/bvec[i];
					iworst = i;
				}
			}
		}
		/* fprintf(ioQQQ,"Step factor %g species %d\n",nrfac,iworst); */
		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
			bvec[i] *= nrfac;
		}
	}
#	endif

	*error = 0;
	/* H0 and H+ rates have been collapsed into one, since that solution
	 * comes from H atom solver.  */
	for( i=MOLMIN; i < N_H_MOLEC; i++ )
	{
		/* Smooth the error mode tailoff */
		etmp = bvec[i]/(ABSLIM+Hmolec_old[i]);

		if (printsol)
			fprintf(ioQQQ,"%15.7g\t",etmp);

		/* square of change in abundance of this species, in this iteration */
		*error += etmp*etmp;
		/* change bvec from being the difference into being the new value
		 * bvec is new soln */
		bvec[i] = Hmolec_old[i]-bvec[i];
	}
	*error = sqrt(*error)/N_H_MAT;

	if (printsol)
	{
		double total=0;
		fprintf(ioQQQ,"err = %15.7g\n",*error);
		/* fprintf(ioQQQ,"derived:\n"); */
		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
				fprintf(ioQQQ,"%15.7g\t",bvec[i]);
				total += nprot[i]*bvec[i];
		}
		fprintf(ioQQQ,"sum = %15.7g\n",total);
	}

	sum = 0.;
	/* check for negative populations */
	lgNegPop = FALSE;
	fracneg = 0.;
	fracnegfac = 0.;
	iworst = -1;
	for( i=MOLMIN; i < N_H_MOLEC; i++ )
	{
		if( bvec[i] < 0. ) {
			lgNegPop = TRUE;
		}
		/* largest relative change in solution for neg soln */
		fracnegtmp = -bvec[i]/SDIV(Hmolec_old[i]);
		/* this can only occur for negative solutions since fracneg starts
		 * as zero */
		if (fracnegtmp > fracneg) 
		{
			fracneg = fracnegtmp;
			fracnegfac = (0.5*Hmolec_old[i]-bvec[i])/(Hmolec_old[i]-bvec[i]);
			iworst = i;
		}
		sum += nprot[i] * bvec[i] ;
	}
	conserve = (protons - sum)/SDIV(protons);
	if(fabs(conserve) > 1e-8 )
		fprintf(ioQQQ,"PROBLEM hmoleee protons %.4e, sum %.4e (h-s)/h %.3e \n", protons , sum , 
			conserve );

#	if 0
	/* NDEBUG is set by the compiler to indicate that a debugging mode
	 * has not been specified.  */
#	ifndef NDEBUG
	/*if (NDEBUG)*/
	{
		fprintf( ioQQQ, "Zone %li\n",nzone);
		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
			fprintf(ioQQQ," %s %.2e", hmi.chLab[i], Hmolec_old[i]);
		}
		fprintf( ioQQQ, " =>\n" );
		for( i=MOLMIN; i < N_H_MOLEC; i++ )
		{
			fprintf(ioQQQ," %s %.2e", hmi.chLab[i], bvec[i]);
		}
		fprintf( ioQQQ, "\n" );
	}
#	endif
#	endif

	if(lgNegPop)
	{
#		ifndef NDEBUG
		/* very common to obtain negative soln on very first try - 
		 * don't print in this case */
		if(*nFixup )
		{
			fprintf( ioQQQ, " PROBLEM  hmole_step finds negative H molecule, in zone %.2f.\n",fnzone );
			fprintf( ioQQQ, " Worst is species %d -ve by fraction %g.\n",iworst,fracneg );
			fprintf( ioQQQ, " The populations follow:\n");
			for( i=MOLMIN; i < N_H_MOLEC; i++ )
			{
				fprintf(ioQQQ," %s %.2e", hmi.chLab[i], bvec[i]);
			}
			fprintf( ioQQQ, "\n" );
		}
#		endif

		/* Fix negative abundance -- assume the new solution is better in some ways */
		{
			double total=0., ntotal=0., ratio;
			enum {FIXUPTYPE = 1};

			++*nFixup;

			if (FIXUPTYPE == 1) {
				for( i=MOLMIN; i < N_H_MOLEC; i++ )
				{
					total += nprot[i]*bvec[i];
					if (bvec[i] < 0) 
					{
						ntotal += nprot[i]*bvec[i];
						bvec[i] = 0.;
					}
				}
				ratio = total/(total-ntotal);
				for( i=MOLMIN; i < N_H_MOLEC; i++ )
				{
					bvec[i] *= ratio;
				}
			}
			else if (FIXUPTYPE == 2) 
			{
				for( i=MOLMIN; i < N_H_MOLEC; i++ )
				{
					bvec[i] = fracnegfac*Hmolec_old[i]+(1-fracnegfac)*bvec[i];
				}
			}
				
#			ifndef NDEBUG
			/*if (NDEBUG)*/
			/* very common to obtain negative soln on very first try - 
			 * don't print in this case */
			if( *nFixup>1 )
			{
				fprintf(ioQQQ," FIXUP taken %i time%s.\n\n", *nFixup, (*nFixup == 1)?"":"s");
				fprintf( ioQQQ, " Initially:\n");
				for( i=MOLMIN; i < N_H_MOLEC; i++ )
				{
					fprintf(ioQQQ," %s %.2e", hmi.chLab[i], Hmolec_old[i]);
				}
				fprintf( ioQQQ, "\n" );
				fprintf( ioQQQ, " Changed to:\n");
				for( i=MOLMIN; i < N_H_MOLEC; i++ )
				{
					fprintf(ioQQQ," %s %.2e", hmi.chLab[i], bvec[i]);
				}
				fprintf( ioQQQ, "\n" );
			}
#			endif
		}
	}

	/* put derived abundances back into appropriate molecular species,
	 * bvec[ipMHp] is still sum of H0 and H+ from chemistry,
	 * now split up using ratio found in H atom solver */
	h1fnd = bvec[ipMHp];
	/* >>chng 04 feb 04, add SMALLFLOAT to protect against fully molecular limit */
	h1rat = h1fnd/SDIV(sum_H0_Hp);
	/* put back into proper places in soln vector so following loops work
	 * as expected */
	bvec[ipMH] = dense.xIonDense[ipHYDROGEN][0] * h1rat;
	bvec[ipMHp] = dense.xIonDense[ipHYDROGEN][1] * h1rat;
	/* ASSERT (fabs(bvec[ipMH]+bvec[ipMHp]-h1fnd) < 1e-12 * h1fnd); */

	if (fabs(bvec[ipMH]+bvec[ipMHp]-h1fnd) >= 1e-12 * h1fnd) 
	{
		fprintf(ioQQQ,"PROBLEM h1fnd %g %g %g %g %g\n",bvec[ipMH],bvec[ipMHp],h1fnd,h1rat,bvec[ipMH]+bvec[ipMHp]-h1fnd);
	}

	for (mol=0;mol<N_H_MOLEC;mol++) 
		hmi.Hmolec[mol] = (float) bvec[mol];

	dense.xIonDense[ipHYDROGEN][0] = (float) bvec[ipMH];
	dense.xIonDense[ipHYDROGEN][1] = (float) bvec[ipMHp];

	/* total H2 - all forms */
	hmi.H2_total = hmi.Hmolec[ipMH2s] + hmi.Hmolec[ipMH2g];
	/* first guess at ortho and para densities */
	h2.ortho_density = 0.75*hmi.H2_total;
	h2.para_density = 0.25*hmi.H2_total;

	/* NB the first index must be kept parallel with nelem and ionstag in
	 * H2Lines EmLine struc,
	 * since that struc expects to find the abundances here */
	/* >>chng 04 feb 19, had been ipMH2g, chng to total */
	dense.xIonDense[LIMELM+2][0] = hmi.H2_total;

	/* identify dominant H2 formation process */
	{
		/* following should be set true to identify H2 formation and destruction processes */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone>400) )
		{
			double createsum ,create_from_Hn2 , create_3body_Ho, create_h2p, 
				create_h3p, create_h3pe, create_grains, create_hminus;
			double destroysum, destroy_hm ,destroy_hm2 ,destroy_soloman ,destroy_2h ,destroy_hp,
				destroy_h,destroy_hp2,destroy_h3p;

			/* H(n=2) + H(n=1) -> H2 */
			create_from_Hn2 = radasc*dense.xIonDense[ipHYDROGEN][0];
			/* 3H => H2 + H */
			create_3body_Ho = bh2dis*dense.xIonDense[ipHYDROGEN][0];
			/* H2+ + H => H2 + H+ */
			create_h2p = hmi.bh2h2p*dense.xIonDense[ipHYDROGEN][0]*Hmolec_old[ipMH2p];
			/* H + H3+ => H2 + H2+ */
			create_h3p = h3ph2p*dense.xIonDense[ipHYDROGEN][0]*hmi.Hmolec[ipMH3p];
			/* e + H3+ => H2 + H */
			create_h3pe = eh3_h2h*dense.eden * hmi.Hmolec[ipMH3p];
			/* from grains */
			create_grains = gv.rate_h2_form_grains_used_total;
			/* from H- */
			create_hminus = Hmolec_old[ipMH]*hmi.assoc_detach*hmi.Hmolec[ipMHm];

			createsum = create_from_Hn2 + create_3body_Ho + create_h2p +
				create_h3p + create_h3pe + create_grains + create_hminus;

			fprintf(ioQQQ,"H2 create zone\t%.2f \tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				fnzone,
				createsum,
				create_hminus   / createsum,
				create_from_Hn2 / createsum,
				create_3body_Ho / createsum,
				create_h2p      / createsum,
				create_h3p      / createsum,
				create_h3pe     / createsum,
				create_grains   / createsum	);

			/* all h2 molecular hydrogen destruction processes */
			/* >>chng 04 jan 28, had wrong Boltzmann factor,
			 * caught by gargi Shaw */
			destroy_hm = assoc_detach_backwards_grnd+assoc_detach_backwards_exct;
			destroy_hm2 = eh2hhm;
			destroy_soloman = hmi.H2_Solomon_dissoc_rate_used;
			destroy_2h = eh2hh;
			destroy_hp = h2hph3p*dense.xIonDense[ipHYDROGEN][1];
			destroy_h = hmi.rh2dis*dense.xIonDense[ipHYDROGEN][0];
			destroy_hp2 = rh2h2p*dense.xIonDense[ipHYDROGEN][1];
			destroy_h3p = h3petc * hmi.Hmolec[ipMH3p];
			destroysum = destroy_hm + destroy_hm2 + destroy_soloman + destroy_2h + 
				destroy_hp+ destroy_h+ destroy_hp2+ destroy_h3p;

			fprintf(ioQQQ,"H2 destroy\t%.3f \t%.2e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				fnzone,
				destroysum,
				destroy_hm / destroysum ,
				destroy_hm2 / destroysum ,
				destroy_soloman / destroysum ,
				destroy_2h / destroysum ,
				destroy_hp / destroysum ,
				destroy_h / destroysum ,
				destroy_hp2 / destroysum ,
				destroy_h3p / destroysum );

		}
	}

	{
		/* following should be set true to identify H- formation and destruction processes */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone>140)/**/ )
		{
			double create_from_Ho,create_3body_Ho,create_batach,create_eh2hhm,destroy_photo,
				destroy_coll_heavies,destroy_coll_electrons,destroy_Hattach,destroy_fhneut,
				destsum , createsum;

			create_from_Ho = (hmi.hminus_rad_attach + hmi.HMinus_induc_rec_rate);
			create_3body_Ho = c3bod;
			/* total formation is sum of g and s attachment */
			create_batach = assoc_detach_backwards_grnd + assoc_detach_backwards_exct;
			create_eh2hhm = eh2hhm;
			destroy_photo = hmi.HMinus_photo_rate;
			destroy_coll_heavies = hmi.hmin_ct_firstions*sum_first_ions;
			destroy_coll_electrons = cionhm;
			destroy_Hattach = Hmolec_old[ipMH]*hmi.assoc_detach;
			destroy_fhneut = fhneut;

			destsum = destroy_photo + destroy_coll_heavies + destroy_coll_electrons + 
				destroy_Hattach + destroy_fhneut;
			fprintf(ioQQQ,"H- destroy zone\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", 
				fnzone,
				phycon.te,
				destsum,
				destroy_photo/destsum , 
				destroy_coll_heavies/destsum,
				destroy_coll_electrons/destsum, 
				destroy_Hattach/destsum,
				destroy_fhneut/destsum );

			createsum = create_from_Ho+create_3body_Ho+create_batach+create_eh2hhm;
			fprintf(ioQQQ,"H- create\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				fnzone,
				phycon.te,
				createsum,
				dense.eden,
				create_from_Ho/createsum,
				create_3body_Ho/createsum,
				create_batach/createsum,
				create_eh2hhm/createsum);
		}
	}

	/* rate H-alpha is created by H- ct */
	hmi.HalphaHmin = (float)(fhneut*hmi.Hmolec[ipMHm]);

	/* heating due to H2 dissociation */
	if( hmi.lgNoH2Mole )
	{
		hmi.HeatH2Dish_TH85 = 0.;
		hmi.HeatH2Dexc_TH85 = 0.;
	}
	else
	{
		/* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
		/*hmi.HeatH2Dish_TH85 = (float)(1.36e-23*hmi.Hmolec[ipMH2g]*h2esc*hmi.UV_Cont_rel2_Habing_TH85_depth);*/
		/* >>chng 04 feb 07, more general to express heating in terms of the assumed
		 * photo rates - the 0.25 was obtained by inverting A8 & A9 of TH85 to find that
		 * there are 0.25 eV per dissociative pumping, ie, 10% of total 
		 * this includes both H2g and H2s - TH85 say just grnd but they include
		 * process for both H2 and H2s - as we did above - both must be in
		 * heating term */
		hmi.HeatH2Dish_TH85 = (float)(0.25*EN1EV * hmi.H2_Solomon_dissoc_rate_used*
			hmi.H2_total);

		/* >>chng 00 nov 25, explictly break out this heat agent */
		/* 2.6 eV of heat per deexcitation, consider difference
		 * between deexcitation (heat) and excitation (cooling) */
		/* >>chng 04 jan 27, code moved here and simplified */
		deexc_htwo = (1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6.;
		deexc_hneut =  (1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.;

		H2star_deexcit = hmi.H2_total*deexc_htwo + hmi.Hmolec[ipMH] * deexc_hneut;
		H2star_excit = Boltz_fac_H2_H2star * H2star_deexcit;
		hmi.HeatH2Dexc_TH85 = (float)( (hmi.Hmolec[ipMH2s]*H2star_deexcit - hmi.Hmolec[ipMH2g]*H2star_excit) * 4.17e-12);
		/* leiden hacks try to turn off H2*, which is all unphysical.  do not includee heating
		 * due to H2 deexcitation since H2* is bogus */
		hmi.HeatH2Dexc_TH85 *= hmi.lgLeiden_Keep_ipMH2s;
		/*fprintf(ioQQQ,
			"DEBUG hmole H2 deex heating:\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			fnzone,
			hmi.HeatH2Dexc_TH85,
			thermal.htot,
			hmi.Hmolec[ipMH2s],
			H2star_deexcit, 
			hmi.Hmolec[ipMH2g],
			H2star_excit,
			Hmolec_old[ipMH2g],
			Hmolec_old[ipMH],
			phycon.te);*/
	}
	/* this is deriv wrt temperature, only if counted as a cooling source */
	hmi.deriv_HeatH2Dexc_TH85 = (float)(MAX2(0.,-hmi.HeatH2Dexc_TH85)* ( 30172. * thermal.tsq1 - thermal.halfte ) );

	{
		/* following should be set true to print populations */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			/* these are the raw results */
			fprintf( ioQQQ, " HMOLE raw; hi\t%.2e" , dense.xIonDense[ipHYDROGEN][0]);
			for( i=0; i < N_H_MOLEC; i++ )
			{
				fprintf( ioQQQ, "\t%s\t%.2e", hmi.chLab[i], bvec[i] );
			}
			fprintf( ioQQQ, " \n" );
		}
	}

	if( trace.lgTrace && trace.lgTr_H2_Mole )
	{
		/* these are the raw results */
		fprintf( ioQQQ, " raw; " );
		for( i=0; i < N_H_MOLEC; i++ )
		{
			fprintf( ioQQQ, " %s:%.2e", hmi.chLab[i], bvec[i] );
		}
		fprintf( ioQQQ, " \n" );

	}

	/* option to print rate H2 forms */
	/* trace.lgTr_H2_Mole is trace molecules option,
	 * punch htwo */
	if( (trace.lgTrace && trace.lgTr_H2_Mole) || punch.lgPunH2 )
	{
		if( punch.lgPunH2 )
		{
			ioFile = punch.ipPunH2;
		}
		else
		{
			ioFile = ioQQQ;
		}

		rate = gv.rate_h2_form_grains_used_total + 
			Hmolec_old[ipMH]*hmi.assoc_detach*hmi.Hmolec[ipMHm] + 
			bh2dis*dense.xIonDense[ipHYDROGEN][0] + 
			hmi.bh2h2p*dense.xIonDense[ipHYDROGEN][0]*hmi.Hmolec[ipMH2p] + 
			radasc*dense.xIonDense[ipHYDROGEN][0] + 
			h3ph2p*hmi.Hmolec[ipMH3p] + 
			h3petc*hmi.Hmolec[ipMH3p];

		if( rate > SMALLFLOAT )
		{
			fprintf( ioFile, 
				" Create H2, rate=%10.2e grain;%5.3f hmin;%5.3f bhedis;%5.3f h2+;%5.3f radasc:%5.3f h3ph2p:%5.3f h3petc:%5.3f\n", 
			  rate, 
			  gv.rate_h2_form_grains_used_total/rate, 
			  Hmolec_old[ipMH]*hmi.assoc_detach*hmi.Hmolec[ipMHm]/rate, 
			  bh2dis/rate, 
			  hmi.bh2h2p*dense.xIonDense[ipHYDROGEN][0]*hmi.Hmolec[ipMH2p]/rate, 
			  radasc*dense.xIonDense[ipHYDROGEN][0]/rate, 
			  h3ph2p*hmi.Hmolec[ipMH3p]/rate, 
			  h3petc*hmi.Hmolec[ipMH3p]/rate );
		}
		else
		{
			fprintf( ioFile, " Create H2, rate=0\n" );
		}
	}

	/* this is H2+ */
	if( trace.lgTrace && trace.lgTr_H2_Mole )
	{
		rate = rh2h2p*hmi.Hmolec[ipMH2g]*dense.xIonDense[ipHYDROGEN][1] + b2pcin*dense.xIonDense[ipHYDROGEN][1]*dense.eden*dense.xIonDense[ipHYDROGEN][0] + 
		  h3ph2p*hmi.Hmolec[ipMH3p] + h3petc*hmi.Hmolec[ipMH3p];
		if( rate > 1e-25 )
		{
			fprintf( ioQQQ, " Create H2+, rate=%10.2e rh2h2p;%5.3f b2pcin;%5.3f h3ph2p;%5.3f h3petc+;%5.3f\n", 
			  rate, rh2h2p*dense.xIonDense[ipHYDROGEN][1]*hmi.Hmolec[ipMH2g]/rate, 
				b2pcin*dense.xIonDense[ipHYDROGEN][1]*dense.xIonDense[ipHYDROGEN][1]*
			  dense.eden/rate, h3ph2p*hmi.Hmolec[ipMH3p]/rate, h3petc*hmi.Hmolec[ipMH3p]/
			  rate );
		}
		else
		{
			fprintf( ioQQQ, " Create H2+, rate=0\n" );
		}
	}

	if( hmi.Hmolec[ipMHm] > 0. && hmi.rel_pop_LTE_Hmin > 0. )
	{
		hmi.hmidep = (double)hmi.Hmolec[ipMHm]/ SDIV( 
			((double)dense.xIonDense[ipHYDROGEN][0]*dense.eden*hmi.rel_pop_LTE_Hmin));
	}
	else
	{
		hmi.hmidep = 1.;
	}

	/* this will be net volume heating rate, photo heat - induc cool */
	hmi.hmihet = hmi.HMinus_photo_heat*hmi.Hmolec[ipMHm] - hmi.HMinus_induc_rec_cooling;
	hmi.h2plus_heat = h2phet*hmi.Hmolec[ipMH2p];

	/* departure coef for H2 defined rel to n(1s)**2
	 * (see Littes and Mihalas Solar Phys 93, 23) */
	plte = (double)dense.xIonDense[ipHYDROGEN][0] * hmi.rel_pop_LTE_H2 * (double)dense.xIonDense[ipHYDROGEN][0];
	if( plte > 0. )
	{
		hmi.h2dep = hmi.Hmolec[ipMH2g]/plte;
	}
	else
	{
		hmi.h2dep = 1.;
	}

	/* departure coef of H2+ defined rel to n(1s) n(p)
	 * sec den was HI before 85.34 */
	plte = (double)dense.xIonDense[ipHYDROGEN][0]*hmi.rel_pop_LTE_H2p*(double)dense.xIonDense[ipHYDROGEN][1];
	if( plte > 0. )
	{
		hmi.h2pdep = hmi.Hmolec[ipMH2p]/plte;
	}
	else
	{
		hmi.h2pdep = 1.;
	}

	/* departure coef of H3+ defined rel to N(H2+) N(p) */
	if( hmi.rel_pop_LTE_H3p > 0. )
	{
		hmi.h3pdep = hmi.Hmolec[ipMH3p]/hmi.rel_pop_LTE_H3p;
	}
	else
	{
		hmi.h3pdep = 1.;
	}


	if( trace.lgTrace && trace.lgTr_H2_Mole )
	{
		fprintf( ioQQQ, " HMOLE, Dep Coef, H-:%10.2e H2:%10.2e H2+:%10.2e\n", 
		  hmi.hmidep, hmi.h2dep, hmi.h2pdep );
		fprintf( ioQQQ, "     H- creat: Rad atch%10.3e Induc%10.3e bHneut%10.2e 3bod%10.2e b=>H2%10.2e N(H-);%10.2e b(H-);%10.2e\n", 
		  hmi.hminus_rad_attach, hmi.HMinus_induc_rec_rate, bhneut, c3bod, assoc_detach_backwards_grnd, hmi.Hmolec[ipMHm], hmi.hmidep );

		fprintf( ioQQQ, "     H- destr: Photo;%10.3e mut neut%10.2e e- coll ion%10.2e =>H2%10.2e x-ray%10.2e p+H-%10.2e\n", 
		  hmi.HMinus_photo_rate, hmi.hmin_ct_firstions*sum_first_ions, cionhm, 
		  Hmolec_old[ipMH]*hmi.assoc_detach, iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s], 
		  fhneut );
		fprintf( ioQQQ, "     H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n", 
		  hmi.hmihet, hmi.HMinus_induc_rec_cooling, hmi.hmicol );
	}

	/* identify creation and destruction processes for H2+ */
	if( trace.lgTrace && trace.lgTr_H2_Mole )
	{
		rate = desh2p;
		if( rate != 0. )
		{
			fprintf( ioQQQ, 
				" Destroy H2+: rate=%10.2e e-;%5.3f phot;%5.3f hard gam;%5.3f H2col;%5.3f h2phhp;%5.3f pion;%5.3f bh2h2p:%5.3f\n", 
			  rate, h2pcin*dense.eden/rate, gamtwo/rate, 2.*iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/
			  rate, h2ph3p/rate, h2phhp/rate, h2pion/rate, hmi.bh2h2p*
			  dense.xIonDense[ipHYDROGEN][0]/rate );

			rate *= hmi.Hmolec[ipMH2p];
			if( rate > 0. )
			{
				fprintf( ioQQQ, 
					" Create  H2+: rate=%.2e HII HI;%.3f Col H2;%.3f HII H2;%.3f HI HI;%.3f\n", 
				  rate, 
				  radath*dense.xIonDense[ipHYDROGEN][1]*dense.xIonDense[ipHYDROGEN][0]/rate, 
				  (gamhd + secondaries.csupra[ipHYDROGEN][0]*2.02)*hmi.Hmolec[ipMH2g]/rate, 
				  rh2h2p*dense.xIonDense[ipHYDROGEN][1]*hmi.Hmolec[ipMH2g]/rate, 
				  b2pcin*dense.xIonDense[ipHYDROGEN][0]*dense.xIonDense[ipHYDROGEN][1]*dense.eden/rate );
			}
			else
			{
				fprintf( ioQQQ, " Create  H2+: rate= is zero\n" );
			}
		}
	}

	{
		/* following should be set true to print populations */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			fprintf(ioQQQ,"hmole bugg\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				fnzone,
				iso.gamnc[ipH_LIKE][ipHYDROGEN][0],
				hmi.HMinus_photo_rate,
				hmi.Hmolec[ipMH2g] , 
				hmi.Hmolec[ipMHm] ,
				dense.xIonDense[ipHYDROGEN][1]);
		}
	}

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