Context Navigation

← Previous Change
Next Change →

source

Timestamp:

12/10/07 12:09:42 (12 months ago)

Author:

rjrw

Message:

Further rationalization of chemistry network

Location:

branches/newmole/source

Files:

: 1 removed
: 18 modified

conv_base.cpp (modified) (1 diff)
cool_oxyg.cpp (modified) (1 diff)
hmi.h (modified) (1 diff)
init_sim_postparse.cpp (modified) (1 diff)
iso_ionize_recombine.cpp (modified) (1 diff)
mole.h (modified) (1 diff)
mole_co_drive.cpp (modified) (23 diffs)
mole_co_etc.cpp (modified) (23 diffs)
mole_co_solve.cpp (modified) (5 diffs)
mole_h2.cpp (modified) (1 diff)
mole_h2_form.cpp (modified) (3 diffs)
mole_h2_io.cpp (modified) (3 diffs)
mole_h_drive.cpp (deleted)
mole_newton_step.cpp (modified) (7 diffs)
mole_priv.h (modified) (1 diff)
mole_reactions.cpp (modified) (4 diffs)
prt_lines_lv1_li_ne.cpp (modified) (2 diffs)
prt_lines_molecules.cpp (modified) (1 diff)
punch_do.cpp (modified) (2 diffs)

Legend:

: Unmodified
: Added
: Removed

branches/newmole/source/conv_base.cpp

r1639	r1658
519	519
520	520	/* do carbon monoxide after oxygen */
521		CO_drive();
	521	mole_drive();
522	522
523	523	ASSERT(lgElemsConserved());

branches/newmole/source/cool_oxyg.cpp

r1610	r1658
192	192	* co.rate_OH_dissoc is rate OH -> O + H [cm-3 s-1],
193	193	* must make it per unit O atom, so this rate is s-1 excitations per O atom */
194		rate_OH_dissoc = CO_findrate("PHOTON,OH=>O,H");
	194	rate_OH_dissoc = mole_findrate("PHOTON,OH=>O,H");
195	195	r12 = rate_OH_dissoc*0.55/SDIV( dense.xIonDense[ipOXYGEN][0] );
196	196	r13 = rate_OH_dissoc*0.05/SDIV( dense.xIonDense[ipOXYGEN][0] );

branches/newmole/source/hmi.h

r1630	r1658
4	4	#ifndef _HMI_H_
5	5	#define _HMI_H_
6
7		~~/*hmole determine populations of hydrogen molecules /~~
8		~~double hmole(bool lgNegPop, bool lgZerPop);~~
9
10		~~/*hmole_reactions - evaluates hydrogen chemistry reactions /~~
11		~~void hmole_reactions(void);~~
12	6
13	7	/**hmirat computes radiative association rate for H-

branches/newmole/source/init_sim_postparse.cpp

r1610	r1658
24	24	thermal.lgUnstable = false;
25	25
26		CO_Init();
27		CO_zero();
	26	mole_Init();
	27	mole_zero();
28	28
29	29	/* >> chng 06 Nov 24 rjrw: Move reaction definitions here so they can be controlled by parsed commands */
30		CO_create_react();
	30	mole_create_react();
31	31
32	32	/* this term applies centrifugal acceleration if DISK option on wind command */

branches/newmole/source/iso_ionize_recombine.cpp

r1636	r1658
204	204
205	205	fixit(); /* Should be generic */
206		//iso.RateCont2Level[ipISO][nelem][ipHe1s1S] += (realnum) CO_sink_rate("He+");
	206	//iso.RateCont2Level[ipISO][nelem][ipHe1s1S] += (realnum) mole_sink_rate("He+");
207	207	/fprintf(ioQQQ,"DEBUG co hep destroy %e\n",co.hep_destroy );/
208	208

branches/newmole/source/mole.h

r1656	r1658
9	9	enum mole_state {MOLE_NULL, MOLE_PASSIVE, MOLE_ACTIVE};
10	10
11		/*CO_drive main driver for heavy molecular equilibrium routines /
12		extern void CO_drive(void);
13
14		/*CO_zero allocate + initialize workspace /
15		extern void CO_zero(void);
16
17		/*CO_create_react build reaction structures /
18		extern void CO_create_react(void);
	11	/*mole_drive main driver for heavy molecular equilibrium routines /
	12	extern void mole_drive(void);
	13
	14	/*mole_zero allocate + initialize workspace /
	15	extern void mole_zero(void);
	16
	17	/*mole_create_react build reaction structures /
	18	extern void mole_create_react(void);
19	19
20	20	/** called from cdInit to initialized co routines */
21		extern void CO_Init(void);
22		/*CO_update_rks update rate coefficients, only temp part /
23		extern void CO_update_rks( void );
24
25		extern void CO_update_species_cache(void);
26
27		extern void CO_return_cached_species(void);
28
29		extern double CO_sink_rate(const char chSpecies[]);
30
31		extern double CO_source_rate(const char chSpecies[]);
32
33		extern double CO_dissoc_rate(const char chSpecies[]);
34
35		extern double CO_chem_heat(void);
36
37		struct mole_reaction_s *CO_findrate_s(const char buf[]);
38
39		extern double CO_findrk(const char buf[]);
40
41		extern double CO_findrate(const char buf[]);
	21	extern void mole_Init(void);
	22	/*mole_update_rks update rate coefficients, only temp part /
	23	extern void mole_update_rks( void );
	24
	25	extern void mole_update_species_cache(void);
	26
	27	extern void mole_return_cached_species(void);
	28
	29	extern double mole_sink_rate(const char chSpecies[]);
	30
	31	extern double mole_source_rate(const char chSpecies[]);
	32
	33	extern double mole_dissoc_rate(const char chSpecies[]);
	34
	35	extern double mole_chem_heat(void);
	36
	37	struct mole_reaction_s *mole_findrate_s(const char buf[]);
	38
	39	extern double mole_findrk(const char buf[]);
	40
	41	extern double mole_findrate(const char buf[]);
42	42
43	43	/** Take one Newton step of the chemical network

branches/newmole/source/mole_co_drive.cpp

r1657	r1658
1	1	/* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland
2	2	* For conditions of distribution and use see copyright notice in license.txt */
3		/~~CO_drive - public routine, calls CO~~_solve to converge molecules /
	3	/mole_drive - public routine, calls mole_solve to converge molecules /
4	4	#include "cddefines.h"
5	5	#include "taulines.h"
…	…
7	7	#include "hmi.h"
8	8	#include "conv.h"
	9	#include "doppvel.h"
9	10	#include "gammas.h"
	11	#include "iso.h"
10	12	#include "mole_co_atom.h"
11	13	#include "opacity.h"
…	…
13	15	#include "physconst.h"
14	16	#include "radius.h"
	17	#include "rfield.h"
15	18	#include "secondaries.h"
16	19	#include "timesc.h"
…	…
42	45	STATIC void mole_effects(void);
43	46
44		static void CO_null_results(void);
45		void grain_form(void);
46		void mole_h_rate_diagnostics(void);
	47	static void mole_null_results(void);
	48	STATIC void mole_h_rate_diagnostics(void);
47	49	void mole_h_fixup(void);
48	50
…	…
50	52
51	53
52		/CO_drive main driver for heavy molecular equilibrium routines /
53		void CO_drive(void)
	54	/mole_drive main driver for heavy molecular equilibrium routines /
	55	void mole_drive(void)
54	56	{
55	57	bool lgNegPop,
…	…
64	66	error;
65	67
66		DEBUG_ENTRY( "CO_drive()" );
	68	DEBUG_ENTRY( "mole_drive()" );
67	69
68	70	/* h2lim is smallest ratio of H2/HDEN for which we will
…	…
84	86	/* we should have an H0 soln at this point
85	87	-- all species for network will be sourced from ionization solvers
86		in CO_update_species_cache below */
	88	in mole_update_species_cache below */
87	89	ASSERT( dense.xIonDense[ipHYDROGEN][0]>SMALLFLOAT );
88	90	}
…	…
136	138	MolecSetup(chReason);
137	139
138		CO_update_species_cache(); /* Update densities of species controlled outside the chemical network */
	140	mole_update_species_cache(); /* Update densities of species controlled outside the chemical network */
139	141
140		grain_form();
141
142		CO_update_rks();
	142	mole_update_rks();
143	143
144	144	solomon_assumed = hmi.H2_Solomon_dissoc_rate_used_H2g/SDIV(hmi.H2_rate_destroy);
…	…
190	190	* substantially below the threshold for solving the CO network.
191	191	* Turn it off */
192		/* >> chng 07 jan 10 rjrw: this was CO_Init(), but the comment suggests
193		* it should really be CO_zero */
194		CO_zero();
	192	/* >> chng 07 jan 10 rjrw: this was mole_Init(), but the comment suggests
	193	* it should really be mole_zero */
	194	mole_zero();
195	195	}
196	196	else
…	…
200	200	conv.lgConvPops = false;
201	201	fprintf( ioQQQ,
202		" PROBLEM CO_drive failed to converge1 due to negative pops, zone=%.2f, CO/C=%.2e OH/C=%.2e H2/H=%.2e\n",
	202	" PROBLEM mole_drive failed to converge1 due to negative pops, zone=%.2f, CO/C=%.2e OH/C=%.2e H2/H=%.2e\n",
203	203	fnzone,
204	204	findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
…	…
233	233	conv.BadConvIoniz[1] = sum_mol;
234	234	if( CODEBUG )
235		fprintf(ioQQQ,"CO_drive not converged\n");
	235	fprintf(ioQQQ,"mole_drive not converged\n");
236	236	}
237	237	sum_ion = dense.xIonDense[ipOXYGEN][0] + dense.xIonDense[ipOXYGEN][1];
…	…
247	247	conv.BadConvIoniz[1] = sum_mol;
248	248	if( CODEBUG )
249		fprintf(ioQQQ,"CO_drive not converged\n");
	249	fprintf(ioQQQ,"mole_drive not converged\n");
250	250	}
251	251	}
…	…
267	267
268	268
269		/* CO_Init set this to -2 when code initialized, so negative
	269	/* mole_Init set this to -2 when code initialized, so negative
270	270	* number shows very first call in this model */
271	271	/* >>chng 05 jul 15, do not init if we have never evaluated CO in this iteration
…	…
341	341	}
342	342
343		static void CO_null_results(void)
	343	static void mole_null_results(void)
344	344	{
345		DEBUG_ENTRY( "CO_null_results()" );
	345	DEBUG_ENTRY( "mole_null_results()" );
346	346	/* these are heavy - heavy charge transfer rate that will still be needed
347	347	* for recombination of Si+, S+, and C+ */
…	…
402	402	/* trace ionization */
403	403	fprintf( ioQQQ,
404		" CO_drive4 do not evaluate CO chemistry.\n");
	404	" mole_drive4 do not evaluate CO chemistry.\n");
405	405	}
406	406
…	…
430	430	mole_eval_balance(mole.num_total,b,c);
431	431
432		co.CODissHeat = (realnum)(CO_findrate("PHOTON,CO=>C,O")*1e-12);
	432	co.CODissHeat = (realnum)(mole_findrate("PHOTON,CO=>C,O")*1e-12);
433	433
434	434	thermal.heating[0][9] = co.CODissHeat;
…	…
490	490	/* local timescale for H2 formation
491	491	* both grains and radiative attachment */
492		timesc.time_H2_Form_here = (~~CO_findrk("H,H,grnh=>H2")+CO~~_findrk("H,H,grnh=>H2"))
	492	timesc.time_H2_Form_here = (mole_findrk("H,H,grnh=>H2")+mole_findrk("H,H,grnh=>H2"))
493	493	/* this corrects for fact that we the timescale for H2 to form from an atomic gas.
494	494	* The rate becomes very small when gas is fully molecular, and ratio of total hydrogen
…	…
550	550
551	551	/* TH85 deexcitation heating */
552		hmi.HeatH2Dexc_TH85 = ((~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))
553		- (~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))) * 4.17e-12;
	552	hmi.HeatH2Dexc_TH85 = ((mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))
	553	- (mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))) * 4.17e-12;
554	554	/* this is derivative wrt temperature, only if counted as a cooling source */
555	555	hmi.deriv_HeatH2Dexc_TH85 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_TH85)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
…	…
563	563
564	564	/* Burton et al. 1990 heating */
565		hmi.HeatH2Dexc_BHT90 = ((~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))
566		- (~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))) * 4.17e-12;
	565	hmi.HeatH2Dexc_BHT90 = ((mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))
	566	- (mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))) * 4.17e-12;
567	567	/* this is derivative wrt temperature, only if counted as a cooling source */
568	568	hmi.deriv_HeatH2Dexc_BHT90 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BHT90)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
…	…
575	575	hmi.H2_Solomon_dissoc_rate_used_H2s * findspecies("H2*")->hevmol);
576	576	/* Bertoldi & Draine heating, same as TH85 */
577		hmi.HeatH2Dexc_BD96 = ((~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))
578		- (~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))) * 4.17e-12;
	577	hmi.HeatH2Dexc_BD96 = ((mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))
	578	- (mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))) * 4.17e-12;
579	579	/* this is derivative wrt temperature, only if counted as a cooling source */
580	580	hmi.deriv_HeatH2Dexc_BD96 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BD96)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
…	…
682	682	f5 = MAX2(1.,0.95 * log_density - 1.45) * 0.2 * log_G0_face;
683	683
684		hmi.HeatH2Dexc_ELWERT = ((~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))
685		- (~~CO_findrate("H2,H2=>H2,H2") + CO~~_findrate("H2,H=>H2,H"))) * 4.17e-12 * f3 +
	684	hmi.HeatH2Dexc_ELWERT = ((mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))
	685	- (mole_findrate("H2,H2=>H2,H2") + mole_findrate("H2,H=>H2,H"))) * 4.17e-12 * f3 +
686	686	f4 * (COmole[element_list[ipHYDROGEN]->ipMl[0]]->hevmol/dense.gas_phase[ipHYDROGEN]) +
687	687	f5 * secondaries.x12tot * EN1EV * hmi.H2_total;
…	…
735	735	* more terms added */
736	736	hmi.H2_rate_create =
737		CO_findrate("H,H,grnh=>H2*")+
738		CO_findrate("H,H,grnh=>H2") +
739		CO_findrate("H,H-=>H2,e-")+
740		CO_findrate("H,H-=>H2*,e-") +
741		CO_findrate("H,H,H=>H2,H") +
742		CO_findrate("H,H2+=>H+,H2*") +
743		CO_findrate("H,H=>H2") +
744		CO_findrate("H,H3+=>H2,H2+") +
745		CO_findrate("H2+,H-=>H2,H") +
746		CO_findrate("H,H,H2=>H2,H2") +
747		CO_findrate("H3+,H-=>H2,H,H") +
748		CO_findrate("H3+,H-=>H2,H2") +
749		CO_findrate("H2,H3+=>H+,H2,H2") +
750		CO_findrate("e-,H3+=>H2,H") +
751		CO_findrate("H3+,PHOTON=>H2,H+") +
752		CO_findrate("H,CH=>C,H2") +
753		CO_findrate("H,CH+=>C+,H2") +
754		CO_findrate("H,CH2=>CH,H2") +
755		CO_findrate("H,CH3+=>CH2+,H2") +
756		CO_findrate("H,OH=>O,H2") +
757		CO_findrate("H-,HCO+=>CO,H2") +
758		CO_findrate("H-,H3O+=>H2O,H2") +
759		CO_findrate("H-,H3O+=>OH,H2,H") +
760		CO_findrate("H+,CH2=>CH+,H2") +
761		CO_findrate("H+,SiH=>Si+,H2") +
762		CO_findrate("H2+,CH=>CH+,H2") +
763		CO_findrate("H2+,CH2=>CH2+,H2") +
764		CO_findrate("H2+,CO=>CO+,H2") +
765		CO_findrate("H2+,H2O=>H2O+,H2") +
766		CO_findrate("H2+,O2=>O2+,H2") +
767		CO_findrate("H2+,OH=>OH+,H2") +
768		CO_findrate("H3+,C=>CH+,H2") +
769		CO_findrate("H3+,CH=>CH2+,H2") +
770		CO_findrate("H3+,CH2=>CH3+,H2") +
771		CO_findrate("H3+,OH=>H2O+,H2") +
772		CO_findrate("H3+,H2O=>H3O+,H2") +
773		CO_findrate("H3+,CO=>HCO+,H2") +
774		CO_findrate("H3+,O=>OH+,H2") +
775		CO_findrate("H3+,SiH=>SiH2+,H2") +
776		CO_findrate("H3+,SiO=>SiOH+,H2") +
777		CO_findrate("H,CH3=>CH2,H2") +
778		CO_findrate("H,CH4+=>CH3+,H2") +
779		CO_findrate("H,CH5+=>CH4+,H2") +
780		CO_findrate("H2+,CH4=>CH3+,H2,H") +
781		CO_findrate("H2+,CH4=>CH4+,H2") +
782		CO_findrate("H3+,CH3=>CH4+,H2") +
783		CO_findrate("H3+,CH4=>CH5+,H2") +
784		CO_findrate("H+,CH4=>CH3+,H2") +
785		CO_findrate("H+,HNO=>NO+,H2") +
786		CO_findrate("H+,HS=>S+,H2") +
787		CO_findrate("H,HS+=>S+,H2") +
788		CO_findrate("H3+,NH=>NH2+,H2") +
789		CO_findrate("H3+,NH2=>NH3+,H2") +
790		CO_findrate("H3+,NH3=>NH4+,H2") +
791		CO_findrate("H3+,CN=>HCN+,H2") +
792		CO_findrate("H3+,NO=>HNO+,H2") +
793		CO_findrate("H3+,S=>HS+,H2") +
794		CO_findrate("H3+,CS=>HCS+,H2") +
795		CO_findrate("H3+,NO2=>NO+,OH,H2") +
796		CO_findrate("H2+,NH=>NH+,H2") +
797		CO_findrate("H2+,NH2=>NH2+,H2") +
798		CO_findrate("H2+,NH3=>NH3+,H2") +
799		CO_findrate("H2+,CN=>CN+,H2") +
800		CO_findrate("H2+,HCN=>HCN+,H2") +
801		CO_findrate("H2+,NO=>NO+,H2") +
802		CO_findrate("H3+,Cl=>HCl+,H2")+
803		CO_findrate("H3+,HCl=>H2Cl+,H2")+
804		CO_findrate("H2+,C2=>C2+,H2")+
805		CO_findrate("H-,NH4+=>NH3,H2")+
806		CO_findrate("H3+,HCN=>HCNH+,H2");
	737	mole_findrate("H,H,grnh=>H2*")+
	738	mole_findrate("H,H,grnh=>H2") +
	739	mole_findrate("H,H-=>H2,e-")+
	740	mole_findrate("H,H-=>H2*,e-") +
	741	mole_findrate("H,H,H=>H2,H") +
	742	mole_findrate("H,H2+=>H+,H2*") +
	743	mole_findrate("H,H=>H2") +
	744	mole_findrate("H,H3+=>H2,H2+") +
	745	mole_findrate("H2+,H-=>H2,H") +
	746	mole_findrate("H,H,H2=>H2,H2") +
	747	mole_findrate("H3+,H-=>H2,H,H") +
	748	mole_findrate("H3+,H-=>H2,H2") +
	749	mole_findrate("H2,H3+=>H+,H2,H2") +
	750	mole_findrate("e-,H3+=>H2,H") +
	751	mole_findrate("H3+,PHOTON=>H2,H+") +
	752	mole_findrate("H,CH=>C,H2") +
	753	mole_findrate("H,CH+=>C+,H2") +
	754	mole_findrate("H,CH2=>CH,H2") +
	755	mole_findrate("H,CH3+=>CH2+,H2") +
	756	mole_findrate("H,OH=>O,H2") +
	757	mole_findrate("H-,HCO+=>CO,H2") +
	758	mole_findrate("H-,H3O+=>H2O,H2") +
	759	mole_findrate("H-,H3O+=>OH,H2,H") +
	760	mole_findrate("H+,CH2=>CH+,H2") +
	761	mole_findrate("H+,SiH=>Si+,H2") +
	762	mole_findrate("H2+,CH=>CH+,H2") +
	763	mole_findrate("H2+,CH2=>CH2+,H2") +
	764	mole_findrate("H2+,CO=>CO+,H2") +
	765	mole_findrate("H2+,H2O=>H2O+,H2") +
	766	mole_findrate("H2+,O2=>O2+,H2") +
	767	mole_findrate("H2+,OH=>OH+,H2") +
	768	mole_findrate("H3+,C=>CH+,H2") +
	769	mole_findrate("H3+,CH=>CH2+,H2") +
	770	mole_findrate("H3+,CH2=>CH3+,H2") +
	771	mole_findrate("H3+,OH=>H2O+,H2") +
	772	mole_findrate("H3+,H2O=>H3O+,H2") +
	773	mole_findrate("H3+,CO=>HCO+,H2") +
	774	mole_findrate("H3+,O=>OH+,H2") +
	775	mole_findrate("H3+,SiH=>SiH2+,H2") +
	776	mole_findrate("H3+,SiO=>SiOH+,H2") +
	777	mole_findrate("H,CH3=>CH2,H2") +
	778	mole_findrate("H,CH4+=>CH3+,H2") +
	779	mole_findrate("H,CH5+=>CH4+,H2") +
	780	mole_findrate("H2+,CH4=>CH3+,H2,H") +
	781	mole_findrate("H2+,CH4=>CH4+,H2") +
	782	mole_findrate("H3+,CH3=>CH4+,H2") +
	783	mole_findrate("H3+,CH4=>CH5+,H2") +
	784	mole_findrate("H+,CH4=>CH3+,H2") +
	785	mole_findrate("H+,HNO=>NO+,H2") +
	786	mole_findrate("H+,HS=>S+,H2") +
	787	mole_findrate("H,HS+=>S+,H2") +
	788	mole_findrate("H3+,NH=>NH2+,H2") +
	789	mole_findrate("H3+,NH2=>NH3+,H2") +
	790	mole_findrate("H3+,NH3=>NH4+,H2") +
	791	mole_findrate("H3+,CN=>HCN+,H2") +
	792	mole_findrate("H3+,NO=>HNO+,H2") +
	793	mole_findrate("H3+,S=>HS+,H2") +
	794	mole_findrate("H3+,CS=>HCS+,H2") +
	795	mole_findrate("H3+,NO2=>NO+,OH,H2") +
	796	mole_findrate("H2+,NH=>NH+,H2") +
	797	mole_findrate("H2+,NH2=>NH2+,H2") +
	798	mole_findrate("H2+,NH3=>NH3+,H2") +
	799	mole_findrate("H2+,CN=>CN+,H2") +
	800	mole_findrate("H2+,HCN=>HCN+,H2") +
	801	mole_findrate("H2+,NO=>NO+,H2") +
	802	mole_findrate("H3+,Cl=>HCl+,H2")+
	803	mole_findrate("H3+,HCl=>H2Cl+,H2")+
	804	mole_findrate("H2+,C2=>C2+,H2")+
	805	mole_findrate("H-,NH4+=>NH3,H2")+
	806	mole_findrate("H3+,HCN=>HCNH+,H2");
807	807
808	808	/>>chng 05 jul 01, GS, h2s ionization by cosmic ray added/
809	809	/* >>chng 05 jun 29, TE, used in new punch H2 destruction file*/
810		hmi.CR_reac_H2g = ~~CO_findrk("H2,CRPHOT=>H2+,e-") + CO~~_findrk("H2,CRPHOT=>H,H+,e-")
811		+ ~~CO_findrk("H2,CRPHOT=>H,H") + CO_findrk("H2,CRPHOT=>H+,H-") + CO~~_findrk("H2,CRPHOT=>H+,H,e-");
812		hmi.CR_reac_H2s = ~~CO_findrk("H2,CRPHOT=>H2+,e-") + CO~~_findrk("H2,CRPHOT=>H,H+,e-") +
813		~~CO_findrk("H2,CRPHOT=>H,H") + CO_findrk("H2,CRPHOT=>H+,H-") + CO~~_findrk("H2*,CRPHOT=>H+,H,e-");
	810	hmi.CR_reac_H2g = mole_findrk("H2,CRPHOT=>H2+,e-") + mole_findrk("H2,CRPHOT=>H,H+,e-")
	811	+ mole_findrk("H2,CRPHOT=>H,H") + mole_findrk("H2,CRPHOT=>H+,H-") + mole_findrk("H2,CRPHOT=>H+,H,e-");
	812	hmi.CR_reac_H2s = mole_findrk("H2,CRPHOT=>H2+,e-") + mole_findrk("H2,CRPHOT=>H,H+,e-") +
	813	mole_findrk("H2,CRPHOT=>H,H") + mole_findrk("H2,CRPHOT=>H+,H-") + mole_findrk("H2*,CRPHOT=>H+,H,e-");
814	814
815	815	if( findspecies("H-")->hevmol > 0. && hmi.rel_pop_LTE_Hmin > 0. )
…	…
869	869	return;
870	870	}
	871	#if defined(__HP_aCC)
	872	#pragma OPTIMIZE OFF
	873	#pragma OPTIMIZE ON
	874	#endif
	875	STATIC void mole_h_rate_diagnostics(void)
	876	{
	877	int i;
	878	double
	879	frac_H2star_hminus(),
	880	rate;
	881
	882	/* identify dominant H2 formation process */
	883	{
	884	/* following should be set true to identify H2 formation and destruction processes */
	885	/@-redef@/
	886	enum {DEBUG_LOC=false};
	887	/@+redef@/
	888	if( DEBUG_LOC && (nzone>50) )
	889	{
	890	double createsum ,create_from_Hn2 , create_3body_Ho, create_h2p,
	891	create_h3p, create_h3pe, create_grains, create_hminus;
	892	double destroysum, destroy_hm ,destroy_solomon ,destroy_2h ,destroy_hp,
	893	destroy_h,destroy_hp2,destroy_h3p;
	894
	895	create_from_Hn2 = mole_findrate("H,H=>H2"); /* H(n=2) + H(n=1) -> H2 */
	896	create_3body_Ho = mole_findrate("H,H,H=>H2,H");
	897	create_h2p = mole_findrate("H,H2+=>H+,H2*");
	898	create_h3p = mole_findrate("H,H3+=>H2,H2+");
	899	create_h3pe = mole_findrate("e-,H3+=>H2,H");
	900	create_grains = mole_findrk("H,H,grnh=>H2*")+mole_findrk("H,H,grnh=>H2");
	901	create_hminus = (mole_findrate("H,H-=>H2,e-")+mole_findrate("H,H-=>H2*,e-"));
	902
	903	createsum = create_from_Hn2 + create_3body_Ho + create_h2p +
	904	create_h3p + create_h3pe + create_grains + create_hminus;
	905
	906	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",
	907	fnzone,
	908	createsum,
	909	create_hminus / createsum,
	910	create_from_Hn2 / createsum,
	911	create_3body_Ho / createsum,
	912	create_h2p / createsum,
	913	create_h3p / createsum,
	914	create_h3pe / createsum,
	915	create_grains / createsum );
	916
	917	/* all h2 molecular hydrogen destruction processes */
	918	/* >>chng 04 jan 28, had wrong Boltzmann factor,
	919	* caught by gargi Shaw */
	920	destroy_hm = (mole_findrk("H2,e-=>H,H-")+mole_findrk("H2,e-=>H,H-"))findspecies("e-")->hevmol;
	921	/destroy_hm2 = eh2hhm;/
	922	destroy_solomon = hmi.H2_Solomon_dissoc_rate_used_H2g;
	923	destroy_2h = (mole_findrk("H2,e-=>H,H,e-")+mole_findrk("H2*,e-=>H,H,e-"));
	924	destroy_hp = mole_findrk("H2,H+=>H3+")*findspecies("H+")->hevmol;
	925	destroy_h = mole_findrk("