/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*CO_photo fills in matrix elements for heavy element photon and cosmic ray processes */
#include "cddefines.h"
#include "secondaries.h"
#include "ionbal.h"
#include "grainvar.h"
#include "hmi.h"
#include "rfield.h"
#include "phycon.h"
#include "mole.h"
/*lint -e778 constant expression evaluatess to 0 in operation '-' */
/*=================================================================*/

/*CO_step fills in matrix for heavy elements molecular routines */
void CO_photo(void )
{

	double
	C2H2_CRPHOT_C2H2P_e,
	C2H2_CRPHOT_C2H_H,
	C3H_CRPHOT_C3_H,
	C2H2_PHOTON_C2H2P_e,
	C2H2_PHOTON_C2H_H,
	C3H_PHOTON_C3_H,
	C3_CRPHOT_C2_C,
	C3_PHOTON_C2_C,
	C2H_CRPHOT_C2HP_e,
	C2H_PHOTON_C2HP_e,
	C2H_PHOTON_C2_H,
	C2H_CRPHOT_C2_H,
	crnu_CH_C_H,
	crnu_CHP_CP_H,
	crnu_H2O_OH_H,
	crnu_O2_O2P_e,
	crnu_O2_O_O,
	crnu_OH_O_H,
	crnu_SiH_Si_H,
	crnu_SiO_Si_O,
	nu_CH_CHP_e,
	nu_CH_C_H,
	nu_CHP_CP_H,
	nu_CH2_CH2P_e,
	nu_CH2_CH_H,
	nu_CH2P_CHP_H,
	nu_CH3P_CH2P_H,
	nu_CH3P_CHP_H2,
	nu_H2O_H2OP_e,
	nu_H2O_OH_H,
	nu_O2_O2P_e,
	nu_O2_O_O,
	nu_OH_OHP_e,
	nu_OH_O_H,
	nu_OHP_O_HP,
	nu_SiH_Si_H,
	nu_SiO_Si_O,
	nu_CH3_CH3P_e,
	nu_CH3_CH2_H,
	nu_CH3_CH_H2,
	crnu_CH3_CH3P_e,
	crnu_CH3_CH2_H,
	crnu_CH3_CH_H2,
	nu_CH4_CH3_H,
	nu_CH4_CH2_H2,
	nu_CH4_CH_H2,
	crnu_CH4_CH2_H2,
	NH_NU_NHP_e   ,
	NH_NU_N_H ,
	NH2_NU_NH_H   ,
	NH2_NU_NH2P_e ,
	NH3_NU_NH3P_e ,
	NH3_NU_NH_H2  ,
	NH3_NU_NH2_H  ,
	CN_NU_N_C ,
	HCN_NU_CN_H   ,
	N2_NU_N_N ,
	NO_NU_O_N ,
	NO_NU_NOP_e   ,
	HNO_NU_NO_H   ,
	HS_NU_S_H ,
	HSP_NU_SP_H   ,
	OCN_NU_CN_O   ,
	CS_NU_S_C ,
	CSP_NU_SP_C   ,
	CS_NU_CSP_e   ,
	NO2_NU_NO_O   ,
	NS_NU_S_N ,
	SO_NU_S_O ,
	SO_NU_SOP_e   ,
	OCS_NU_OCSP_e ,
	NH_CRP_NHP_e,
	NH_CRP_N_H   ,     
	NH2_CRP_NH_H  ,    
	NH2_CRP_NH2P_e ,   
	NH3_CRP_NH3P_e  ,  
	NH3_CRP_NH_H2    , 
	NH3_CRP_NH2_H     ,
	CN_CRP_N_C        ,
	HCN_CRP_CN_H      ,
	N2_CRP_N_N        ,
	NO_CRP_O_N        ,
	NO_CRP_NOP_e      ,
	HNO_CRP_NO_H      ,
	HS_CRP_S_H        ,
	OCN_CRP_CN_O      ,
	CS_CRP_S_C        ,
	CS_CRP_CSP_e      ,
	NO2_CRP_NO_O      ,
	NS_CRP_S_N        ,
	SO_CRP_S_O        ,
	SO_CRP_SOP_e      ,
	OCS_CRP_OCSP_e    ,
	C2_nu_C_C,
	C2_nu_C2P_e,
	C2_crnu_C_C,
	C2P_nu_CP_C,
	HCl_NU_Cl_H,
	CCl_NU_Cl_C,
	ClO_NU_Cl_O,
	HCl_crnu_Cl_H,
	CCl_crnu_Cl_C,
	ClO_crnu_Cl_O,
	HNC_PHOTON_CN_H   ,
	HNC_CRPHOT_CN_H   ,
	OH,
	crnu_rate, 
	TH85_field;


	TH85_field = hmi.UV_Cont_rel2_Habing_TH85_depth/1.66f;
	OH = co.hevmol[ipOH];

	/* >> chng aug 24, 05 NPA This is the cosmic ray ionization rate used in the molecular network.  TH85 and the e-mail from Amiel Sternberg has each cosmic ray ionization rate as a leading coefficient multiplied by a scale factor.
	The leading coefficient is defined as the cosmic ray rate for H + CR = > H+
	+ e- .  For molecules in the heavy element molecular network, this scale
	factor is derived by taking the rate for:

		X + CRPHOT => Y + Z 
		and dividing it by the rate:
		H + CRP => H+ + e-

	This scale factor is 2.17 for all cosmic ray reactions in this network 
			crnu_rate = secondaries.csupra[ipHYDROGEN][0];*/

	crnu_rate = 2.17*secondaries.csupra[ipHYDROGEN][0];


	
	/* >>chng 02 jul 04 -- The following are dissociation rates for various molecular species
	For right now we will calculate this rate by the standard form:
	   
			(alpha)*Go*exp(-Beta*AV)


    when the command "set Leiden hack UMIST rates" is used.  Otherwise we
	will just let cloudy calculate the value of the UV radiation field */

	if (!co.lgUMISTrates)
	{

		nu_CH_CHP_e             = HMRATE((1e-10)*7.6,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.8*rfield.extin_mag_V_point))/1.66);
		nu_CH_C_H               = HMRATE((1e-10)*8.6,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.2*rfield.extin_mag_V_point))/1.66);
		nu_CHP_CP_H             = HMRATE((1e-10)*2.5,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.5*rfield.extin_mag_V_point))/1.66);
		nu_CH2_CH2P_e           = HMRATE(1.0e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.3*rfield.extin_mag_V_point))/1.66);
		nu_CH2_CH_H             = HMRATE((1e-10)*7.2,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.7*rfield.extin_mag_V_point))/1.66);
		nu_CH2P_CHP_H           = HMRATE((1e-9)*1.7,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.7*rfield.extin_mag_V_point))/1.66);
		nu_CH3P_CH2P_H          = HMRATE(1.0e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.7*rfield.extin_mag_V_point))/1.66);
		nu_CH3P_CHP_H2          = HMRATE(1.0e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.7*rfield.extin_mag_V_point))/1.66);
		nu_H2O_H2OP_e           = HMRATE(3.3e-11,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(3.9*rfield.extin_mag_V_point))/1.66);
		nu_H2O_OH_H             = HMRATE((1e-10)*5.9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.7*rfield.extin_mag_V_point))/1.66);
		nu_O2_O2P_e             = HMRATE((1e-11)*5.6,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(3.7*rfield.extin_mag_V_point))/1.66);
		nu_O2_O_O               = HMRATE((1e-10)*6.9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.8*rfield.extin_mag_V_point))/1.66);
		nu_OH_OHP_e             = HMRATE((1e-12)*1.6,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(3.1*rfield.extin_mag_V_point))/1.66);
		nu_OH_O_H               = HMRATE(3.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.7*rfield.extin_mag_V_point))/1.66);
		nu_OHP_O_HP             = HMRATE((1e-12)*1.,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.8*rfield.extin_mag_V_point))/1.66);
		nu_SiH_Si_H             = HMRATE((1e-9)*2.8,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.1*rfield.extin_mag_V_point))/1.66);
		nu_SiO_Si_O             = HMRATE(1.0e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.3*rfield.extin_mag_V_point))/1.66);
		nu_CH3_CH3P_e           = HMRATE(1.0e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.1*rfield.extin_mag_V_point))/1.66);
		nu_CH3_CH2_H            = HMRATE((1e-10)*2.5,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.9*rfield.extin_mag_V_point))/1.66);
		nu_CH3_CH_H2            = HMRATE((1e-10)*2.5,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(1.9*rfield.extin_mag_V_point))/1.66);
		nu_CH4_CH3_H            = HMRATE((1e-10)*2.2,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.2*rfield.extin_mag_V_point))/1.66);
		nu_CH4_CH2_H2           = HMRATE((1e-10)*9.8,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.2*rfield.extin_mag_V_point))/1.66);
		nu_CH4_CH_H2            = HMRATE((1e-10)*2.2,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*(float)exp(-(2.2*rfield.extin_mag_V_point))/1.66);
		
	}
	/*  This is what will usually be done, calculating the photodissociation rate by calculating the UV
	 * radiation field self-consistently. */
	else
	{
		
		nu_CH_CHP_e             = HMRATE((1e-10)*7.6,0,0)*(TH85_field);
		nu_CH_C_H               = HMRATE((1e-10)*8.6,0,0)*(TH85_field);
		nu_CHP_CP_H             = HMRATE((1e-10)*2.5,0,0)*(TH85_field);
		nu_CH2_CH2P_e           = HMRATE(1.0e-9,0,0)*(TH85_field);
		nu_CH2_CH_H             = HMRATE((1e-10)*7.2,0,0)*(TH85_field);
		nu_CH2P_CHP_H           = HMRATE((1e-9)*1.7,0,0)*(TH85_field);
		nu_CH3P_CH2P_H          = HMRATE(1.0e-9,0,0)*(TH85_field);
		nu_CH3P_CHP_H2          = HMRATE(1.0e-9,0,0)*(TH85_field);
		nu_H2O_H2OP_e           = HMRATE(3.3e-11,0,0)*(TH85_field);
		nu_H2O_OH_H             = HMRATE((1e-10)*5.9,0,0)*(TH85_field);
		nu_O2_O2P_e             = HMRATE((1e-11)*5.6,0,0)*(TH85_field);
		nu_O2_O_O               = HMRATE((1e-10)*6.9,0,0)*(TH85_field);
		nu_OH_OHP_e             = HMRATE((1e-12)*1.6,0,0)*(TH85_field);
		nu_OH_O_H               = HMRATE(3.5e-10,0,0)*(TH85_field);
		nu_OHP_O_HP             = HMRATE((1e-12)*1.,0,0)*(TH85_field);
		nu_SiH_Si_H             = HMRATE((1e-9)*2.8,0,0)*(TH85_field);
		nu_SiO_Si_O             = HMRATE(1.0e-10,0,0)*(TH85_field);
		nu_CH3_CH3P_e           = HMRATE(1.0e-10,0,0)*(TH85_field);
		nu_CH3_CH2_H            = HMRATE((1e-10)*2.5,0,0)*(TH85_field);
		nu_CH3_CH_H2            = HMRATE((1e-10)*2.5,0,0)*(TH85_field);
		nu_CH4_CH3_H            = HMRATE((1e-10)*2.2,0,0)*(TH85_field);
		nu_CH4_CH2_H2           = HMRATE((1e-10)*9.8,0,0)*(TH85_field);
		nu_CH4_CH_H2            = HMRATE((1e-10)*2.2,0,0)*(TH85_field);
	}
	NH_NU_NHP_e       = HMRATE(1.0e-11,0,0)*(TH85_field);
	NH_NU_N_H         = HMRATE(5.0e-10,0,0)*(TH85_field);
	NH2_NU_NH_H       = HMRATE((1e-10)*3.9,0,0)*(TH85_field);
	NH2_NU_NH2P_e     = HMRATE((1e-10)*1.73,0,0)*(TH85_field);
	NH3_NU_NH3P_e     = HMRATE((1e-10)*1.8,0,0)*(TH85_field);
	NH3_NU_NH_H2      = HMRATE((1e-10)*3.3,0,0)*(TH85_field);
	NH3_NU_NH2_H      = HMRATE(1.1e-9,0,0)*(TH85_field);
	CN_NU_N_C         = HMRATE(1.1e-9,0,0)*(TH85_field);
	HCN_NU_CN_H       = HMRATE(1.3e-9,0,0)*(TH85_field);
	N2_NU_N_N         = HMRATE((1e-10)*2.3,0,0)*(TH85_field);
	NO_NU_O_N         = HMRATE((1e-10)*4.3,0,0)*(TH85_field);
	NO_NU_NOP_e       = HMRATE((1e-10)*2.6,0,0)*(TH85_field);
	HNO_NU_NO_H       = HMRATE((1e-10)*1.7,0,0)*(TH85_field);
	HS_NU_S_H         = HMRATE((1e-10)*9.7,0,0)*(TH85_field);
	HSP_NU_SP_H       = HMRATE((1e-10)*3.,0,0)*(TH85_field);
	OCN_NU_CN_O       = HMRATE(1.0e-11,0,0)*(TH85_field);
	CS_NU_S_C         = HMRATE((1e-10)*9.7,0,0)*(TH85_field);
	CSP_NU_SP_C       = HMRATE(2.0e-10,0,0)*(TH85_field);
	CS_NU_CSP_e       = HMRATE(2.0e-10,0,0)*(TH85_field);
	NO2_NU_NO_O       = HMRATE(1.2e-9,0,0)*(TH85_field);
	NS_NU_S_N         = HMRATE(1.0e-11,0,0)*(TH85_field);
	SO_NU_S_O         = HMRATE((1e-9)*3.7,0,0)*(TH85_field);
	SO_NU_SOP_e       = HMRATE((1e-10)*6.,0,0)*(TH85_field);
	OCS_NU_OCSP_e     = HMRATE((1e-10)*4.2,0,0)*(TH85_field);
	NH_CRP_NHP_e      = crnu_rate * 500 * 2;
	NH_CRP_N_H        = crnu_rate * 500 * 2;
	NH2_CRP_NH_H      = crnu_rate * 81 * 2;
	NH2_CRP_NH2P_e    = crnu_rate * 610 * 2;
	NH3_CRP_NH3P_e    = crnu_rate * 543 * 2;
	NH3_CRP_NH_H2     = crnu_rate * 541 * 2;
	NH3_CRP_NH2_H     = crnu_rate * 1325 * 2;
	CN_CRP_N_C        = crnu_rate * 10580 * 2;
	HCN_CRP_CN_H      = crnu_rate * 3114 * 2;
	N2_CRP_N_N        = crnu_rate * 50 * 2;
	NO_CRP_O_N        = crnu_rate * 427 * 2;
	NO_CRP_NOP_e      = crnu_rate * 430 * 2;
	HNO_CRP_NO_H      = crnu_rate * 1000 * 2;
	HS_CRP_S_H        = crnu_rate * 500 * 2;
	OCN_CRP_CN_O      = crnu_rate * 1500 * 2;
	CS_CRP_S_C        = crnu_rate * 500 * 2;
	CS_CRP_CSP_e      = crnu_rate * 500 * 2;
	NO2_CRP_NO_O      = crnu_rate * 1500 * 2;
	NS_CRP_S_N        = crnu_rate * 500 * 2;
	SO_CRP_S_O        = crnu_rate * 500 * 2;
	SO_CRP_SOP_e      = crnu_rate * 500 * 2;
	OCS_CRP_OCSP_e    = crnu_rate * 1444 * 2;
	HCl_NU_Cl_H = HMRATE(0.00000000011,0,0)*TH85_field;
	CCl_NU_Cl_C = HMRATE(0.0000000001,0,0)*TH85_field;
	ClO_NU_Cl_O = HMRATE(0.0000000001,0,0)*TH85_field;
	HCl_crnu_Cl_H = crnu_rate*2*610;
	CCl_crnu_Cl_C = crnu_rate*2*500;
	ClO_crnu_Cl_O = crnu_rate*2*500;
	nu_CH_CHP_e = HMRATE((1e-10)*4.2,0,0)*(TH85_field);
	nu_CH_C_H = HMRATE((1e-9)*1.3,0,0)*(TH85_field);
	/*>>chng 05 dec 17 The rate for this reaction comes from Eric Herbst's website */
	C2_nu_C_C = HMRATE((1e-11)*4.7,0,0)*(TH85_field);
	CN_NU_N_C = HMRATE(1.0e-10,0,0)*(TH85_field);
	HNC_PHOTON_CN_H = HMRATE(1.5e-9,0,0)*(TH85_field);
	HNC_CRPHOT_CN_H = crnu_rate*2*2986;
	/*>>chng 05 dec 17 The rate for this reaction comes from Eric Herbst's website */
	C2_nu_C2P_e = HMRATE((1e-10),0,0)*(TH85_field);
	C2_crnu_C_C = crnu_rate * 2. *500;
	C2P_nu_CP_C = HMRATE(1.0e-11,0,0)*(TH85_field);
	crnu_CH_C_H             = crnu_rate * 2. * 756;
	crnu_CHP_CP_H           = crnu_rate * 2. * 183;
	crnu_H2O_OH_H           = crnu_rate * 2. * 979;
	crnu_O2_O2P_e           = crnu_rate * 2. *88;
	crnu_O2_O_O             = crnu_rate * 2. *730;
	crnu_OH_O_H             = crnu_rate * 2. *522;
	crnu_SiH_Si_H           = crnu_rate * 2. *500;
	crnu_SiO_Si_O           = crnu_rate * 2. *500;
	crnu_CH3_CH3P_e         = crnu_rate * 2.*500;
	crnu_CH3_CH2_H          = crnu_rate * 2.*500;
	crnu_CH3_CH_H2          = crnu_rate * 2.*500;
	crnu_CH4_CH2_H2         = crnu_rate * 2.*2272;
	C2H_CRPHOT_C2_H = crnu_rate*2*5000;
	C2H_PHOTON_C2_H = 0.00000000051*TH85_field;
	C2H_CRPHOT_C2HP_e = crnu_rate*2*5000;
	C2H_PHOTON_C2HP_e = 0.00000000001*TH85_field;
	C3_CRPHOT_C2_C = crnu_rate*2*1119;
	C3_PHOTON_C2_C = 0.0000000038*TH85_field;
	C2H2_CRPHOT_C2H2P_e = crnu_rate*2*1309;
	C2H2_CRPHOT_C2H_H = crnu_rate*2*5155;
	C3H_CRPHOT_C3_H = crnu_rate*2*5000;
	C2H2_PHOTON_C2H2P_e = 0.00000000046*TH85_field;
	C2H2_PHOTON_C2H_H = 0.0000000073*TH85_field;
	C3H_PHOTON_C3_H = 0.000000001*TH85_field;

	/*  >>refer   mol  NH     +  NU    =>   NH+   +     e     UMIST  */

	mole.c[ipNH][ipNH] -= NH_NU_NHP_e;
	mole.c[ipNH][ipNHP] += NH_NU_NHP_e;

	/*  >>refer   mol  NH     +  NU    =>   N     +     H     UMIST  */

	mole.c[ipNH][ipNH] -= NH_NU_N_H;
	mole.c[ipNH][ipATN] += NH_NU_N_H;

	/*  >>refer   mol  NH2    +  NU    =>   NH    +     H     UMIST  */

	mole.c[ipNH2][ipNH2] -= NH2_NU_NH_H;
	mole.c[ipNH2][ipNH] += NH2_NU_NH_H;

	/*  >>refer   mol  NH2    +  NU    =>   NH2+  +     e     UMIST  */

	mole.c[ipNH2][ipNH2] -= NH2_NU_NH2P_e;
	mole.c[ipNH2][ipNH2P] += NH2_NU_NH2P_e;

	/*  >>refer   mol  NH3    +  NU    =>   NH3+  +     e     UMIST  */

	mole.c[ipNH3][ipNH3] -= NH3_NU_NH3P_e;
	mole.c[ipNH3][ipNH3P] += NH3_NU_NH3P_e;

	/*  >>refer   mol  NH3    +  NU    =>   NH    +     H2    UMIST  */

	mole.c[ipNH3][ipNH3] -= NH3_NU_NH_H2;
	mole.c[ipNH3][ipNH] += NH3_NU_NH_H2;

	/*  >>refer   mol  NH3    +  NU    =>   NH2   +     H     UMIST  */

	mole.c[ipNH3][ipNH3] -= NH3_NU_NH2_H;
	mole.c[ipNH3][ipNH2] += NH3_NU_NH2_H;

	/*  >>refer   mol  CN     +  NU    =>   N     +     C     UMIST  */

	mole.c[ipCN][ipCN] -= CN_NU_N_C;
	mole.c[ipCN][ipATN] += CN_NU_N_C;
	mole.c[ipCN][ipATC] += CN_NU_N_C;

	/*  >>refer   mol  HCN    +  NU    =>   CN    +     H     UMIST  */

	mole.c[ipHCN][ipHCN] -= HCN_NU_CN_H;
	mole.c[ipHCN][ipCN] += HCN_NU_CN_H;

	/*  >>refer   mol  N2     +  NU    =>   N     +     N     UMIST  */

	mole.c[ipN2][ipN2] -= N2_NU_N_N;
	mole.c[ipN2][ipATN] += 2*N2_NU_N_N;

	/*  >>refer   mol  NO     +  NU    =>   O     +     N     UMIST  */

	mole.c[ipNO][ipNO] -= NO_NU_O_N;
	mole.c[ipNO][ipATN] += NO_NU_O_N;
	mole.c[ipNO][ipATO] += NO_NU_O_N;

	/*  >>refer   mol  NO     +  NU    =>   NO+   +     e     UMIST  */

	mole.c[ipNO][ipNO] -= NO_NU_NOP_e;
	mole.c[ipNO][ipNOP] += NO_NU_NOP_e;

	/*  >>refer   mol  HNO    +  NU    =>   NO    +     H     UMIST  */

	mole.c[ipHNO][ipHNO] -= HNO_NU_NO_H;
	mole.c[ipHNO][ipNO] += HNO_NU_NO_H;

	/*  >>refer   mol  HS     +  NU    =>   S     +     H     UMIST  */

	mole.c[ipHS][ipHS] -= HS_NU_S_H;
	mole.c[ipHS][ipATS] += HS_NU_S_H;

	/*  >>refer   mol  HS+    +  NU    =>   S+    +     H     UMIST  */

	mole.c[ipHSP][ipHSP] -= HSP_NU_SP_H;
	mole.c[ipHSP][ipSP] += HSP_NU_SP_H;

	/*  >>refer   mol  OCN    +  NU    =>   CN    +     O     UMIST  */

	mole.c[ipOCN][ipOCN] -= OCN_NU_CN_O;
	mole.c[ipOCN][ipCN] += OCN_NU_CN_O;
	mole.c[ipOCN][ipATO] += OCN_NU_CN_O;

	/*  >>refer   mol  CS     +  NU    =>   S     +     C     UMIST  */

	mole.c[ipCS][ipCS] -= CS_NU_S_C;
	mole.c[ipCS][ipATS] += CS_NU_S_C;
	mole.c[ipCS][ipATO] += CS_NU_S_C;

	/*  >>refer   mol  CS+    +  NU    =>   S+    +     C     UMIST  */

	mole.c[ipCSP][ipCSP] -= CSP_NU_SP_C;
	mole.c[ipCSP][ipSP] += CSP_NU_SP_C;
	mole.c[ipCSP][ipATC] += CSP_NU_SP_C;

	/*  >>refer   mol  CS     +  NU    =>   CS+   +     e     UMIST  */

	mole.c[ipCS][ipCS] -= CS_NU_CSP_e;
	mole.c[ipCS][ipCSP] += CS_NU_CSP_e;

	/*  >>refer   mol  NO2    +  NU    =>   NO    +     O     UMIST  */

	mole.c[ipNO2][ipNO2] -= NO2_NU_NO_O;
	mole.c[ipNO2][ipNO] += NO2_NU_NO_O;
	mole.c[ipNO2][ipATO] += NO2_NU_NO_O;

	/*  >>refer   mol  NS     +  NU    =>   S     +     N     UMIST  */

	mole.c[ipNS][ipNS] -= NS_NU_S_N;
	mole.c[ipNS][ipATS] += NS_NU_S_N;
	mole.c[ipNS][ipATN] += NS_NU_S_N;

	/*  >>refer   mol  SO     +  NU    =>   S     +     O     UMIST  */

	mole.c[ipSO][ipSO] -= SO_NU_S_O;
	mole.c[ipSO][ipATO] += SO_NU_S_O;
	mole.c[ipSO][ipATS] += SO_NU_S_O;

	/*  >>refer   mol  SO     +  NU    =>   SO+   +     e     UMIST  */

	mole.c[ipSO][ipSO] -= SO_NU_SOP_e;
	mole.c[ipSO][ipSOP] += SO_NU_SOP_e;

	/*  >>refer   mol  OCS    +  NU    =>   OCS+  +     e     UMIST  */

	mole.c[ipOCS][ipOCS] -= OCS_NU_OCSP_e;
	mole.c[ipOCS][ipOCSP] += OCS_NU_OCSP_e;

	/*  >>refer   mol  NH     +  CRP   =>   NH+   +     e     UMIST  */

	mole.c[ipNH][ipNH] -= NH_CRP_NHP_e;
	mole.c[ipNH][ipNHP] += NH_CRP_NHP_e;

	/*  >>refer   mol  NH     +  CRP   =>   N     +     H     UMIST  */

	mole.c[ipNH][ipNH] -= NH_CRP_N_H;
	mole.c[ipNH][ipATN] += NH_CRP_N_H;

	/*  >>refer   mol  NH2    +  CRP   =>   NH    +     H     UMIST  */

	mole.c[ipNH2][ipNH2] -= NH2_CRP_NH_H;
	mole.c[ipNH2][ipNH] += NH2_CRP_NH_H;

	/*  >>refer   mol  NH2    +  CRP   =>   NH2+  +     e     UMIST  */

	mole.c[ipNH2][ipNH2] -= NH2_CRP_NH2P_e;
	mole.c[ipNH2][ipNH2P] += NH2_CRP_NH2P_e;

	/*  >>refer   mol  NH3    +  CRP   =>   NH3+  +     e     UMIST  */

	mole.c[ipNH3][ipNH3] -= NH3_CRP_NH3P_e;
	mole.c[ipNH3][ipNH3P] += NH3_CRP_NH3P_e;

	/*  >>refer   mol  NH3    +  CRP   =>   NH    +     H2    UMIST  */

	mole.c[ipNH3][ipNH3] -= NH3_CRP_NH_H2;
	mole.c[ipNH3][ipNH] += NH3_CRP_NH_H2;

	/*  >>refer   mol  NH3    +  CRP   =>   NH2   +     H     UMIST  */

	mole.c[ipNH3][ipNH3] -= NH3_CRP_NH2_H;
	mole.c[ipNH3][ipNH2] += NH3_CRP_NH2_H;

	/*  >>refer   mol  CN     +  CRP   =>   N     +     C     UMIST  */

	mole.c[ipCN][ipCN] -= CN_CRP_N_C;
	mole.c[ipCN][ipATN] += CN_CRP_N_C;
	mole.c[ipCN][ipATC] += CN_CRP_N_C;

	/*  >>refer   mol  HCN    +  CRP   =>   CN    +     H     UMIST  */

	mole.c[ipHCN][ipHCN] -= HCN_CRP_CN_H;
	mole.c[ipHCN][ipCN] += HCN_CRP_CN_H;

	/*  >>refer   mol  N2     +  CRP   =>   N     +     N     UMIST  */

	mole.c[ipN2][ipN2] -= N2_CRP_N_N;
	mole.c[ipN2][ipATN] += 2*N2_CRP_N_N;

	/*  >>refer   mol  NO     +  CRP   =>   O     +     N     UMIST  */
		
	mole.c[ipNO][ipNO] -= NO_CRP_O_N;
	mole.c[ipNO][ipATO] += NO_CRP_O_N;
	mole.c[ipNO][ipATN] += NO_CRP_O_N;

	/*  >>refer   mol  NO     +  CRP   =>   NO+   +     e     UMIST  */

	mole.c[ipNO][ipNO] -= NO_CRP_NOP_e;
	mole.c[ipNO][ipNOP] += NO_CRP_NOP_e;

	/*  >>refer   mol  HNO    +  CRP   =>   NO    +     H     UMIST  */

	mole.c[ipHNO][ipHNO] -= HNO_CRP_NO_H;
	mole.c[ipHNO][ipNO] += HNO_CRP_NO_H;

	/*  >>refer   mol  HS     +  CRP   =>   S     +     H     UMIST  */

	mole.c[ipHS][ipHS] -= HS_CRP_S_H;
	mole.c[ipHS][ipATS] += HS_CRP_S_H;

	/*  >>refer   mol  OCN    +  CRP   =>   CN    +     O     UMIST  */

	mole.c[ipOCN][ipOCN] -= OCN_CRP_CN_O;
	mole.c[ipOCN][ipCN] += OCN_CRP_CN_O;
	mole.c[ipOCN][ipATO] += OCN_CRP_CN_O;

	/*  >>refer   mol  CS     +  CRP   =>   S     +     C     UMIST  */

	mole.c[ipCS][ipCS] -= CS_CRP_S_C;
	mole.c[ipCS][ipATS] += CS_CRP_S_C;
	mole.c[ipCS][ipATC] += CS_CRP_S_C;

	/*  >>refer   mol  CS     +  CRP   =>   CS+   +     e     UMIST  */

	mole.c[ipCS][ipCS] -= CS_CRP_CSP_e;
	mole.c[ipCS][ipCSP] += CS_CRP_CSP_e;

	/*  >>refer   mol  NO2    +  CRP   =>   NO    +     O     UMIST  */

	mole.c[ipNO2][ipNO2] -= NO2_CRP_NO_O;
	mole.c[ipNO2][ipNO] += NO2_CRP_NO_O;
	mole.c[ipNO2][ipATO] += NO2_CRP_NO_O;

	/*  >>refer   mol  NS     +  CRP   =>   S     +     N     UMIST  */

	mole.c[ipNS][ipNS] -= NS_CRP_S_N;
	mole.c[ipNS][ipATS] += NS_CRP_S_N;
	mole.c[ipNS][ipATN] += NS_CRP_S_N;

	/*  >>refer   mol  SO     +  CRP   =>   S     +     O     UMIST  */

	mole.c[ipSO][ipSO] -= SO_CRP_S_O;
	mole.c[ipSO][ipATS] += SO_CRP_S_O;
	mole.c[ipSO][ipATO] += SO_CRP_S_O;

	/*  >>refer   mol  SO     +  CRP   =>   SO+   +     e     UMIST  */

	mole.c[ipSO][ipSO] -= SO_CRP_SOP_e;
	mole.c[ipSO][ipSOP] += SO_CRP_SOP_e;

	/*  >>refer   mol  OCS    +  CRP   =>   OCS+  +     e     UMIST  */

	mole.c[ipOCS][ipOCS] -= OCS_CRP_OCSP_e;
	mole.c[ipOCS][ipOCSP] += OCS_CRP_OCSP_e;

	/* >>refer	mol	PHOTON	CH	=>	CH+	e	UMIST	*/
		
	mole.c[ipCH][ipCH] -= nu_CH_CHP_e;
	mole.c[ipCH][ipCHP] += nu_CH_CHP_e;
	
	/* >>refer	mol	PHOTON	CH	=>	C	H	UMIST	*/
		
	mole.c[ipCH][ipCH] -= nu_CH_C_H;
	mole.c[ipCH][ipATC] += nu_CH_C_H;
	
	/* >>refer	mol	PHOTON	CH+	=>	C+	H	UMIST	*/
		
	mole.c[ipCHP][ipCHP] -= nu_CHP_CP_H;
	mole.c[ipCHP][ipCP]  += nu_CHP_CP_H;
	
	/* >>refer	mol	PHOTON	CH2	=>	CH2+	e	UMIST	*/
		
	mole.c[ipCH2][ipCH2] -= nu_CH2_CH2P_e;
	mole.c[ipCH2][ipCH2P] += nu_CH2_CH2P_e;
	
	/* >>refer	mol	PHOTON	CH2	=>	CH	H	UMIST	*/
		
	mole.c[ipCH2][ipCH2] -= nu_CH2_CH_H;
	mole.c[ipCH2][ipCH] += nu_CH2_CH_H;
	
	/* >>refer	mol	PHOTON	CH2+	=>	CH+	H	UMIST	*/
		
	mole.c[ipCH2P][ipCH2P] -= nu_CH2P_CHP_H;
	mole.c[ipCH2P][ipCHP] += nu_CH2P_CHP_H;
	
	/* >>refer	mol	PHOTON	CH3+	=>	CH2+	H	UMIST	*/
		
	mole.c[ipCH3P][ipCH3P] -= nu_CH3P_CH2P_H;
	mole.c[ipCH3P][ipCH2P] += nu_CH3P_CH2P_H;
	
	/* >>refer	mol	PHOTON	CH3+	=>	CH+	H2	UMIST	*/
		
	mole.c[ipCH3P][ipCH3P] -= nu_CH3P_CHP_H2;
	mole.c[ipCH3P][ipCHP] += nu_CH3P_CHP_H2;
	/* >>refer	mol	PHOTON	H2O	=>	H2O+	e	UMIST	*/
		
	mole.c[ipH2O][ipH2O] -= nu_H2O_H2OP_e;
	mole.c[ipH2O][ipH2OP] += nu_H2O_H2OP_e;
	
	/* >>refer	mol	PHOTON	H2O	=>	OH	H	UMIST	*/
		
	mole.c[ipH2O][ipH2O] -= nu_H2O_OH_H;
	mole.c[ipH2O][ipOH] += nu_H2O_OH_H;
	
	/* >>refer	mol	PHOTON	O2	=>	O2+	e	UMIST	*/
		
	mole.c[ipO2][ipO2] -= nu_O2_O2P_e;
	mole.c[ipO2][ipO2P] += nu_O2_O2P_e;
	
	/* >>refer	mol	PHOTON	O2	=>	O	O	UMIST	*/
	
	mole.c[ipO2][ipO2] -= nu_O2_O_O;
	mole.c[ipO2][ipATO] += 2*nu_O2_O_O;
	
	/* >>refer	mol	PHOTON	OH	=>	OH+	e	UMIST	*/
		
	mole.c[ipOH][ipOH] -= nu_OH_OHP_e;
	mole.c[ipOH][ipOHP] += nu_OH_OHP_e;
	
	/* >>refer	mol	PHOTON	OH	=>	O	H	UMIST	*/
		

	/* A note here is that approx. 55% of the O produced in this reaction goes into the 1D
	   state, which is the upper level of the 6300 Angstrom line, while 5% forms in the 1S
	   state and can create the 5577 line */
	
	/* >>chng 03 oct 20, following had been deleted so OH formation was always zero,
	 * this is used to get OH photo formation of [OI] 6300
	 *>>refer	OI	photoexcitation	Storzer, H., & Hollenbach, D. 2000, ApJ, 539, 751-759 
	 * discussion on bottom left side of page 752 */

	co.rate_OH_dissoc = (float)(nu_OH_O_H* OH);

	mole.c[ipOH][ipOH] -= nu_OH_O_H;
	mole.c[ipOH][ipATO] += nu_OH_O_H;
	
	/* >>refer	mol	PHOTON	OH+	=>	O	H+	UMIST	*/
		
	mole.c[ipOHP][ipOHP] -= nu_OHP_O_HP;
	mole.c[ipOHP][ipATO] += nu_OHP_O_HP;
	
	/* >>refer	mol	PHOTON	SiH	=>	Si	H	UMIST	*/
		
	mole.c[ipSIH][ipSIH] -= nu_SiH_Si_H;
	mole.c[ipSIH][ipATSI] += nu_SiH_Si_H;
	
	/* >>refer	mol	PHOTON	SiO	=>	Si	O	UMIST	*/
		
	mole.c[ipSIO][ipSIO] -= nu_SiO_Si_O;
	mole.c[ipSIO][ipATO] += nu_SiO_Si_O;
	mole.c[ipSIO][ipATSI] += nu_SiO_Si_O;
	

	/*  >>refer  mol  HNC   +   PHOTON   =   CN   +   H   */

	mole.c[ipHNC][ipHNC] -= HNC_PHOTON_CN_H; 
	mole.c[ipHNC][ipCN] += HNC_PHOTON_CN_H; 

	/*  >>refer  mol  HNC   +   CRPHOT   =   CN   +   H   */

	mole.c[ipHNC][ipHNC] -= HNC_CRPHOT_CN_H; 
	mole.c[ipHNC][ipCN] += HNC_CRPHOT_CN_H; 

	/*  >>refer   mol  HCl    +  NU     =>   Cl    +     H   UMIST  */

	mole.c[ipHCl][ipHCl] -= HCl_NU_Cl_H;
	mole.c[ipHCl][ipATCl] += HCl_NU_Cl_H;

	/*  >>refer   mol  CCl    +  NU     =>   Cl    +     C    UMIST  */

	mole.c[ipCCl][ipCCl] -= CCl_NU_Cl_C;
	mole.c[ipCCl][ipATCl] += CCl_NU_Cl_C;
	mole.c[ipCCl][ipATC] += CCl_NU_Cl_C;

	/*  >>refer   mol  ClO    +  NU     =>   Cl    +     O    UMIST  */

	mole.c[ipClO][ipClO] -= ClO_NU_Cl_O;
	mole.c[ipClO][ipATCl] += ClO_NU_Cl_O;
	mole.c[ipClO][ipATO] += ClO_NU_Cl_O;

	/*  >>refer   mol  HCl    +  crnu    =>   Cl    +     H    UMIST  */

	mole.c[ipHCl][ipHCl] -= HCl_crnu_Cl_H;
	mole.c[ipHCl][ipATCl] += HCl_crnu_Cl_H;

	/*  >>refer   mol  CCl    +  crnu    =>   Cl    +     C    UMIST  */

	mole.c[ipCCl][ipCCl] -= CCl_crnu_Cl_C;
	mole.c[ipCCl][ipATCl] += CCl_crnu_Cl_C;
	mole.c[ipCCl][ipATC] += CCl_crnu_Cl_C;

	/*  >>refer   mol  ClO    +  crnu    =>   Cl    +     O    UMIST  */

	mole.c[ipClO][ipClO] -= ClO_crnu_Cl_O;
	mole.c[ipClO][ipATCl] += ClO_crnu_Cl_O;
	mole.c[ipClO][ipATO] += ClO_crnu_Cl_O;

		/* >>refer	mol	cr	CH	=>	C	H	UMIST	*/
		
	mole.c[ipCH][ipCH] -= crnu_CH_C_H;
	mole.c[ipCH][ipATC] += crnu_CH_C_H;
	
	/* >>refer	mol	cr	CH+	=>	C+	H	UMIST	*/
		
	mole.c[ipCHP][ipCHP] -= crnu_CHP_CP_H;
	mole.c[ipCHP][ipCP]  += crnu_CHP_CP_H;
	
	/* >>refer	mol	cr	H2O	=>	OH	H	UMIST	*/
		
	mole.c[ipH2O][ipH2O] -= crnu_H2O_OH_H;
	mole.c[ipH2O][ipOH] += crnu_H2O_OH_H;
	
	/* >>refer	mol	cr	O2	=>	O2+	e	UMIST	*/
		
	mole.c[ipO2][ipO2] -= crnu_O2_O2P_e;
	mole.c[ipO2][ipO2P] += crnu_O2_O2P_e;
	
	/* >>refer	mol	cr	O2	=>	O	O	UMIST	*/
		
	mole.c[ipO2][ipO2] -= crnu_O2_O_O;
	mole.c[ipO2][ipATO] += 2*crnu_O2_O_O;
	
	/* >>refer	mol	cr	OH	=>	O	H	UMIST	*/
		
	mole.c[ipOH][ipOH] -= crnu_OH_O_H;
	mole.c[ipOH][ipATO] += crnu_OH_O_H;
	
	/* >>refer	mol	cr	SiH	=>	Si	H	UMIST	*/
		
	mole.c[ipSIH][ipSIH] -= crnu_SiH_Si_H;
	mole.c[ipSIH][ipATSI] += crnu_SiH_Si_H;
	
	/* >>refer	mol	cr	SiO	=>	Si	O	UMIST	*/
		
	mole.c[ipSIO][ipSIO] -= crnu_SiO_Si_O;
	mole.c[ipSIO][ipATSI] += crnu_SiO_Si_O;
	mole.c[ipSIO][ipATO] += crnu_SiO_Si_O;
	
	/* >>refer mol   nu CH3  =>  CH3P e   UMIST    */
	mole.c[ipCH3][ipCH3]  -=  nu_CH3_CH3P_e;
	mole.c[ipCH3][ipCH3P]  +=  nu_CH3_CH3P_e;
	
	/* >>refer mol   nu CH3  =>  CH2 H   UMIST    */
	mole.c[ipCH3][ipCH3]  -=  nu_CH3_CH2_H;
	mole.c[ipCH3][ipCH2]  +=  nu_CH3_CH2_H;
	
	/* >>refer mol   nu CH3  =>  CH H2   UMIST    */
	mole.c[ipCH3][ipCH3]  -=  nu_CH3_CH_H2;
	mole.c[ipCH3][ipCH]  +=  nu_CH3_CH_H2;
	
	/* >>refer mol   crnu CH3  =>  CH3P e   UMIST    */
	mole.c[ipCH3][ipCH3]  -=  crnu_CH3_CH3P_e;
	mole.c[ipCH3][ipCH3P]  +=  crnu_CH3_CH3P_e;
	
	/* >>refer mol   crnu CH3  =>  CH2 H   UMIST    */
	mole.c[ipCH3][ipCH3]  -=  crnu_CH3_CH2_H;
	mole.c[ipCH3][ipCH2]  +=  crnu_CH3_CH2_H;
	
	/* >>refer mol   crnu CH3  =>  CH H2   UMIST    */
	mole.c[ipCH3][ipCH3]  -=  crnu_CH3_CH_H2;
	mole.c[ipCH3][ipCH]  +=  crnu_CH3_CH_H2;

	/* >>refer	mol	C + crnu => C+	UMIST */
	/* >>chng 04 sep 19, rm this reaction as per NA email */
	/* This reaction should only be used when the UMIST rate file is wanted */ 
	/* >>chng 05 oct 01, following two turned on as per NA email 
	 This should be uncommented for the comparison.  I know the comment says that a previous e-mail of mine says to comment this out.  I had originally said that this should be commented out because the cosmic rays were included in the code.  However, the enhancement factor of 2*510 is not included.
	 Therefore, for the comparison, this should be in.  This makes a minor difference in two of the models we did not compare in the group email.
	 However, since a new model might be chosen, we should put this in.  The code should now be:
	 Where the extra term is added because the commented code predates the inclusion of sources and sinks in the heavy element chemistry.  I have no current idea of what line number this is, but if you grep on "as per NA" in mole_co_step.c, then you will find it.  If we merge the code I have changed tomorrow, then this change will also be included.  I will run the hacks with the changes and submit the data.
	 */
	if(!co.lgUMISTrates)
	{
		mole.c[ipATC][ipCP] += 5e-17* 2.*510 ;
		mole.c[ipATC][ipATC] -= 5e-17* 2.*510 ;
	}

	/* >>refer mol   nu CH4  =>  CH3 H   UMIST    */
	mole.c[ipCH4][ipCH4]  -=  nu_CH4_CH3_H;
	mole.c[ipCH4][ipCH3]  +=  nu_CH4_CH3_H;	

	/* >>refer mol   nu CH4  =>  CH2 H2   UMIST    */
	mole.c[ipCH4][ipCH4]  -=  nu_CH4_CH2_H2;
	mole.c[ipCH4][ipCH2]  +=  nu_CH4_CH2_H2;
	
	/* >>refer mol   nu CH4  =>  CH H2 H  UMIST    */
	mole.c[ipCH4][ipCH4]  -=  nu_CH4_CH_H2;
	mole.c[ipCH4][ipCH]  +=  nu_CH4_CH_H2;
	
	/* >>refer mol   crnu CH4  =>  CH2 H2   UMIST    */
	mole.c[ipCH4][ipCH4]  -=  crnu_CH4_CH2_H2;
	mole.c[ipCH4][ipCH2]  +=  crnu_CH4_CH2_H2;
	
	/* >>refer	mol	C2		nu	=>		C		C		UMIST	*/
	
	mole.c[ipC2][ipC2] -= C2_nu_C_C;
	mole.c[ipC2][ipATC]  += C2_nu_C_C;
	mole.c[ipC2][ipATC]  += C2_nu_C_C;

	/* >>refer	mol	C2		nu	=>		C2+		e		UMIST	*/
	
	mole.c[ipC2][ipC2] -= C2_nu_C2P_e;
	mole.c[ipC2][ipC2P]  += C2_nu_C2P_e;

	/* >>refer	mol	C2		cr	=>		C		C		UMIST	*/
	
	mole.c[ipC2][ipC2] -= C2_crnu_C_C;
	mole.c[ipC2][ipATC]  += C2_crnu_C_C;
	mole.c[ipC2][ipATC]  += C2_crnu_C_C;

	/* >>refer	mol	C2+		nu	=>		C+		C		UMIST	*/
	
	mole.c[ipC2P][ipC2P] -= C2P_nu_CP_C;
	mole.c[ipC2P][ipCP]  += C2P_nu_CP_C;
	mole.c[ipC2P][ipATC]  += C2P_nu_CP_C;
	

	/* >>refer	mol	C + hnu => C+	e	 Done internally by Cloudy*/

	/* Photoionization and grain ct ionization of Co */
	mole.c[ipATC][ipCP] += (ionbal.RateIonizTot[ipCARBON][0]+gv.GrainChTrRate[ipCARBON][0][1]);
	mole.c[ipATC][ipATC] -= (ionbal.RateIonizTot[ipCARBON][0]+gv.GrainChTrRate[ipCARBON][0][1]);

	/* >>refer	mol	O + hnu => O+	e	Done internally by Cloudy*/

	/* Photoionization and grain ct ionization of Oo */
	mole.c[ipATO][ipOP] += (ionbal.RateIonizTot[ipOXYGEN][0]+gv.GrainChTrRate[ipOXYGEN][0][1]);
	mole.c[ipATO][ipATO] -= (ionbal.RateIonizTot[ipOXYGEN][0]+gv.GrainChTrRate[ipOXYGEN][0][1]);

	/* >>refer	mol	Si + hnu => Si+	e	Done internally by Cloudy*/

	/* Photoionization of Si */

	mole.c[ipATSI][ipSIP] += (ionbal.RateIonizTot[ipSILICON][0]+gv.GrainChTrRate[ipSILICON][0][1]);
	mole.c[ipATSI][ipATSI] -= (ionbal.RateIonizTot[ipSILICON][0]+gv.GrainChTrRate[ipSILICON][0][1]);

	/* >>refer	mol	Cl + hnu => Cl+	e	 Done internally by Cloudy*/

	/* Photoionization of S */
	mole.c[ipATCl][ipClP] += (ionbal.RateIonizTot[ipCHLORINE][0]+gv.GrainChTrRate[ipCHLORINE][0][1]);
	mole.c[ipATCl][ipATCl] -= (ionbal.RateIonizTot[ipCHLORINE][0]+gv.GrainChTrRate[ipCHLORINE][0][1]);


	/* >>refer	mol	N + hnu => N+	e	 Done internally by Cloudy*/

	/* Photoionization of N */
	mole.c[ipATN][ipNP] += (ionbal.RateIonizTot[ipNITROGEN][0]+gv.GrainChTrRate[ipNITROGEN][0][1]);
	mole.c[ipATN][ipATN] -= (ionbal.RateIonizTot[ipNITROGEN][0]+gv.GrainChTrRate[ipNITROGEN][0][1]);


	/* >>refer	mol	S + hnu => S+	e	 Done internally by Cloudy*/

	/* Photoionization of S */
	mole.c[ipATS][ipSP] += (ionbal.RateIonizTot[ipSULPHUR][0]+gv.GrainChTrRate[ipSULPHUR][0][1]);
	mole.c[ipATS][ipATS] -= (ionbal.RateIonizTot[ipSULPHUR][0]+gv.GrainChTrRate[ipSULPHUR][0][1]);

	/* >>refer  mole   C2H   CRPHOT   -->   C2   H        */				
	
	mole.c[ipC2H][ipC2H] -= C2H_CRPHOT_C2_H;				
	mole.c[ipC2H][ipC2] += C2H_CRPHOT_C2_H;				

	/* >>refer  mole   C2H   PHOTON   -->   C2   H        */

	mole.c[ipC2H][ipC2H] -= C2H_PHOTON_C2_H;	
	mole.c[ipC2H][ipC2] += C2H_PHOTON_C2_H;	

	/* >>refer  mole   C2H   PHOTON   -->   C2H+   e-        */	

	mole.c[ipC2H][ipC2H] -= C2H_PHOTON_C2HP_e;		
	mole.c[ipC2H][ipC2HP] += C2H_PHOTON_C2HP_e;		

	/* >>refer  mole   C2H   CRPHOT   -->   C2H+   e-        */					
	
	mole.c[ipC2H][ipC2H] -= C2H_CRPHOT_C2HP_e;					
	mole.c[ipC2H][ipC2HP] += C2H_CRPHOT_C2HP_e;					

	/* >>refer  mole   C3   PHOTON   -->   C2   C        */

	mole.c[ipC3][ipC3] -= C3_PHOTON_C2_C;	
	mole.c[ipC3][ipC2] += C3_PHOTON_C2_C;	
	mole.c[ipC3][ipATC] += C3_PHOTON_C2_C;	

	/* >>refer  mole   C3   CRPHOT   -->   C2   C        */	
	
	mole.c[ipC3][ipC3] -= C3_CRPHOT_C2_C;	
	mole.c[ipC3][ipATC] += C3_CRPHOT_C2_C;	
	mole.c[ipC3][ipC2] += C3_CRPHOT_C2_C;	

	/* >>refer  mole   C2H2   CRPHOT   -->   C2H2+   e-        */			
	
	mole.c[ipC2H2][ipC2H2] -= C2H2_CRPHOT_C2H2P_e;			
	mole.c[ipC2H2][ipC2H2P] += C2H2_CRPHOT_C2H2P_e;			

	/* >>refer  mole   C2H2   CRPHOT   -->   C2H   H        */		
	
	mole.c[ipC2H2][ipC2H2] -= C2H2_CRPHOT_C2H_H;		
	mole.c[ipC2H2][ipC2H] += C2H2_CRPHOT_C2H_H;		

	/* >>refer  mole   C3H   CRPHOT   -->   C3   H        */
	
	mole.c[ipC3H][ipC3H] -= C3H_CRPHOT_C3_H;					
	mole.c[ipC3H][ipC3] += C3H_CRPHOT_C3_H;					

	/* >>refer  mole   C2H2   PHOTON   -->   C2H2+   e-        */			

	mole.c[ipC2H2][ipC2H2] -= C2H2_PHOTON_C2H2P_e;				
	mole.c[ipC2H2][ipC2H2P] += C2H2_PHOTON_C2H2P_e;				

	/* >>refer  mole   C2H2   PHOTON   -->   C2H   H        */		

	mole.c[ipC2H2][ipC2H2] -= C2H2_PHOTON_C2H_H;			
	mole.c[ipC2H2][ipC2H] += C2H2_PHOTON_C2H_H;			

	/* >>refer  mole   C3H   PHOTON   -->   C3   H        */				

	mole.c[ipC3H][ipC3H] -= C3H_PHOTON_C3_H;
	mole.c[ipC3H][ipC3] += C3H_PHOTON_C3_H;


}
