00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "cddefines.h"
00020 #include "physconst.h"
00021 #include "iso.h"
00022 #include "atmdat.h"
00023 #include "grainvar.h"
00024 #include "ionbal.h"
00025 #include "dense.h"
00026 #include "secondaries.h"
00027 #include "mole.h"
00028 #include "mole_co_priv.h"
00029 #include "opacity.h"
00030 #include "rfield.h"
00031 #include "thermal.h"
00032 #include "timesc.h"
00033 #include "trace.h"
00034 #include "phycon.h"
00035 #include "doppvel.h"
00036 #include "thirdparty.h"
00037 #include "gammas.h"
00038 #include "h2.h"
00039 #include "dynamics.h"
00040 #include "conv.h"
00041 #include "radius.h"
00042 #include "hextra.h"
00043 #include "hmi.h"
00044
00045
00046 #if defined(__HP_aCC)
00047 # pragma OPT_LEVEL 1
00048 #endif
00049
00050 #define ABSLIM 1e-12
00051
00052
00053 #define INTSZ(a) (sizeof(a)/sizeof(int))
00054
00055 struct Hmole_rate_s {
00056 int index;
00057 int nreactants, nrates, nproducts;
00058 int reactants[MAXREACTANTS];
00059 int rate_species[MAXREACTANTS];
00060 int products[MAXPRODUCTS];
00061 double rk;
00062 struct Hmole_rate_s *next;
00063 };
00064 typedef struct Hmole_rate_s reaction;
00065
00066
00067 reaction *newreaction(
00068
00069 int rindex,
00070
00071 int *in,
00072
00073 int nin,
00074
00075 int *out,
00076
00077 int nout,
00078
00079
00080
00081
00082 int *rate,
00083
00084 int nrate)
00085 {
00086 static reaction *list = NULL, *r;
00087 static int poolsize=1, index = 0;
00088 int i;
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 if(rate == NULL)
00106 {
00107 rate=in;
00108 nrate=nin;
00109 }
00110
00111
00112
00113 if(list == NULL || index == poolsize)
00114 {
00115 poolsize <<=1;
00116 list = ((reaction *)MALLOC( (size_t)poolsize*sizeof(reaction) ));
00117 index = 0;
00118 }
00119
00120
00121
00122 r = list+index;
00123 index++;
00124 r->next = NULL;
00125 r->index = rindex;
00126 ASSERT(nin <= MAXREACTANTS && nout <= MAXPRODUCTS && nrate <= MAXREACTANTS);
00127
00128 r->nreactants = nin;
00129 r->nrates = nrate;
00130 r->nproducts = nout;
00131
00132
00133 for(i=0; i<r->nreactants; i++)
00134 r->reactants[i] = in[i];
00135
00136
00137 for(i=0; i<r->nrates; i++)
00138 r->rate_species[i] = rate[i];
00139
00140
00141 for(i=0; i<r->nproducts; i++)
00142 r->products[i] = out[i];
00143
00144 return r;
00145 }
00146
00147
00148 #if defined (__ICC) && defined(__amd64)
00149 #pragma optimization_level 1
00150 #endif
00151
00152
00153 void hmole_step(int *nFixup, double *error)
00154 {
00155 enum {PRINTSOL = false};
00156
00157 int32 ipiv[N_H_MOLEC];
00158
00159 long int i,
00160 ipConserve,
00161 j,
00162 limit ,
00163 nd,
00164 mol;
00165
00166 int printsol = PRINTSOL;
00167
00168 bool lgNegPop;
00169 int iworst;
00170
00171 double frac_H2star_grains,
00172 frac_H2star_hminus;
00173 double sum_first_ions;
00174
00175 double
00176 bhneut,
00177 Hneut,
00178 c3bod,
00179 cionhm,
00180 corr,
00181 H2star_deexcit,
00182 deexc_htwo,
00183 deexc_hneut,
00184 desh2p,
00185 etmp,
00186 eh3p_3h,
00187 Boltz_fac_H2_H2star,
00188 fhneut,
00189 gamheh,
00190 h1fnd,
00191 h1rat,
00192 h2pcin,
00193 h2phhp,
00194 h2pion,
00195 H2star_excit ,
00196 radath,
00197 hmihph2p,
00198 h2phmhhh,
00199 h2crphh,
00200 h2scrphh,
00201 h2crphphm,
00202 h2scrphphm,
00203 h2crphpeh,
00204 h2scrphpeh,
00205 h2crh2pe,
00206 h2crhphe,
00207 h2scrhphe,
00208 h2scrh2pe,
00209 h2pehh,
00210 h3ph2ph,
00211 hphmhhpe,
00212 h2hmhh2e,
00213 hehmeheh,
00214 hephmhhe,
00215 fracneg,
00216 fracnegtmp,
00217 fracnegfac,
00218 sum_H0_Hp,
00219 conserve,
00220 rate,
00221 rk,
00222 rated,
00223 rate_deriv[MAXREACTANTS],
00224 sinkrate[MAXREACTANTS],
00225 T_ortho_para_crit ,
00226 TStancil;
00227
00228
00229 static double
00230 gamtwo,
00231 h2phet,
00232 proton_sum_old,
00233 proton_sum_new;
00234
00235 static double
00236 b2pcin,
00237 *amat=NULL,
00238 *bvec=NULL,
00239 *Hmolec_old=NULL,
00240 **c=NULL,
00241 plte;
00242 static double oatomic = -1., oion = -1.;
00243
00244
00245
00246 static bool lgMustMalloc = true;
00247
00248 static reaction *rlist = NULL;
00249 reaction *r;
00250 long int rindex, ratei, ratej;
00251
00252 DEBUG_ENTRY( "hmole_step()" );
00253
00254 if( lgMustMalloc )
00255 {
00256
00257 lgMustMalloc = false;
00258
00259 bvec = ((double*)MALLOC( (size_t)N_H_MOLEC*sizeof(double) ));
00260 Hmolec_old = ((double*)MALLOC( (size_t)N_H_MOLEC*sizeof(double) ));
00261 amat = ((double*)MALLOC( (size_t)(N_H_MOLEC*N_H_MOLEC)*sizeof(double) ));
00262 c = ((double**)MALLOC( (size_t)N_H_MOLEC*sizeof(double *) ));
00263 for( i=0; i<N_H_MOLEC; ++i )
00264 {
00265
00266
00267
00268 c[i] = ((double*)MALLOC( (size_t)N_H_MOLEC*sizeof(double) ));
00269 }
00270 }
00271
00272
00273 *error = 0;
00274
00275
00276
00277
00278
00279 if( hmi.lgNoH2Mole )
00280 {
00281 dense.xMolecules[ipHYDROGEN] = 0.;
00282
00283
00284 for(mol=0; mol<N_H_MOLEC; ++mol)
00285 {
00286 hmi.Hmolec[mol] = 0.;
00287 }
00288 hmi.Hmolec[ipMH] = dense.xIonDense[ipHYDROGEN][0];
00289 hmi.Hmolec[ipMHp] = dense.xIonDense[ipHYDROGEN][1];
00290 hmi.H2_total = 0.;
00291
00292 dense.xIonDense[LIMELM+2][0] = 0.;
00293 hmi.hmihet = 0.;
00294 hmi.h2plus_heat = 0.;
00295 hmi.H2Opacity = 0.;
00296 hmi.hmicol = 0.;
00297 hmi.hmidep = 1.;
00298 hmi.rh2dis = 0.;
00299 hmi.HalphaHmin = 0.;
00300
00301 hmi.HeatH2Dish_used = 0.;
00302 hmi.HeatH2Dish_BigH2 = 0.;
00303 hmi.HeatH2Dish_TH85 = 0.;
00304 hmi.HeatH2Dish_BD96 = 0.;
00305 hmi.HeatH2Dish_BHT90 = 0.;
00306 hmi.HeatH2Dish_ELWERT = 0.;
00307
00310 hmi.HeatH2Dexc_used = 0.;
00311 hmi.HeatH2Dexc_BigH2 = 0.;
00312 hmi.HeatH2Dexc_TH85 = 0.;
00313 hmi.HeatH2Dexc_BD96 = 0.;
00314 hmi.HeatH2Dexc_BHT90 = 0.;
00315 hmi.HeatH2Dexc_ELWERT = 0.;
00316
00318 hmi.deriv_HeatH2Dexc_used = 0.;
00319 hmi.deriv_HeatH2Dexc_BigH2 = 0.;
00320 hmi.deriv_HeatH2Dexc_TH85 = 0.;
00321 hmi.deriv_HeatH2Dexc_BD96 = 0.;
00322 hmi.deriv_HeatH2Dexc_BHT90 = 0.;
00323 hmi.deriv_HeatH2Dexc_ELWERT = 0.;
00324 return;
00325 }
00326
00327
00328
00329
00330 if( hmi.H2_frac_abund_set>0.)
00331 {
00332 for(mol=0;mol<N_H_MOLEC;mol++)
00333 {
00334 hmi.Hmolec[mol] = 0.;
00335 }
00336
00337
00338 dense.xIonDense[ipHYDROGEN][0] = dense.xIonDense[ipHYDROGEN][1] =
00339 2.f*SMALLFLOAT*dense.gas_phase[ipHYDROGEN];
00340
00341 hmi.Hmolec[ipMH2g] = (realnum)(dense.gas_phase[ipHYDROGEN] * hmi.H2_frac_abund_set);
00342 hmi.Hmolec[ipMH2s] = 0.;
00343
00344 hmi.H2_total = hmi.Hmolec[ipMH2g] + hmi.Hmolec[ipMH2s];
00345
00346 h2.ortho_density = 0.75*hmi.H2_total;
00347 h2.para_density = 0.25*hmi.H2_total;
00348
00349 hmi.hmihet = 0.;
00350 hmi.h2plus_heat = 0.;
00351 hmi.H2Opacity = 0.;
00352 hmi.hmicol = 0.;
00353 hmi.HeatH2Dish_TH85 = 0.;
00354 hmi.HeatH2Dexc_TH85 = 0.;
00355 hmi.deriv_HeatH2Dexc_TH85 = 0.;
00356 hmi.hmidep = 1.;
00357 hmi.HalphaHmin = 0.;
00358
00359 for( nd=0; nd < gv.nBin; nd++ )
00360 {
00361 gv.bin[nd]->rate_h2_form_grains_used = 0.;
00362 }
00363 return;
00364 }
00365
00366
00367 hmi.Hmolec[ipMH] = dense.xIonDense[ipHYDROGEN][0];
00368 hmi.Hmolec[ipMHp] = dense.xIonDense[ipHYDROGEN][1];
00369
00370 for(mol=0;mol<N_H_MOLEC;mol++)
00371 {
00372 Hmolec_old[mol] = hmi.Hmolec[mol];
00373 }
00374 for(i=0; i<MAXREACTANTS; ++i )
00375 {
00376 rate_deriv[i] = 0.;
00377 sinkrate[i] = 0.;
00378 }
00379
00380
00381 if( phycon.te < 3074. )
00382 {
00383 cionhm = 1.46e-32*(powi(phycon.te,6))*phycon.sqrte*hmi.exphmi;
00384 }
00385 else if( phycon.te >= 3074. && phycon.te < 30000. )
00386 {
00387 cionhm = 5.9e-19*phycon.tesqrd*phycon.sqrte*phycon.te05;
00388 }
00389 else
00390 {
00391 cionhm = 3e-7;
00392 }
00393
00394
00395
00396
00397 if( gv.lgDustOn )
00398 {
00399
00400 # ifndef IGNORE_QUANTUM_HEATING
00401
00402 GrainDrive();
00403 # endif
00404
00405
00406
00407
00408 hmi.rate_grain_h2_op_conserve = 0.;
00409
00410 hmi.rate_grain_h2_J1_to_J0 = 0.;
00411
00412
00413 for( nd=0; nd < gv.nBin; nd++ )
00414 {
00415 # ifndef IGNORE_QUANTUM_HEATING
00416 long k, qnbin;
00417 double *qtemp, *qprob;
00418 bool lgUseQHeat = gv.lgGrainPhysicsOn && gv.bin[nd]->lgQHeat;
00419 # endif
00420
00421
00422
00423
00424
00425
00426
00427 double sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) +
00428 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
00429
00430 # ifndef IGNORE_QUANTUM_HEATING
00431
00432 if( lgUseQHeat )
00433 {
00434 qtemp = (double*)MALLOC((size_t)(NQGRID*sizeof(double)));
00435 qprob = (double*)MALLOC((size_t)(NQGRID*sizeof(double)));
00436
00437 qheat(qtemp,qprob,&qnbin,nd);
00438
00439 if( gv.bin[nd]->lgUseQHeat )
00440 {
00441 ASSERT( qnbin > 0 );
00442 }
00443 else
00444 {
00445 qnbin = 1;
00446 qprob[0] = 1.;
00447 qtemp[0] = gv.bin[nd]->tedust;
00448 }
00449
00450 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;
00451
00452 for( k=0; k < qnbin; k++ )
00453 {
00454
00455
00456
00457
00458
00459
00460
00461
00462 double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./qtemp[k]));
00463 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) +
00464 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
00465
00466 gv.bin[nd]->rate_h2_form_grains_HM79 += qprob[k] * sticking_probability_H *
00467 conversion_efficiency_HM79;
00468 }
00469
00470
00471
00472
00473 gv.bin[nd]->rate_h2_form_grains_HM79 *= 0.5 * DoppVel.AveVel[ipHYDROGEN]*
00474 gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;
00475
00476 ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
00477 }
00478 else
00479 # endif
00480 {
00481
00482
00483
00484
00485
00486
00487
00488
00489 double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./gv.bin[nd]->tedust));
00490
00491
00492
00493
00494 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.5 * DoppVel.AveVel[ipHYDROGEN]* gv.bin[nd]->IntArea/4. *
00495
00496 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * conversion_efficiency_HM79;
00497 ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
00498 }
00499
00500 # ifndef IGNORE_QUANTUM_HEATING
00501 if( lgUseQHeat )
00502 {
00503
00504
00505
00506 double f = 1e-10;
00507
00508
00509 double sqrt_term = 35.399494936611667;
00510
00511 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;
00512
00513 for( k=0; k < qnbin; k++ )
00514 {
00515 double beta_alpha = 0.25 * sqrt_term *sexp(200./qtemp[k] );
00516
00517 double xi = 1./ (1. + 1.3e13*sexp(1.5*1e4/qtemp[k])*sqrt_term/(2.*f) );
00518
00519 double beta = 3e12 * sexp( 320. / qtemp[k] );
00520
00521
00522 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );
00523 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) +
00524 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
00525
00526
00527
00528
00529 gv.bin[nd]->rate_h2_form_grains_CT02 += qprob[k] * sticking_probability_H *
00530 recombination_efficiency_CT02;
00531 }
00532
00533
00534
00535
00536
00537 gv.bin[nd]->rate_h2_form_grains_CT02 *= 0.5 * DoppVel.AveVel[ipHYDROGEN]*
00538 gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;
00539
00540 ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );
00541
00542 free(qtemp);
00543 free(qprob);
00544 }
00545 else
00546 # endif
00547 {
00548
00549
00550
00551 double f = 1e-10;
00552
00553
00554 double sqrt_term = 35.399494936611667;
00555 double beta_alpha = 0.25 * sqrt_term *sexp(200./gv.bin[nd]->tedust );
00556
00557 double xi = 1./ (1. + 1.3e13*sexp(1.5*1e4/ gv.bin[nd]->tedust )*sqrt_term/(2.*f) );
00558
00559 double beta = 3e12 * sexp( 320. / gv.bin[nd]->tedust );
00560
00561
00562 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );
00563
00564
00565
00566
00567
00568 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.5 * DoppVel.AveVel[ipHYDROGEN]* gv.bin[nd]->IntArea/4. *
00569 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_CT02;
00570 ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );
00571 }
00572
00573 # ifndef IGNORE_QUANTUM_HEATING
00574
00575 sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) +
00576 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
00577 # endif
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589 hmi.rate_grain_h2_op_conserve += DoppVel.AveVel[LIMELM+2]* gv.bin[nd]->IntArea/4. *
00590 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H;
00591
00592
00593
00594
00595
00596
00597
00598
00599
00608
00609
00610
00611
00612
00613
00614
00615 T_ortho_para_crit = 2. * hmi.Tad / log( POW2(60. *1.1e11)*hmi.Tad);
00616 if( gv.bin[nd]->tedust < T_ortho_para_crit )
00617 {
00618 double efficiency_opr = sexp(60.*1.1e11*sqrt(hmi.Tad)*sexp(hmi.Tad/gv.bin[nd]->tedust));
00619
00620
00621 hmi.rate_grain_h2_J1_to_J0 += DoppVel.AveVel[LIMELM+2]* gv.bin[nd]->IntArea/4. *
00622 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * efficiency_opr;
00623 }
00624 }
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634 hmi.rate_grain_h2_op_conserve *= mole.lgH2_grain_deexcitation;
00635 hmi.rate_grain_h2_J1_to_J0 *= mole.lgH2_grain_deexcitation;
00636
00637 }
00638 else
00639 {
00640
00641 for( nd=0; nd < gv.nBin; nd++ )
00642 {
00643 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;
00644 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;
00645 }
00646
00647 hmi.rate_grain_h2_op_conserve = 0.;
00648
00649 hmi.rate_grain_h2_J1_to_J0 = 0.;
00650 }
00651
00652
00653
00654
00655
00656 gv.rate_h2_form_grains_used_total = 0.;
00657 for( nd=0; nd < gv.nBin; nd++ )
00658 {
00659 if( hmi.chJura == 'C' )
00660 {
00661
00662
00663
00664 gv.bin[nd]->rate_h2_form_grains_used =
00665 gv.bin[nd]->rate_h2_form_grains_CT02*hmi.ScaleJura;
00666 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
00667 }
00668 else if( hmi.chJura == 'T' )
00669 {
00670
00671 gv.bin[nd]->rate_h2_form_grains_used =
00672 gv.bin[nd]->rate_h2_form_grains_HM79*hmi.ScaleJura;
00673 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
00674 }
00675 else if( hmi.chJura == 'S' )
00676 {
00677
00678
00679 gv.bin[nd]->rate_h2_form_grains_used =
00680 3e-18 * phycon.sqrte / gv.nBin * dense.gas_phase[ipHYDROGEN]*hmi.ScaleJura;
00681
00682 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
00683 }
00684
00685
00686
00687 else if( hmi.chJura == 'F' )
00688 {
00689
00690
00691 gv.bin[nd]->rate_h2_form_grains_used = hmi.rate_h2_form_grains_set*dense.gas_phase[ipHYDROGEN] / gv.nBin;
00692 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
00693 }
00694 }
00695 ASSERT( gv.rate_h2_form_grains_used_total >= 0. );
00696
00697 # ifndef IGNORE_QUANTUM_HEATING
00698 printf( " fnzone %.2f H2 rate %.4e\n", fnzone, gv.rate_h2_form_grains_used_total );
00699 # endif
00700
00701
00702 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
00703 {
00704 frac_H2star_grains = hmi.H2star_forms_grains /
00705 SDIV(hmi.H2star_forms_grains+hmi.H2_forms_grains);
00706
00707 frac_H2star_hminus = hmi.H2star_forms_hminus /
00708 SDIV(hmi.H2star_forms_hminus+hmi.H2_forms_hminus);
00709
00710
00711
00712
00713
00714
00715 }
00716 else
00717 {
00718
00719
00720
00721
00722
00723
00724 frac_H2star_grains = 0.9416;
00725 frac_H2star_hminus = 1. - 4.938e-6;
00726 }
00727
00728
00729
00730
00731
00732
00733
00734
00735 corr = MIN2(6.,14.44-phycon.alogte*3.08);
00736
00737 if(corr > 0.)
00738 corr = pow(10.,corr*Hmolec_old[ipMH]/(Hmolec_old[ipMH]+1.6e4));
00739 else
00740 corr = 1.;
00741
00742 hmi.rh2dis = (realnum)(1.55e-8/phycon.sqrte*sexp(65107./phycon.te)* corr)*co.lgUMISTrates;
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767 hmi.bh2h2p = 6.4e-10f;
00768 # if 0
00769
00770
00771 if(hmi.rel_pop_LTE_H2g != 0.)
00772 hmi.rh2h2p = hmi.bh2h2p*hmi.rel_pop_LTE_H2p/hmi.rel_pop_LTE_H2g;
00773 else
00774 hmi.rh2h2p = 0;
00775 # endif
00776
00777
00778 if(phycon.te<=3.e4 )
00779 {
00780
00781 double teused = MAX2(100., phycon.te );
00782 double telog = log(teused);
00783 hmi.rh2h2p = sexp(2.123715e4/teused)*(-3.3232183e-7 + 3.3735382e-7*log(teused) -
00784 1.4491368e-7*pow(telog,2) +3.4172805e-8*pow(telog,3) -
00785 4.781372e-9*pow(telog,4) + 3.9731542e-10*pow(telog,5) -
00786 1.8171411e-11*pow(telog,6) +3.5311932e-13*pow(telog,7));
00787
00788 hmi.rh2h2p = hmi.rh2h2p*co.lgUMISTrates;
00789 }
00790 else
00791 hmi.rh2h2p= 0;
00792
00793
00794
00795
00796
00797
00798
00799 gamtwo = GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.);
00800
00801
00802 if(!co.lgUMISTrates)
00803 gamtwo = 5.7e-10*hmi.UV_Cont_rel2_Habing_TH85_face*(realnum)sexp((1.9*rfield.extin_mag_V_point))/1.66f;
00804
00805
00806
00807 h2phet = thermal.HeatNet;
00808
00809
00810
00811
00812
00813 sum_H0_Hp = Hmolec_old[ipMH]+Hmolec_old[ipMHp];
00814
00815 rindex = 0;
00816 r = rlist;
00817
00818 if(r == NULL)
00819 {
00820 int in[]={-1},out[]={-1};
00821 r = rlist = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00822 }
00823 rindex++;
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840 if(r->next == NULL) {
00841 int in[]={ipMH},out[]={ipMHm};
00842 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00843 }
00844 r = r->next;
00845 rindex++;
00846
00847
00848
00849
00850 r->rk = (hmi.hminus_rad_attach + hmi.HMinus_induc_rec_rate)*co.lgUMISTrates;
00851
00852
00853
00854
00855
00856
00857
00858
00859 if(phycon.te <= 7891.)
00860 {
00861
00862 hmihph2p = 6.9e-9 / (phycon.te30 * phycon.te05);
00863 }
00864 else
00865 {
00866
00867
00868 hmihph2p = 9.6e-7 / phycon.te90;
00869 }
00870
00871 if(r->next == NULL) {
00872 int in[]={ipMHm,ipMHp},out[]={ipMH2p};
00873 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00874 }
00875 r = r->next;
00876 rindex++;
00877
00878
00879
00880
00881 hmihph2p = hmihph2p*co.lgUMISTrates;
00882 r->rk = hmihph2p;
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893 TStancil = MIN2(phycon.te, 1e4 );
00894 TStancil = MAX2( TStancil , 10. );
00895
00897
00898
00899
00900 hmi.h2phmh2h = 1.4e-7*co.lgUMISTrates*17.305/phycon.sqrte;
00901 if(r->next == NULL) {
00902 int in[]={ipMH2p,ipMHm},out[]={ipMH2g,ipMH};
00903 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00904 }
00905 r = r->next;
00906 rindex++;
00907 r->rk = hmi.h2phmh2h;
00908
00909
00910
00911
00912
00913
00914
00915
00916 h2phmhhh = 1.4e-7f*17.3205/phycon.sqrte;
00917 if(r->next == NULL) {
00918 int in[]={ipMH2p,ipMHm},out[]={ipMH,ipMH,ipMH};
00919 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00920 }
00921 r = r->next;
00922 rindex++;
00923
00924 r->rk = h2phmhhh*co.lgUMISTrates;
00925
00926
00927
00928
00929
00930 hphmhhpe = 4.75e-30*pow(phycon.te,3.1);
00931 if(r->next == NULL) {
00932 int in[]={ipMHp,ipMHm},out[]={ipMH,ipMHp};
00933 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00934 }
00935 r = r->next;
00936 rindex++;
00937
00938 r->rk =hphmhhpe*co.lgUMISTrates;
00939
00941
00942
00943
00944
00945
00946 h2hmhh2e = 6.74e-17*co.lgUMISTrates*phycon.tesqrd*sexp(19870/phycon.te);
00947 if(r->next == NULL) {
00948 int in[]={ipMH2g,ipMHm},out[]={ipMH,ipMH2g};
00949 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00950 }
00951 r = r->next;
00952 rindex++;
00953 r->rk =h2hmhh2e;
00954
00955
00956
00957
00958 if(r->next == NULL) {
00959 int in[]={ipMHm},out[]={ipMH};
00960 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
00961 }
00962 r = r->next;
00963 rindex++;
00964 r->rk = hmi.HMinus_photo_rate;
00965
00966
00967
00968
00969
00970
00971 sum_first_ions = 0.;
00972 for( i=ipLITHIUM; i < LIMELM; i++ )
00973 {
00974 sum_first_ions += dense.xIonDense[i][1];
00975 }
00976
00977 {
00978
00979 enum {DEBUG_LOC=false};
00980 if( DEBUG_LOC && nzone>140 )
00981 {
00982 fprintf(ioQQQ,"DEBUG sumfirstions\t%.2f\t%.4e\t%.4e\t%.4e",
00983 fnzone,phycon.te,
00984 sum_first_ions,
00985 sum_first_ions/dense.eden);
00986 for( i=ipLITHIUM; i < LIMELM; i++ )
00987 {
00988 if( dense.xIonDense[i][1]/sum_first_ions >0.1 )
00989 fprintf(ioQQQ,"\t%li\t%.3e",
00990 i,dense.xIonDense[i][1]/sum_first_ions);
00991 }
00992 fprintf(ioQQQ,"\n");
00993 }
00994 }
00995
00996 hmi.hmin_ct_firstions = 4e-6f/(realnum)phycon.sqrte;
00997
00998 if(r->next == NULL) {
00999 int in[]={ipMHm},out[]={ipMH};
01000 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01001 }
01002 r = r->next;
01003 rindex++;
01004 r->rk = hmi.hmin_ct_firstions*sum_first_ions*co.lgUMISTrates;
01005
01006
01007 cionhm *= dense.eden;
01008
01009 if(r->next == NULL) {
01010 int in[]={ipMHm},out[]={ipMH};
01011 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01012 }
01013 r = r->next;
01014 rindex++;
01015 r->rk = cionhm*co.lgUMISTrates;
01016
01017
01018 c3bod = cionhm*(hmi.rel_pop_LTE_Hmin*dense.eden);
01019
01020 if(r->next == NULL) {
01021 int in[]={ipMH},out[]={ipMHm};
01022 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01023 }
01024 r = r->next;
01025 rindex++;
01026 r->rk = c3bod*co.lgUMISTrates;
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038 {
01039 double y , x;
01040 x = MAX2(10., phycon.te );
01041 x = MIN2(1e4, x );
01042 y=545969508.1323510+x*71239.23653059864;
01043 hmi.assoc_detach = 1./y;
01044 }
01045
01046
01047
01048 if(r->next == NULL) {
01049 int in[]={ipMH, ipMHm},out[]={ipMH2g};
01050 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01051 }
01052 r = r->next;
01053 rindex++;
01054 r->rk = hmi.assoc_detach*co.lgUMISTrates*(1.-frac_H2star_hminus);
01055
01056
01057
01058
01059 if(r->next == NULL) {
01060 int in[]={ipMH, ipMHm},out[]={ipMH2s};
01061 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01062 }
01063 r = r->next;
01064 rindex++;
01065 r->rk = hmi.assoc_detach*frac_H2star_hminus*co.lgUMISTrates;
01066
01067 {
01068
01069 enum {DEBUG_LOC=false};
01070 if( DEBUG_LOC && nzone>140 )
01071 {
01072 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",
01073 fnzone ,
01074 gv.rate_h2_form_grains_used_total ,
01075 hmi.H2_forms_grains+hmi.H2star_forms_grains ,
01076 frac_H2star_grains,
01077 hmi.assoc_detach*dense.xIonDense[ipHYDROGEN][0]*hmi.Hmolec[ipMHm],
01078 hmi.H2star_forms_hminus+hmi.H2_forms_hminus,
01079 frac_H2star_hminus,
01080 hmi.assoc_detach,dense.xIonDense[ipHYDROGEN][0],hmi.Hmolec[ipMHm]
01081 );
01082 }
01083 }
01084
01085
01086
01087 if( hmi.rel_pop_LTE_H2g > 0. )
01088 {
01089 hmi.assoc_detach_backwards_grnd = hmi.assoc_detach*hmi.rel_pop_LTE_Hmin/hmi.rel_pop_LTE_H2g*
01090 dense.eden*co.lgUMISTrates;
01091 }
01092 else
01093 {
01094 hmi.assoc_detach_backwards_grnd = 0.;
01095 }
01096
01097
01098
01099 if( hmi.rel_pop_LTE_H2s > 0. )
01100 {
01101
01102 hmi.assoc_detach_backwards_exct = hmi.assoc_detach*hmi.rel_pop_LTE_Hmin/hmi.rel_pop_LTE_H2s*
01103 dense.eden*co.lgUMISTrates;
01104 }
01105 else
01106 {
01107 hmi.assoc_detach_backwards_exct = 0.;
01108 }
01109
01110 {
01111
01112
01113
01114 enum {DEBUG_LOC=false};
01115 if( DEBUG_LOC && nzone>140)
01116 {
01117
01118 fprintf(ioQQQ,"hmi.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",
01119
01120 fnzone,
01121 phycon.te,
01122 dense.eden,
01123
01124 hmi.assoc_detach,
01125 hmi.assoc_detach_backwards_grnd,
01126 hmi.assoc_detach_backwards_exct,
01127 hmi.hminus_rad_attach,
01128 hmi.HMinus_induc_rec_rate,
01129
01130 hmi.Hmolec[ipMH],
01131
01132 hmi.Hmolec[ipMHp],
01133
01134 hmi.Hmolec[ipMHm],
01135 hmi.H2_total,
01136 hmi.rel_pop_LTE_Hmin,
01137 hmi.rel_pop_LTE_H2g,
01138 hmi.rel_pop_LTE_H2s
01139 );
01140 }
01141 }
01142
01143
01144 if(r->next == NULL) {
01145 int in[]={ipMH2g},out[]={ipMH,ipMHm};
01146 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01147 }
01148 r = r->next;
01149 rindex++;
01150
01151
01152
01153
01154
01155 hmi.assoc_detach_backwards_grnd *= ((1.-frac_H2star_hminus) * co.lgUMISTrates);
01156 r->rk = hmi.assoc_detach_backwards_grnd;
01157
01158
01159 if(r->next == NULL) {
01160 int in[]={ipMH2s},out[]={ipMH,ipMHm};
01161 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01162 }
01163 r = r->next;
01164 rindex++;
01165
01166
01167
01168
01169
01170
01171
01172 hmi.assoc_detach_backwards_exct *= (frac_H2star_hminus * co.lgUMISTrates);
01173 r->rk = hmi.assoc_detach_backwards_exct;
01174
01175
01176
01177
01178
01179 if( phycon.te < 14125. )
01180 {
01181
01182
01183 Hneut = 1.4e-7*pow(phycon.te/300,-0.487)*exp(phycon.te/29300);
01184 }
01185 else
01186 {
01187 Hneut = 3.4738192887404660e-008;
01188 }
01189
01190
01192 fhneut = Hmolec_old[ipMHp]*Hneut;
01193
01194 if(r->next == NULL) {
01195 int in[]={ipMHm,ipMHp},out[]={ipMH,ipMH};
01196 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01197 }
01198 r = r->next;
01199 rindex++;
01200
01201
01202
01203
01204 r->rk = Hneut*co.lgUMISTrates;
01205
01206
01207 if( phycon.te > 1000. )
01208 {
01209
01210 bhneut = (Hneut*hmi.rel_pop_LTE_Hmin*dense.eden)*iso.DepartCoef[ipH_LIKE][ipHYDROGEN][3];
01211 }
01212 else
01213 {
01214 bhneut = 0.;
01215 }
01216
01217
01218
01220
01221
01222 if(r->next == NULL) {
01223 int in[]={ipMH,ipMH},out[]={ipMHm,ipMHp}, ratesp[]={ipMH,ipMHp};
01224 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01225 }
01226 r = r->next;
01227 rindex++;
01228 r->rk = bhneut*co.lgUMISTrates;
01229 bhneut *= Hmolec_old[ipMHp];
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244 # define CATALYST true
01245 if( CATALYST )
01246 {
01247
01248
01249
01250
01251 if(r->next == NULL) {
01252 int in[]={ipMH,ipMH},out[]={ipMH2s},ratesp[]={ipMH};
01253 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01254 }
01255 r = r->next;
01256 rindex++;
01257 r->rk = gv.rate_h2_form_grains_used_total*frac_H2star_grains;
01258
01259
01260
01261
01262
01263 if(r->next == NULL) {
01264 int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH};
01265 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01266 }
01267 r = r->next;
01268 rindex++;
01269 r->rk = gv.rate_h2_form_grains_used_total*(1. - frac_H2star_grains);
01270 }
01271
01272 else
01273 {
01274
01275
01276
01277
01278
01279
01280
01281 if(r->next == NULL) {
01282 int in[]={ipMH,ipMH},out[]={ipMH2s},ratesp[]={ipMH};
01283 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01284 }
01285 r = r->next;
01286 rindex++;
01287 r->rk = gv.rate_h2_form_grains_used_total*frac_H2star_grains*
01288 dense.gas_phase[ipHYDROGEN]/SDIV(dense.xIonDense[ipHYDROGEN][0]);
01289
01290
01291
01292
01293
01294 if(r->next == NULL) {
01295 int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH};
01296 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01297 }
01298 r = r->next;
01299 rindex++;
01300 r->rk = gv.rate_h2_form_grains_used_total*(1. - frac_H2star_grains)*
01301 dense.gas_phase[ipHYDROGEN]/SDIV(dense.xIonDense[ipHYDROGEN][0]);
01302 }
01303
01304
01305
01306
01307
01308
01309
01310
01311
01313 hmi.radasc = ((StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop + StatesElem[ipH_LIKE][ipHYDROGEN][ipH2s].Pop))*3e-14;
01314
01315
01316
01317
01318 if(r->next == NULL) {
01319 int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH,ipMHp};
01320 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01321 }
01322 r = r->next;
01323 rindex++;
01324 r->rk = hmi.radasc*co.lgUMISTrates;
01325 hmi.radasc *= Hmolec_old[ipMHp];
01326
01327
01328
01329 if(r->next == NULL) {
01330 int in[]={ipMH2g},out[]={ipMH,ipMH};
01331 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01332 }
01333 r = r->next;
01334 rindex++;
01335
01336
01337
01338 r->rk = hmi.H2_Solomon_dissoc_rate_used_H2g;
01339
01340
01341
01342
01343
01344 if(r->next == NULL) {
01345 int in[]={ipMH2s},out[]={ipMH,ipMH};
01346 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01347 }
01348 r = r->next;
01349 rindex++;
01350
01351
01352
01353 r->rk = hmi.H2_Solomon_dissoc_rate_used_H2s;
01354
01355
01356
01357
01358
01359
01360
01361
01364 hmi.h2hph3p = 1.0e-16f*co.lgUMISTrates;
01365
01366 if(r->next == NULL) {
01367 int in[]={ipMH2g,ipMHp},out[]={ipMH3p};
01368 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01369 }
01370 r = r->next;
01371 rindex++;
01372 r->rk = hmi.h2hph3p;
01373
01374
01375
01376
01377
01378
01379
01380 if(r->next == NULL) {
01381 int in[]={ipMH2g},out[]={ipMH,ipMH},ratesp[]={ipMH,ipMH2g};
01382 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01383 }
01384 r = r->next;
01385 rindex++;
01386 hmi.rh2dis *= co.lgUMISTrates;
01387 r->rk = hmi.rh2dis;
01388
01389
01390
01391
01392
01393
01395 hmi.bh2h22hh2 = 5.5e-29*co.lgUMISTrates/(8.*phycon.te);
01396
01397 if(r->next == NULL) {
01398 int in[]={ipMH,ipMH,ipMH2g},out[]={ipMH2g,ipMH2g};
01399 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01400 }
01401 r = r->next;
01402 rindex++;
01403 r->rk = hmi.bh2h22hh2;
01404
01405
01406
01407
01408
01409
01411 if( hmi.rel_pop_LTE_H2g > 0. )
01412 {
01413 hmi.h2h22hh2 = hmi.bh2h22hh2/hmi.rel_pop_LTE_H2g;
01414 }
01415 else
01416 {
01417 hmi.h2h22hh2 =0.;
01418 }
01419
01420 if(r->next == NULL) {
01421 int in[]={ipMH2g,ipMH2g},out[]={ipMH,ipMH,ipMH2g};
01422 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01423 }
01424 r = r->next;
01425 rindex++;
01426 hmi.h2h22hh2 *= co.lgUMISTrates;
01427 r->rk = hmi.h2h22hh2;
01428
01429
01430
01431
01432
01433
01435 hmi.bh2dis = hmi.rh2dis*hmi.rel_pop_LTE_H2g*co.lgUMISTrates;
01436
01437 if(r->next == NULL) {
01438 int in[]={ipMH,ipMH},out[]={ipMH2g},ratesp[]={ipMH,ipMH,ipMH};
01439 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
01440 }
01441 r = r->next;
01442 rindex++;
01443 r->rk = hmi.bh2dis;
01444
01445 hmi.bh2dis = hmi.rh2dis*hmi.rel_pop_LTE_H2g*Hmolec_old[ipMH]*Hmolec_old[ipMH]*co.lgUMISTrates;
01446
01447
01448
01449
01450
01457
01458
01459
01460
01461
01463
01464
01465 {
01466 static long int nzone_eval = -1, iteration_evaluated=-1;
01467
01468 if( ( nzone_eval!=nzone || iteration_evaluated!=iteration ) || !nzone )
01469 {
01470
01471 hmi.H2_photoionize_rate =
01472 GammaK(opac.ipH2_photo_thresh ,
01473 rfield.nupper,
01474 opac.ipH2_photo_opac_offset,1.)*
01475 ionbal.lgPhotoIoniz_On +
01476
01477
01478
01479 2.*ionbal.CompRecoilIonRate[ipHYDROGEN][0];
01480
01481
01482
01483 hmi.H2_photo_heat_soft = thermal.HeatLowEnr * ionbal.lgPhotoIoniz_On;
01484 hmi.H2_photo_heat_hard = thermal.HeatHiEnr * ionbal.lgPhotoIoniz_On;
01485 nzone_eval = nzone;
01486 iteration_evaluated = iteration;
01487 }
01488 }
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504 if(r->next == NULL) {
01505 int in[]={ipMH2g},out[]={ipMH2p};
01506 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01507 }
01508 r = r->next;
01509 rindex++;
01510
01511
01512
01513
01514
01515
01516 if(co.lgUMISTrates)
01517 {
01518 h2crh2pe = hmi.H2_photoionize_rate + secondaries.csupra[ipHYDROGEN][0]*2.02;
01519 }
01520
01521 else
01522 {
01523 h2crh2pe = 4.4e-17;
01524 }
01525
01526 r->rk = h2crh2pe;
01527
01528
01529 if(r->next == NULL) {
01530 int in[]={ipMH2g},out[]={ipMH,ipMHp};
01531 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01532 }
01533 r = r->next;
01534 rindex++;
01535
01536
01537
01538
01539
01540
01541
01542 if(co.lgUMISTrates)
01543 {
01544 h2crhphe = secondaries.csupra[ipHYDROGEN][0]*0.0478;
01545 }
01546 else
01547 {
01548 h2crhphe = 1e-19;
01549 }
01550 r->rk = h2crhphe;
01551
01552
01553
01554 if(r->next == NULL) {
01555 int in[]={ipMH2s},out[]={ipMH,ipMHp};
01556 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01557 }
01558 r = r->next;
01559 rindex++;
01560 if(co.lgUMISTrates)
01561 {
01562 h2scrhphe = secondaries.csupra[ipHYDROGEN][0]*0.0478;
01563 }
01564 else
01565 {
01566 h2scrhphe = 1e-19;
01567 }
01568 r->rk = h2scrhphe;
01569
01570
01571
01572
01573
01574
01575 if(r->next == NULL) {
01576 int in[]={ipMH2s},out[]={ipMH2p};
01577 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01578 }
01579 r = r->next;
01580 rindex++;
01581
01582 if( co.lgUMISTrates )
01583 {
01584
01585 h2scrh2pe = hmi.H2_photoionize_rate + secondaries.csupra[ipHYDROGEN][0]*2.02;
01586 }
01587 else
01588 {
01589
01590 h2scrh2pe = 4.4e-17;
01591 }
01592 r->rk = h2scrh2pe;
01593
01594
01595
01596
01597
01598
01599 if(r->next == NULL) {
01600 int in[]={ipMH2g},out[]={ipMH,ipMH};
01601 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01602 }
01603 r = r->next;
01604 rindex++;
01605
01606
01607 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
01608 {
01609 h2crphh = hmi.H2_tripletdissoc_H2g;
01610 }
01611 else
01612 {
01613 h2crphh = secondaries.x12tot*3.;
01614 }
01615
01616
01617
01618
01619
01620
01621 if(!co.lgUMISTrates)
01622 h2crphh = 5e-18;
01623
01624 r->rk = h2crphh;
01625
01626
01627 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
01628 {
01629 h2scrphh = hmi.H2_tripletdissoc_H2s;
01630 }
01631 else
01632 {
01633 h2scrphh = secondaries.x12tot*3.;
01634 }
01635
01636 if(r->next == NULL) {
01637 int in[]={ipMH2s},out[]={ipMH,ipMH};
01638 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01639 }
01640 r = r->next;
01641 rindex++;
01642 r->rk = h2scrphh;
01643
01644
01645
01646
01647
01648
01649
01650 h2crphphm = 3.9e-21 * hextra.cryden_ov_background * co.lgUMISTrates;
01651
01652 if(r->next == NULL) {
01653 int in[]={ipMH2g},out[]={ipMHp,ipMHm};
01654 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01655 }
01656 r = r->next;
01657 rindex++;
01658
01659 r->rk = h2crphphm;
01660
01661
01662
01663 h2scrphphm = 3.9e-21 * hextra.cryden_ov_background * co.lgUMISTrates;
01664
01665 if(r->next == NULL) {
01666 int in[]={ipMH2s},out[]={ipMHp,ipMHm};
01667 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01668 }
01669 r = r->next;
01670 rindex++;
01671
01672 r->rk = h2scrphphm;
01673
01674
01675
01676
01677
01678
01679
01680
01681 h2crphpeh = 2.2e-19 * hextra.cryden_ov_background * co.lgUMISTrates;
01682
01683 if(r->next == NULL) {
01684 int in[]={ipMH2g},out[]={ipMHp,ipMH};
01685 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01686 }
01687 r = r->next;
01688 rindex++;
01689 r->rk = h2crphpeh;
01690
01691
01692
01693 h2scrphpeh = 2.2e-19 * hextra.cryden_ov_background * co.lgUMISTrates;
01694
01695 if(r->next == NULL) {
01696 int in[]={ipMH2s},out[]={ipMHp,ipMH};
01697 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01698 }
01699 r = r->next;
01700 rindex++;
01701 r->rk = h2scrphpeh;
01702
01703
01704
01705 hmi.CR_reac_H2g = h2crh2pe + h2crhphe + h2crphh + h2crphphm + h2crphpeh;
01706 hmi.CR_reac_H2s = h2scrh2pe + h2scrhphe + h2scrphh + h2scrphphm + h2scrphpeh;
01707
01708
01709
01710
01711
01712
01714
01715 hmi.h3phm2h2 = 1.3e-7 / (phycon.sqrte/31.62278) * co.lgUMISTrates;
01716 if(r->next == NULL) {
01717 int in[]={ipMH3p,ipMHm},out[]={ipMH2g,ipMH2g};
01718 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01719 }
01720 r = r->next;
01721 rindex++;
01722 r->rk = hmi.h3phm2h2;
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739 if(co.lgUMISTrates)
01740 {
01741 h3ph2ph = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_depth/1.66f;
01742 }
01743 else
01744 {
01745 h3ph2ph = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_face*(realnum)sexp((2.3*rfield.extin_mag_V_point))/1.66f;
01746 }
01747
01748 if(r->next == NULL) {
01749 int in[]={ipMH3p},out[]={ipMH2p,ipMH};
01750 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01751 }
01752 r = r->next;
01753 rindex++;
01754 r->rk =h3ph2ph;
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01770 if(co.lgUMISTrates)
01771 {
01772 hmi.h3ph2hp = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_depth/1.66f;
01773 }
01774 else
01775 {
01776 hmi.h3ph2hp = 5.0e-13*hmi.UV_Cont_rel2_Habing_TH85_face*sexp((1.8*rfield.extin_mag_V_point))/1.66;
01777 }
01778
01779 if(r->next == NULL) {
01780 int in[]={ipMH3p},out[]={ipMH2g,ipMHp};
01781 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01782 }
01783 r = r->next;
01784 rindex++;
01785 r->rk =hmi.h3ph2hp;
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795 if(r->next == NULL) {
01796 int in[]={ipMH2g,ipMHp},out[]={ipMH,ipMH2p};
01797 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01798 }
01799 r = r->next;
01800 rindex++;
01801 r->rk = hmi.rh2h2p;
01802
01803
01804
01805
01806
01807
01808
01809 if(r->next == NULL)
01810 {
01811
01812
01813 int in[]={ipMH,ipMH2p},out[]={ipMHp,ipMH2s};
01814 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01815 }
01816 r = r->next;
01817 rindex++;
01818 r->rk = hmi.bh2h2p;
01819
01821
01822
01823
01824
01825
01826
01828
01829
01830
01831 hmi.h3ph2p = HMRATE(2.08e-9,0.,1.88e4)*co.lgUMISTrates;
01832
01833 if(r->next == NULL)
01834 {
01835 int in[]={ipMH,ipMH3p},out[]={ipMH2g,ipMH2p};
01836 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01837 }
01838 r = r->next;
01839 rindex++;
01840 r->rk = hmi.h3ph2p;
01841
01842
01843
01844
01845
01846
01847
01848
01850 hmi.h3phmh2hh = 2.3e-7f*pow(phycon.te/300 , -0.5)*co.lgUMISTrates;
01851 if(r->next == NULL)
01852 {
01853 int in[]={ipMH3p,ipMHm},out[]={ipMH2g,ipMH,ipMH};
01854 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01855 }
01856 r = r->next;
01857 rindex++;
01858 r->rk = hmi.h3phmh2hh;
01859
01860
01861
01862
01864 hmi.h3petc = HMRATE(3.41e-11,0.5,7.16e4)*co.lgUMISTrates;
01865
01866 if(r->next == NULL)
01867 {
01868 int in[]={ipMH2g,ipMH3p},out[]={ipMH2g,ipMH2p,ipMH};
01869 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01870 }
01871 r = r->next;
01872 rindex++;
01873 r->rk = hmi.h3petc;
01874
01875
01876
01877
01879 hmi.h32h2 = HMRATE(3.41e-11,0.5,5.04e4)*co.lgUMISTrates;
01880
01881 if(r->next == NULL) {
01882 int in[]={ipMH2g,ipMH3p},out[]={ipMHp,ipMH2g,ipMH2g};
01883 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01884 }
01885 r = r->next;
01886 rindex++;
01887 r->rk = hmi.h32h2;
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01900
01901
01902
01903 # define USE_MCCALL
01904 # ifdef USE_MCCALL
01905 # define FACTOR 2.25
01906 # else
01907 # define FACTOR 1.0
01908 # endif
01909 hmi.eh3_h2h = HMRATE(4.00e-8/FACTOR,-0.5,0.)*dense.eden;
01910
01911
01912
01913 if(!co.lgUMISTrates)
01914 hmi.eh3_h2h = HMRATE(2.5e-8,-0.3,0.)*dense.eden;
01915
01916
01917 if(r->next == NULL) {
01918 int in[]={ipMH3p},out[]={ipMH,ipMH2g};
01919 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01920 }
01921 r = r->next;
01922 rindex++;
01923 r->rk = hmi.eh3_h2h;
01924
01925
01926
01927 eh3p_3h = HMRATE(1.6e-7/FACTOR,-0.5,0.)*dense.eden;
01928 # undef FACTOR
01929 # ifdef USE_MCCALL
01930 # undef USE_MCCALL
01931 # endif
01932
01933
01934
01935 if(!co.lgUMISTrates)
01936 eh3p_3h = HMRATE(7.5e-8,-0.3,0.)*dense.eden;
01937
01938
01939
01940
01941
01942 if(r->next == NULL) {
01943 int in[]={ipMH3p},out[]={ipMH,ipMH,ipMH};
01944 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
01945 }
01946 r = r->next;
01947 rindex++;
01948 r->rk = eh3p_3h;
01949
01950 if( (trace.lgTrace && trace.lgTr_H2_Mole) )
01951 {
01952 if( hmi.H2_rate_destroy > SMALLFLOAT )
01953 {
01954 fprintf( ioQQQ,
01955 " H2 destroy rate=%.2e DIS;%.3f bat;%.3f h2dis;%.3f hmi.H2_photoionize_rate;%.3f h2h2p;%.3f E-h;%.3f hmi.h2hph3p;%.3f sec;%.3f\n",
01956 hmi.H2_rate_destroy,
01957 hmi.H2_Solomon_dissoc_rate_used_H2g / hmi.H2_rate_destroy,
01958 hmi.assoc_detach_backwards_grnd / hmi.H2_rate_destroy,
01959 hmi.rh2dis*dense.xIonDense[ipHYDROGEN][0] / hmi.H2_rate_destroy,
01960 hmi.H2_photoionize_rate / hmi.H2_rate_destroy,
01961 hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1] / hmi.H2_rate_destroy,
01962 hmi.eh2hh /hmi.H2_rate_destroy,
01963 hmi.h2hph3p / hmi.H2_rate_destroy ,
01964 secondaries.csupra[ipHYDROGEN][0]*2.02 / hmi.H2_rate_destroy
01965 );
01966 }
01967 else
01968 {
01969 fprintf( ioQQQ, " Destroy H2: rate=0\n" );
01970 }
01971 }
01972
01973
01974
01975
01976
01978
01979
01980
01981
01983
01984
01985
01986
01987
01988
01989
01990
01991 radath = MAX2(0.,2.325*MIN2(5000.,phycon.te)-1375.)*1e-20;
01992
01993
01994
01995
01996 if( !co.lgUMISTrates)
01997 radath = HMRATE(5.3e-19,1.85,0);
01998
01999 if( r->next == NULL) {
02000 int in[]={ipMH,ipMHp},out[]={ipMH2p};
02001 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02002 }
02003 r = r->next;
02004
02005
02006 rindex++;
02007 r->rk = radath;
02008
02009
02010
02011
02012
02013
02014
02015 h2pion = 2.4e-27*POW3(phycon.te)*co.lgUMISTrates;
02016
02017 if(r->next == NULL) {
02018 int in[]={ipMH2p},out[]={ipMH,ipMHp},ratesp[]={ipMHp,ipMH2p};
02019 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
02020 }
02021 r = r->next;
02022 rindex++;
02023 r->rk = h2pion;
02024
02025
02026
02027
02028
02029
02030 h2pcin = 2e-7*sexp(30720./phycon.te)*dense.eden*co.lgUMISTrates;
02031
02032 if(r->next == NULL) {
02033 int in[]={ipMH2p},out[]={ipMH,ipMHp};
02034 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02035 }
02036 r = r->next;
02037 rindex++;
02038 r->rk = h2pcin;
02039
02040
02041
02042
02043 if( iteration==1 && conv.lgSearch )
02044 {
02045
02046 oatomic = 0.;
02047 oion = 0.;
02048 }
02049 else
02050 {
02051
02052 # define OLD_FRAC 0.0
02053 oatomic = oatomic*OLD_FRAC + dense.xIonDense[ipOXYGEN][0]*(1.-OLD_FRAC);
02054 oion = oion*OLD_FRAC + dense.xIonDense[ipOXYGEN][1]*(1.-OLD_FRAC);
02055
02056
02057
02058 oatomic *= ionbal.lgHO_ct_chem;
02059 oion *= ionbal.lgHO_ct_chem;
02060 }
02061
02062
02063
02064 if(r->next == NULL) {
02065 int in[]={ipMHp},out[]={ipMH};
02066 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02067 }
02068
02069 r = r->next;
02070 rindex++;
02071
02072 r->rk = atmdat.HCharExcIonOf[ipOXYGEN][0]*oatomic;
02073
02074
02075 if(r->next == NULL) {
02076 int in[]={ipMH},out[]={ipMHp};
02077 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02078 }
02079 r = r->next;
02080 rindex++;
02081
02082 r->rk = atmdat.HCharExcRecTo[ipOXYGEN][0]*oion;
02083
02084
02085
02086
02087
02088
02089
02090
02091 h2pehh = 2.8e-8*12.882/(phycon.te30*phycon.te07)*dense.eden;
02092
02093
02094
02095
02096 if(!co.lgUMISTrates)
02097 h2pehh = HMRATE(1.6e-8,-0.43,0);
02098
02099 if(r->next == NULL) {
02100 int in[]={ipMH2p},out[]={ipMH,ipMH};
02101 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02102 }
02103 r = r->next;
02104 rindex++;
02105 r->rk =h2pehh;
02106
02107
02108
02109
02110
02111
02112
02113
02114 b2pcin = h2pcin*hmi.rel_pop_LTE_H2p;
02115
02116
02117 if(r->next == NULL) {
02118 int in[]={ipMH,ipMHp},out[]={ipMH2p};
02119 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02120 }
02121 r = r->next;
02122 rindex++;
02123 r->rk = b2pcin;
02124
02125
02126
02127 if(r->next == NULL) {
02128 int in[]={ipMH2p},out[]={ipMH,ipMHp};
02129 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02130 }
02131 r = r->next;
02132 rindex++;
02133 r->rk = gamtwo;
02134
02135
02136
02137
02138
02139 if(r->next == NULL) {
02140 int in[]={ipMH2p},out[]={ipMH,ipMHp};
02141 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02142 }
02143 r = r->next;
02144 rindex++;
02145
02146
02147
02148
02149
02150
02151
02152
02153 r->rk = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*co.lgUMISTrates;
02154
02155
02158 hmi.h2ph3p = 1.40e-9*(1. - sexp(9940./phycon.te));
02159
02160 if(!co.lgUMISTrates)
02161 hmi.h2ph3p = 2.08e-9;
02162
02163 if(r->next == NULL) {
02164 int in[]={ipMH2g,ipMH2p},out[]={ipMH,ipMH3p};
02165 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02166 }
02167 r = r->next;
02168 rindex++;
02169 r->rk = hmi.h2ph3p;
02170
02171
02173
02174
02175
02176 h2phhp = 2.41e-12*phycon.sqrte*sexp(30720./phycon.te)*co.lgUMISTrates;
02177
02178 if(r->next == NULL) {
02179 int in[]={ipMH2g,ipMH2p},out[]={ipMH,ipMHp,ipMH2g};
02180 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02181 }
02182 r = r->next;
02183 rindex++;
02184 r->rk = h2phhp;
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194 if(r->next == NULL) {
02195 int in[]={ipMH3p},out[]={ipMH2p,ipMHp};
02196 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02197 }
02198 r = r->next;
02199 rindex++;
02200
02201
02202
02203 r->rk = 2.*iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]*co.lgUMISTrates;
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215 Boltz_fac_H2_H2star = 1.*sexp( 30172./phycon.te);
02216 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02217 {
02218 deexc_htwo = hmi.Average_collH2_deexcit;
02219 deexc_hneut = hmi.Average_collH_deexcit;
02220 }
02221 else
02222 {
02223 deexc_htwo = (1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6.;
02224 deexc_hneut = (1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.;
02225 }
02226
02227 H2star_deexcit = hmi.H2_total*deexc_htwo + hmi.Hmolec[ipMH] * deexc_hneut;
02228
02229
02230 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02231 {
02232 H2star_excit = hmi.Average_collH2_excit *hmi.H2_total +
02233 hmi.Average_collH_excit*hmi.Hmolec[ipMH];
02234 }
02235 else
02236 {
02237 H2star_excit = Boltz_fac_H2_H2star * H2star_deexcit;
02238 }
02239
02240
02241
02242
02243
02244 if(r->next == NULL) {
02245 int in[]={ipMH2s},out[]={ipMH2g};
02246 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02247 }
02248 r = r->next;
02249 rindex++;
02250
02251
02252
02253 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02254 {
02255 hmi.h2s_sp_decay = hmi.Average_A;
02256 }
02257 else
02258 {
02259 hmi.h2s_sp_decay = 2e-7;
02260 }
02261 r->rk = hmi.h2s_sp_decay;
02262
02263
02264 if(r->next == NULL) {
02265 int in[]={ipMH2s},out[]={ipMH2g},ratesp[]={ipMH2s,ipMH2g};
02266 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
02267 }
02268 r = r->next;
02269 rindex++;
02270 r->rk = deexc_htwo;
02271
02272 if(r->next == NULL) {
02273 int in[]={ipMH2s},out[]={ipMH2g},ratesp[]={ipMH2s,ipMH};
02274 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
02275 }
02276 r = r->next;
02277 rindex++;
02278 r->rk = deexc_hneut;
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291 if(r->next == NULL) {
02292 int in[]={ipMH2g},out[]={ipMH2s},ratesp[]={ipMH2g,ipMH2g};
02293 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
02294 }
02295 r = r->next;
02296 rindex++;
02297
02298 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02299 {
02300 r->rk = hmi.Average_collH2_excit;
02301 }
02302 else
02303 {
02304 r->rk = deexc_htwo*Boltz_fac_H2_H2star;
02305 }
02306
02307 if(r->next == NULL) {
02308 int in[]={ipMH2g},out[]={ipMH2s},ratesp[]={ipMH,ipMH2g};
02309 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),ratesp,INTSZ(ratesp));
02310 }
02311 r = r->next;
02312 rindex++;
02313
02314 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02315 {
02316 r->rk = hmi.Average_collH_excit;
02317 }
02318 else
02319 {
02320 r->rk = deexc_hneut*Boltz_fac_H2_H2star;
02321 }
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331 hmi.h2sh = HMRATE(4.67e-7,-1.,5.5e4);
02332
02333 if(r->next == NULL) {
02334 int in[]={ipMH2s,ipMH},out[]={ipMH,ipMH,ipMH};
02335 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02336 }
02337 r = r->next;
02338 rindex++;
02339 r->rk = hmi.h2sh;
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412 hmi.h2sh2g = HMRATE(1e-11,0.,2.18e4);
02413 if(r->next == NULL) {
02414 int in[]={ipMH2s,ipMH2g},out[]={ipMH2g,ipMH,ipMH};
02415 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02416 }
02417 r = r->next;
02418 rindex++;
02419
02420 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02421 {
02422 hmi.h2sh2g = hmi.Average_collH2s_dissoc;
02423 }
02424
02425 r->rk = hmi.h2sh2g;
02426
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436 hmi.h2sh2sh2g2h = HMRATE(1e-11,0.,0.);
02437 if(r->next == NULL) {
02438 int in[]={ipMH2s,ipMH2s},out[]={ipMH2g,ipMH,ipMH};
02439 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02440 }
02441 r = r->next;
02442 rindex++;
02443
02444 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02445 {
02446 hmi.h2sh2sh2g2h = hmi.Average_collH2s_dissoc;
02447 }
02448
02449 r->rk = hmi.h2sh2sh2g2h;
02450
02451
02452
02453
02454
02455
02456 hmi.h2sh2sh2s2h = HMRATE(1e-11,0.,2.18e4);
02457 if(r->next == NULL) {
02458 int in[]={ipMH2s,ipMH2s},out[]={ipMH2s,ipMH,ipMH};
02459 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02460 }
02461 r = r->next;
02462 rindex++;
02463 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
02464 {
02465 hmi.h2sh2sh2s2h = hmi.Average_collH2s_dissoc;
02466 }
02467
02468 r->rk = hmi.h2sh2sh2s2h;
02469
02470
02471
02472
02473
02474
02475
02476
02477 if(r->next == NULL) {
02478 int in[]={ipMH2g},out[]={ipMH2s};
02479 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02480 }
02481 r = r->next;
02482 rindex++;
02483
02484
02485
02486
02487 r->rk = hmi.H2_H2g_to_H2s_rate_used;
02488
02489
02490 if(r->next == NULL) {
02491 int in[]={ipMH2s},out[]={ipMH,ipMH};
02492 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02493 }
02494 r = r->next;
02495 rindex++;
02496 r->rk = hmi.H2_photodissoc_used_H2s;
02497
02498
02499
02500 if(r->next == NULL) {
02501 int in[]={ipMH2g},out[]={ipMH,ipMH};
02502 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02503 }
02504 r = r->next;
02505 rindex++;
02506 r->rk = hmi.H2_photodissoc_used_H2g;
02507
02508
02509
02510 if(r->next == NULL) {
02511 int in[]={ipMH2g},out[]={ipMH,ipMH};
02512 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02513 }
02514 r = r->next;
02515 rindex++;
02516 hmi.h2ge2h = 1e-14*sexp(4.478*EVDEGK/phycon.te)*dense.eden;
02517 r->rk = hmi.h2ge2h;
02518
02519
02520
02521 if(r->next == NULL) {
02522 int in[]={ipMH2s},out[]={ipMH,ipMH};
02523 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02524 }
02525 r = r->next;
02526 rindex++;
02527 hmi.h2se2h = 1e-14*sexp(1.978*EVDEGK/phycon.te)*dense.eden;
02528 r->rk = hmi.h2se2h;
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538 if(r->next == NULL) {
02539 int in[]={ipMH},out[]={ipMHeHp};
02540 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02541 }
02542 r = r->next;
02543 rindex++;
02544 r->rk = 1e-15*dense.xIonDense[ipHELIUM][1];
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562
02563
02564 hehmeheh = 4.1e-17*phycon.tesqrd*sexp(19870/phycon.te)*dense.xIonDense[ipHELIUM][0];
02565 if(r->next == NULL) {
02566 int in[]={ipMHm},out[]={ipMH};
02567 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02568 }
02569 r = r->next;
02570 rindex++;
02571 r->rk =hehmeheh;
02572
02573
02574
02575
02576
02577
02578
02579
02580
02581
02582
02583
02584
02586 hmi.rheph2hpheh = (realnum)HMRATE(3.7e-14, 0., 35);
02587 if(r->next == NULL) {
02588 int in[]={ipMH2g},out[]={ipMHp,ipMH};
02589 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02590 }
02591 r = r->next;
02592 rindex++;
02593 r->rk = hmi.rheph2hpheh*dense.xIonDense[ipHELIUM][1];
02594
02595
02596
02597
02598
02600 hmi.heph2heh2p = (realnum)HMRATE(7.2e-15, 0., 0);
02601 if(r->next == NULL) {
02602 int in[]={ipMH2g},out[]={ipMH2p};
02603 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02604 }
02605 r = r->next;
02606 rindex++;
02607 r->rk = hmi.heph2heh2p*dense.xIonDense[ipHELIUM][1];
02608
02609
02610 if(r->next == NULL) {
02611 int in[]={ipMHp},out[]={ipMHeHp};
02612 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02613 }
02614 r = r->next;
02615 rindex++;
02616 r->rk = 1e-20*dense.xIonDense[ipHELIUM][0];
02617
02618
02619 if(r->next == NULL) {
02620 int in[]={ipMH2p},out[]={ipMH,ipMHeHp};
02621 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02622 }
02623 r = r->next;
02624 rindex++;
02625
02626
02627
02628
02629
02630 r->rk = 3e-10*exp(-6717./phycon.te)*dense.xIonDense[ipHELIUM][0]*co.lgUMISTrates;
02631
02632
02633
02634
02635
02636 gamheh = 0.;
02637 limit = MIN2(hmi.iheh2-1 , rfield.nflux );
02638 for( i=hmi.iheh1-1; i < limit; i++ )
02639 {
02640 gamheh += rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i] + rfield.outlin_noplot[i];
02641 }
02642 gamheh *= 4e-18;
02643
02644
02645 gamheh += 3.*iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s];
02646
02647
02648 gamheh += dense.eden*1e-9;
02649
02650 if(r->next == NULL) {
02651 int in[]={ipMHeHp},out[]={ipMH};
02652 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02653 }
02654 r = r->next;
02655 rindex++;
02656 r->rk = gamheh;
02657
02658
02659 if(r->next == NULL) {
02660 int in[]={ipMH,ipMHeHp},out[]={ipMH2p};
02661 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02662 }
02663 r = r->next;
02664 rindex++;
02665
02666
02667
02668
02669
02670 r->rk = 1e-10*co.lgUMISTrates;
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02684 hmi.hehph2h3phe = 1.3e-9*co.lgUMISTrates;
02685 if(r->next == NULL) {
02686 int in[]={ipMH2g,ipMHeHp},out[]={ipMH3p};
02687 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02688 }
02689 r = r->next;
02690 rindex++;
02691 r->rk = hmi.hehph2h3phe;
02692
02693
02694
02695
02696
02697
02698
02699 hephmhhe = 2.32e-7*pow(phycon.te/300,-.52)*exp(TStancil/22400.)*dense.xIonDense[ipHELIUM][1];
02700 if(r->next == NULL) {
02701 int in[]={ipMHm},out[]={ipMH};
02702 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02703 }
02704 r = r->next;
02705 rindex++;
02706 r->rk =hephmhhe;
02707
02708 # define CO_ON true
02709
02710
02711
02712 if( CO_ON && !dynamics.lgAdvection )
02713 {
02714
02715
02716
02717
02718
02719
02720 co.H_CH_C_H_H = CO_findrk("H,CH=>C,H,H")*findspecies("CH")->hevmol;
02721 co.H_OH_O_H_H = CO_findrk("H,OH=>O,H,H")*findspecies("OH")->hevmol;
02722 co.H_H2O_OH_H_H = CO_findrk("H,H2O=>OH,H,H")*findspecies("H2O")->hevmol;
02723 co.H_COP_CO_HP = CO_findrk("H,CO+=>CO,H+")*findspecies("CO+")->hevmol;
02724 co.H_CH_C_H2 = CO_findrk("H,CH=>C,H2")*findspecies("CH")->hevmol;
02725 co.H_CHP_CP_H2 = CO_findrk("H,CH+=>C+,H2")*findspecies("CH+")->hevmol;
02726 co.H_CH2_CH_H2 = CO_findrk("H,CH2=>CH,H2")*findspecies("CH2")->hevmol;
02727 co.H_CH3P_CH2P_H2 = CO_findrk("H,CH3+=>CH2+,H2")*findspecies("CH3+")->hevmol;
02728 co.H_OH_O_H2 = CO_findrk("H,OH=>O,H2")*findspecies("OH")->hevmol;
02729 co.H_H2O_OH_H2 = CO_findrk("H,H2O=>OH,H2")*findspecies("H2O")->hevmol;
02730 co.Hminus_HCOP_CO_H2 = CO_findrk("H-,HCO+=>CO,H2")*findspecies("HCO+")->hevmol;
02731 co.Hminus_H3OP_H2O_H2 = CO_findrk("H-,H3O+=>H2O,H2")*findspecies("H3O+")->hevmol;
02732 co.Hminus_H3OP_OH_H2_H = CO_findrk("H-,H3O+=>OH,H2,H")*findspecies("H3O+")->hevmol;
02733 co.HP_CH_CHP_H = CO_findrk("H+,CH=>CH+,H")*findspecies("CH")->hevmol;
02734 co.HP_CH2_CH2P_H = CO_findrk("H+,CH2=>CH2+,H")*findspecies("CH2")->hevmol;
02735 co.HP_H2O_H2OP_H = CO_findrk("H+,H2O=>H2O+,H")*findspecies("H2O")->hevmol;
02736 co.HP_O2_O2P_H = CO_findrk("H+,O2=>O2+,H")*findspecies("O2")->hevmol;
02737 co.HP_OH_OHP_H = CO_findrk("H+,OH=>OH+,H")*findspecies("OH")->hevmol;
02738 co.HP_SiO_SiOP_H = CO_findrk("H+,SiO=>SiO+,H")*findspecies("SiO")->hevmol;
02739 co.HP_CH2_CHP_H2 = CO_findrk("H+,CH2=>CH+,H2")*findspecies("CH2")->hevmol;
02740 co.HP_SiH_SiP_H2 = CO_findrk("H+,SiH=>Si+,H2")*findspecies("SiH")->hevmol;
02741 co.H2_CHP_CH2P_H = CO_findrk("H2,CH+=>CH2+,H")*findspecies("CH+")->hevmol;
02742 co.H2_CH2P_CH3P_H = CO_findrk("H2,CH2+=>CH3+,H")*findspecies("CH2+")->hevmol;
02743 co.H2_OHP_H2OP_H = CO_findrk("H2,OH+=>H2O+,H")*findspecies("OH+")->hevmol;
02744 co.H2_H2OP_H3OP_H = CO_findrk("H2,H2O+=>H3O+,H")*findspecies("H2O+")->hevmol;
02745 co.H2_COP_HCOP_H = CO_findrk("H2,CO+=>HCO+,H")*findspecies("CO+")->hevmol;
02746 co.H2_OP_OHP_H = CO_findrk("H2,O+=>OH+,H")*findspecies("O+")->hevmol;
02747 co.H2_SiOP_SiOHP_H = CO_findrk("H2,SiO+=>SiOH+,H")*findspecies("SiO+")->hevmol;
02748 co.H2_C_CH_H = CO_findrk("H2,C=>CH,H")*findspecies("C")->hevmol;
02749 co.H2_CP_CHP_H = CO_findrk("H2,C+=>CH+,H")*findspecies("C+")->hevmol;
02750 co.H2_CH_CH2_H = CO_findrk("H2,CH=>CH2,H")*findspecies("CH")->hevmol;
02751 co.H2_OH_H2O_H = CO_findrk("H2,OH=>H2O,H")*findspecies("OH")->hevmol;
02752 co.H2_O_OH_H = CO_findrk("H2,O=>OH,H")*findspecies("O")->hevmol;
02753 co.H2_CH_C_H2_H = CO_findrk("H2,CH=>C,H2,H")*findspecies("CH")->hevmol;
02754 co.H2_OH_O_H2_H = CO_findrk("H2,OH=>O,H2,H")*findspecies("OH")->hevmol;
02755 co.H2_H2O_OH_H2_H = CO_findrk("H2,H2O=>OH,H2,H")*findspecies("H2O")->hevmol;
02756 co.H2_O2_O_O_H2 = CO_findrk("H2,O2=>O,O,H2")*findspecies("O2")->hevmol;
02757 co.H2s_CH_C_H2_H = CO_findrk("H2*,CH=>C,H2,H")*findspecies("CH")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02758 co.H2s_OH_O_H2_H = CO_findrk("H2*,OH=>O,H2,H")*findspecies("OH")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02759 co.H2s_H2O_OH_H2_H = CO_findrk("H2*,H2O=>OH,H2,H")*findspecies("H2O")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02760 co.H2s_O2_O_O_H2 = CO_findrk("H2*,O2=>O,O,H2")*findspecies("O2")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02761 co.H2P_C_CHP_H = CO_findrk("H2+,C=>CH+,H")*findspecies("C")->hevmol;
02762 co.H2P_CH_CH2P_H = CO_findrk("H2+,CH=>CH2+,H")*findspecies("CH")->hevmol;
02763 co.H2P_CH2_CH3P_H = CO_findrk("H2+,CH2=>CH3+,H")*findspecies("CH2")->hevmol;
02764 co.H2P_OH_H2OP_H = CO_findrk("H2+,OH=>H2O+,H")*findspecies("OH")->hevmol;
02765 co.H2P_H2O_H3OP_H = CO_findrk("H2+,H2O=>H3O+,H")*findspecies("H2O")->hevmol;
02766 co.H2P_CO_HCOP_H = CO_findrk("H2+,CO=>HCO+,H")*findspecies("CO")->hevmol;
02767 co.H2P_O_OHP_H = CO_findrk("H2+,O=>OH+,H")*findspecies("O")->hevmol;
02768 co.H2P_CH_CHP_H2 = CO_findrk("H2+,CH=>CH+,H2")*findspecies("CH")->hevmol;
02769 co.H2P_CH2_CH2P_H2 = CO_findrk("H2+,CH2=>CH2+,H2")*findspecies("CH2")->hevmol;
02770 co.H2P_CO_COP_H2 = CO_findrk("H2+,CO=>CO+,H2")*findspecies("CO")->hevmol;
02771 co.H2P_H2O_H2OP_H2 = CO_findrk("H2+,H2O=>H2O+,H2")*findspecies("H2O")->hevmol;
02772 co.H2P_O2_O2P_H2 = CO_findrk("H2+,O2=>O2+,H2")*findspecies("O2")->hevmol;
02773 co.H2P_OH_OHP_H2 = CO_findrk("H2+,OH=>OH+,H2")*findspecies("OH")->hevmol;
02774 co.H3P_C_CHP_H2 = CO_findrk("H3+,C=>CH+,H2")*findspecies("C")->hevmol;
02775 co.H3P_CH_CH2P_H2 = CO_findrk("H3+,CH=>CH2+,H2")*findspecies("CH")->hevmol;
02776 co.H3P_CH2_CH3P_H2 = CO_findrk("H3+,CH2=>CH3+,H2")*findspecies("CH2")->hevmol;
02777 co.H3P_OH_H2OP_H2 = CO_findrk("H3+,OH=>H2O+,H2")*findspecies("OH")->hevmol;
02778 co.H3P_H2O_H3OP_H2 = CO_findrk("H3+,H2O=>H3O+,H2")*findspecies("H2O")->hevmol;
02779 co.H3P_CO_HCOP_H2 = CO_findrk("H3+,CO=>HCO+,H2")*findspecies("CO")->hevmol;
02780 co.H3P_O_OHP_H2 = CO_findrk("H3+,O=>OH+,H2")*findspecies("O")->hevmol;
02781 co.H3P_SiH_SiH2P_H2 = CO_findrk("H3+,SiH=>SiH2+,H2")*findspecies("SiH")->hevmol;
02782 co.H3P_SiO_SiOHP_H2 = CO_findrk("H3+,SiO=>SiOH+,H2")*findspecies("SiO")->hevmol;
02783 co.H2s_CH_CH2_H = CO_findrk("H2*,CH=>CH2,H")*findspecies("CH")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02784 co.H2s_O_OH_H = CO_findrk("H2*,O=>OH,H")*findspecies("O")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02785 co.H2s_OH_H2O_H = CO_findrk("H2*,OH=>H2O,H")*findspecies("OH")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02786 co.H2s_C_CH_H = CO_findrk("H2*,C=>CH,H")*findspecies("C")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02787 co.H2s_CP_CHP_H = CO_findrk("H2*,C+=>CH+,H")*findspecies("C+")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02788 co.H_CH3_CH2_H2 = CO_findrk("H,CH3=>CH2,H2")*findspecies("CH3")->hevmol;
02789 co.H_CH4P_CH3P_H2 = CO_findrk("H,CH4+=>CH3+,H2")*findspecies("CH4+")->hevmol;
02790 co.H_CH5P_CH4P_H2 = CO_findrk("H,CH5+=>CH4+,H2")*findspecies("CH5+")->hevmol;
02791 co.H2_CH2_CH3_H = CO_findrk("H2,CH2=>CH3,H")*findspecies("CH2")->hevmol;
02792 co.H2_CH3_CH4_H = CO_findrk("H2,CH3=>CH4,H")*findspecies("CH3")->hevmol;
02793 co.H2_CH4P_CH5P_H = CO_findrk("H2,CH4+=>CH5+,H")*findspecies("CH4+")->hevmol;
02794 co.H2s_CH2_CH3_H = CO_findrk("H2*,CH2=>CH3,H")*findspecies("CH2")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02795 co.H2s_CH3_CH4_H = CO_findrk("H2*,CH3=>CH4,H")*findspecies("CH3")->hevmol*hmi.lgLeiden_Keep_ipMH2s;
02796 co.H2P_CH4_CH3P_H2 = CO_findrk("H2+,CH4=>CH3+,H2,H")*findspecies("CH4")->hevmol;
02797 co.H2P_CH4_CH4P_H2 = CO_findrk("H2+,CH4=>CH4+,H2")*findspecies("CH4")->hevmol;
02798 co.H2P_CH4_CH5P_H = CO_findrk("H2+,CH4=>CH5+,H")*findspecies("CH4")->hevmol;
02799
02800 co.H3P_CH3_CH4P_H2 = CO_findrk("H3+,CH3=>CH4+,H2")*findspecies("CH3")->hevmol;
02801 co.H3P_CH4_CH5P_H2 = CO_findrk("H3+,CH4=>CH5+,H2")*findspecies("CH4")->hevmol;
02802 co.HP_CH3_CH3P_H = CO_findrk("H+,CH3=>CH3+,H")*findspecies("CH3")->hevmol;
02803 co.HP_CH4_CH3P_H2 = CO_findrk("H+,CH4=>CH3+,H2")*findspecies("CH4")->hevmol;
02804 co.HP_CH4_CH4P_H = CO_findrk("H+,CH4=>CH4+,H")*findspecies("CH4")->hevmol;
02805
02806
02807
02808
02809
02810 co.H2_N_NH_H = CO_findrk("H2,N=>NH,H")*findspecies("N")->hevmol;
02811 co.H2_NH_NH2_H = CO_findrk("H2,NH=>NH2,H")*findspecies("NH")->hevmol;
02812 co.H2_NH2_NH3_H = CO_findrk("H2,NH2=>NH3,H")*findspecies("NH2")->hevmol;
02813 co.H2_CN_HCN_H = CO_findrk("H2,CN=>HCN,H")*findspecies("CN")->hevmol;
02814 co.HP_HNO_NOP_H2 = CO_findrk("H+,HNO=>NO+,H2")*findspecies("HNO")->hevmol;
02815 co.HP_HS_SP_H2 = CO_findrk("H+,HS=>S+,H2")*findspecies("HS")->hevmol;
02816 co.H_HSP_SP_H2 = CO_findrk("H,HS+=>S+,H2")*findspecies("HS+")->hevmol;
02817 co.H2P_N_NHP_H = CO_findrk("H2+,N=>NH+,H")*findspecies("N")->hevmol;
02818 co.H2_NP_NHP_H = CO_findrk("H2,N+=>NH+,H")*findspecies("N+")->hevmol;
02819 co.H2_NHP_N_H3P = CO_findrk("H2,NH+=>N,H3+")*findspecies("NH+")->hevmol;
02820 co.H2P_NH_NH2P_H = CO_findrk("H2+,NH=>NH2+,H")*findspecies("NH")->hevmol;
02821 co.H2_NHP_NH2P_H = CO_findrk("H2,NH+=>NH2+,H")*findspecies("NH+")->hevmol;
02822 co.H2_NH2P_NH3P_H = CO_findrk("H2,NH2+=>NH3+,H")*findspecies("NH2+")->hevmol;
02823 co.H2_NH3P_NH4P_H = CO_findrk("H2,NH3+=>NH4+,H")*findspecies("NH3+")->hevmol;
02824 co.H2P_CN_HCNP_H = CO_findrk("H2+,CN=>HCN+,H")*findspecies("CN")->hevmol;
02825 co.H2_CNP_HCNP_H = CO_findrk("H2,CN+=>HCN+,H")*findspecies("CN+")->hevmol;
02826 co.H2P_NO_HNOP_H = CO_findrk("H2+,NO=>HNO+,H")*findspecies("NO")->hevmol;
02827 co.H2_SP_HSP_H = CO_findrk("H2,S+=>HS+,H")*findspecies("S+")->hevmol;
02828 co.H2_CSP_HCSP_H = CO_findrk("H2,CS+=>HCS+,H")*findspecies("CS+")->hevmol;
02829 co.H3P_NH_NH2P_H2 = CO_findrk("H3+,NH=>NH2+,H2")*findspecies("NH")->hevmol;
02830 co.H3P_NH2_NH3P_H2 = CO_findrk("H3+,NH2=>NH3+,H2")*findspecies("NH2")->hevmol;
02831 co.H3P_NH3_NH4P_H2 = CO_findrk("H3+,NH3=>NH4+,H2")*findspecies("NH3")->hevmol;
02832 co.H3P_CN_HCNP_H2 = CO_findrk("H3+,CN=>HCN+,H2")*findspecies("CN")->hevmol;
02833 co.H3P_NO_HNOP_H2 = CO_findrk("H3+,NO=>HNO+,H2")*findspecies("NO")->hevmol;
02834 co.H3P_S_HSP_H2 = CO_findrk("H3+,S=>HS+,H2")*findspecies("S")->hevmol;
02835 co.H3P_CS_HCSP_H2 = CO_findrk("H3+,CS=>HCS+,H2")*findspecies("CS")->hevmol;
02836 co.H3P_NO2_NOP_OH_H2 = CO_findrk("H3+,NO2=>NO+,OH,H2")*findspecies("NO2")->hevmol;
02837 co.HP_NH_NHP_H = CO_findrk("H+,NH=>NH+,H")*findspecies("NH")->hevmol;
02838 co.HP_NH2_NH2P_H = CO_findrk("H+,NH2=>NH2+,H")*findspecies("NH2")->hevmol;
02839 co.HP_NH3_NH3P_H = CO_findrk("H+,NH3=>NH3+,H")*findspecies("NH3")->hevmol;
02840 co.H_CNP_CN_HP = CO_findrk("H,CN+=>CN,H+")*findspecies("CN+")->hevmol;
02841 co.HP_HCN_HCNP_H = CO_findrk("H+,HCN=>HCN+,H")*findspecies("HCN")->hevmol;
02842 co.H_HCNP_HCN_HP = CO_findrk("H,HCN+=>HCN,H+")*findspecies("HCN+")->hevmol;
02843 co.H_N2P_N2_HP = CO_findrk("H,N2+=>N2,H+")*findspecies("N2+")->hevmol;
02844 co.HP_NO_NOP_H = CO_findrk("H+,NO=>NO+,H")*findspecies("NO")->hevmol;
02845 co.HP_HS_HSP_H = CO_findrk("H+,HS=>HS+,H")*findspecies("HS")->hevmol;
02846 co.HP_SiN_SiNP_H = CO_findrk("H+,SiN=>SiN+,H")*findspecies("SiN")->hevmol;
02847 co.HP_CS_CSP_H = CO_findrk("H+,CS=>CS+,H")*findspecies("CS")->hevmol;
02848 co.HP_NS_NSP_H = CO_findrk("H+,NS=>NS+,H")*findspecies("NS")->hevmol;
02849 co.HP_SO_SOP_H = CO_findrk("H+,SO=>SO+,H")*findspecies("SO")->hevmol;
02850 co.HP_OCS_OCSP_H = CO_findrk("H+,OCS=>OCS+,H")*findspecies("OCS")->hevmol;
02851 co.HP_S2_S2P_H = CO_findrk("H+,S2=>S2+,H")*findspecies("S2")->hevmol;
02852 co.H2P_NH_NHP_H2 = CO_findrk("H2+,NH=>NH+,H2")*findspecies("NH")->hevmol;
02853 co.H2P_NH2_NH2P_H2 = CO_findrk("H2+,NH2=>NH2+,H2")*findspecies("NH2")->hevmol;
02854 co.H2P_NH3_NH3P_H2 = CO_findrk("H2+,NH3=>NH3+,H2")*findspecies("NH3")->hevmol;
02855 co.H2P_CN_CNP_H2 = CO_findrk("H2+,CN=>CN+,H2")*findspecies("CN")->hevmol;
02856 co.H2P_HCN_HCNP_H2 = CO_findrk("H2+,HCN=>HCN+,H2")*findspecies("HCN")->hevmol;
02857 co.H2P_NO_NOP_H2 = CO_findrk("H2+,NO=>NO+,H2")*findspecies("NO")->hevmol;
02858
02859
02860 co.HP_C2_C2P_H = CO_findrk("H+,C2=>C2+,H")*findspecies("C2")->hevmol;
02861 co.H2P_C2_C2P_H2 = CO_findrk("H2+,C2=>C2+,H2")*findspecies("C2")->hevmol;
02862 co.Hminus_NH4P_NH3_H2 = CO_findrk("H-,NH4+=>NH3,H2")*findspecies("NH4+")->hevmol;
02863 co.Hminus_NP_N_H = CO_findrk("H-,N+=>N,H")*findspecies("N+")->hevmol;
02864
02865
02866
02867
02868
02869
02870
02871
02872 co.H2_ClP_HClP_H = CO_findrk("H2,Cl+=>HCl+,H")*findspecies("Cl+")->hevmol;
02873 co.H2_HClP_H2ClP_H = CO_findrk("H2,HCl+=>H2Cl+,H")*findspecies("HCl+")->hevmol;
02874 co.H3P_Cl_HClP_H2 = CO_findrk("H3+,Cl=>HCl+,H2")*findspecies("Cl")->hevmol;
02875 co.H3P_HCl_H2ClP_H2 = CO_findrk("H3+,HCl=>H2Cl+,H2")*findspecies("HCl")->hevmol;
02876 co.HP_HCl_HClP_H = CO_findrk("H+,HCl=>HCl+,H")*findspecies("HCl")->hevmol;
02877
02878 co.H2_S_HS_H = CO_findrk("H2,S=>HS,H")*findspecies("S")->hevmol;
02879
02880 co.HP_HNC_HCN_HP = CO_findrk("H+,HNC=>HCN,H+")*findspecies("HNC")->hevmol;
02881 co.H_HNC_HCN_H = CO_findrk("H,HNC=>HCN,H")*findspecies("HNC")->hevmol;
02882
02883
02884 co.H2_HCNP_HCNHP_H = CO_findrk("H2,HCN+=>HCNH+,H")*findspecies("HCN+")->hevmol;
02885
02886 co.H3P_HCN_HCNHP_H2 = CO_findrk("H3+,HCN=>HCNH+,H2")*findspecies("HCN")->hevmol;
02887
02888 co.H2s_OP_OHP_H = CO_findrk("H2*,O+=>OH+,H")*findspecies("O+")->hevmol;
02889
02890
02891
02892
02893
02894 co.HP_C2H2_C2H2P_H = CO_findrk("H+,C2H2=>C2H2+,H")*findspecies("C2H2")->hevmol;
02895 co.HP_C2H2_C2HP_H2 = CO_findrk("H+,C2H2=>C2H+,H2")*findspecies("C2H2")->hevmol;
02896 co.HP_C3H_C3HP_H = CO_findrk("H+,C3H=>C3H+,H")*findspecies("C3H")->hevmol;
02897 co.HP_C3H_C3P_H2 = CO_findrk("H+,C3H=>C3+,H2")*findspecies("C3H")->hevmol;
02898 co.H2P_C2H_C2H2P_H = CO_findrk("H2+,C2H=>C2H2+,H")*findspecies("C2H")->hevmol;
02899 co.H2P_C2H2_C2H2P_H2 = CO_findrk("H2+,C2H2=>C2H2+,H2")*findspecies("C2H2")->hevmol;
02900 co.H3P_C2H_C2H2P_H2 = CO_findrk("H3+,C2H=>C2H2+,H2")*findspecies("C2H")->hevmol;
02901 co.H3P_C3_C3HP_H2 = CO_findrk("H3+,C3=>C3H+,H2")*findspecies("C3")->hevmol;
02902 co.H2_C2HP_C2H2P_H = CO_findrk("H2,C2H+=>C2H2+,H")*findspecies("C2H+")->hevmol;
02903 co.H2_C3P_C3HP_H = CO_findrk("H2,C3+=>C3H+,H")*findspecies("C3+")->hevmol;
02904 co.H_C2H3P_C2H2P_H2 = CO_findrk("H,C2H3+=>C2H2+,H2")*findspecies("C2H3+")->hevmol;
02905 co.H3P_C2H2_C2H3P_H2 = CO_findrk("H3+,C2H2=>C2H3+,H2")*findspecies("C2H2")->hevmol;
02906 co.H2P_C2H2_C2H3P_H = CO_findrk("H2+,C2H2=>C2H3+,H")*findspecies("C2H2")->hevmol;
02907 co.HP_C3_C3P_H = CO_findrk("H+,C3=>C3+,H")*findspecies("C3")->hevmol;
02908 co.HP_C2H_C2HP_H = CO_findrk("H+,C2H=>C2H+,H")*findspecies("C2H")->hevmol;
02909 co.H2P_C2_C2HP_H = CO_findrk("H2+,C2=>C2H+,H")*findspecies("C2")->hevmol;
02910 co.H2P_C2H_C2HP_H2 = CO_findrk("H2+,C2H=>C2H+,H2")*findspecies("C2H")->hevmol;
02911 co.H3P_C2_C2HP_H2 = CO_findrk("H3+,C2=>C2H+,H2")*findspecies("C2")->hevmol;
02912 co.H2_C2P_C2HP_H = CO_findrk("H2,C2+=>C2H+,H")*findspecies("C2+")->hevmol;
02913 co.HP_C2H_C2P_H2 = CO_findrk("H+,C2H=>C2+,H2")*findspecies("C2H")->hevmol;
02914
02915
02916
02917 co.N2_H3P_N2HP_H2 = CO_findrk("N2,H3+=>N2H+,H2")*findspecies("N2")->hevmol;
02918 if(r->next == NULL)
02919 {
02920 int in[]={ipMH},out[]={ipMH,ipMH};
02921 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02922 }
02923 r = r->next;
02924 rindex++;
02925 r->rk = co.H_CH_C_H_H;
02926
02927 if(r->next == NULL)
02928 {
02929 int in[]={ipMH},out[]={ipMH,ipMH};
02930 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02931 }
02932 r = r->next;
02933 rindex++;
02934 r->rk =co.H_OH_O_H_H;
02935
02936 if(r->next == NULL)
02937 {
02938 int in[]={ipMH},out[]={ipMH,ipMH};
02939 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02940 }
02941 r = r->next;
02942 rindex++;
02943 r->rk =co.H_H2O_OH_H_H;
02944
02945 if(r->next == NULL)
02946 {
02947 int in[]={ipMH},out[]={ipMHp};
02948 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02949 }
02950 r = r->next;
02951 rindex++;
02952 r->rk =co.H_COP_CO_HP;
02953
02954 if(r->next == NULL)
02955 {
02956 int in[]={ipMH},out[]={ipMH2g};
02957 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02958 }
02959 r = r->next;
02960 rindex++;
02961 r->rk =co.H_CH_C_H2;
02962
02963 if(r->next == NULL)
02964 {
02965 int in[]={ipMH},out[]={ipMH2g};
02966 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02967 }
02968 r = r->next;
02969 rindex++;
02970 r->rk =co.H_CHP_CP_H2;
02971
02972 if(r->next == NULL)
02973 {
02974 int in[]={ipMH},out[]={ipMH2g};
02975 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02976 }
02977 r = r->next;
02978 rindex++;
02979 r->rk =co.H_CH2_CH_H2;
02980
02981 if(r->next == NULL)
02982 {
02983 int in[]={ipMH},out[]={ipMH2g};
02984 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02985 }
02986 r = r->next;
02987 rindex++;
02988 r->rk =co.H_CH3P_CH2P_H2;
02989
02990 if(r->next == NULL)
02991 {
02992 int in[]={ipMH},out[]={ipMH2g};
02993 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
02994 }
02995 r = r->next;
02996 rindex++;
02997 r->rk =co.H_OH_O_H2;
02998
02999 if(r->next == NULL)
03000 {
03001 int in[]={ipMH},out[]={ipMH2g};
03002 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03003 }
03004 r = r->next;
03005 rindex++;
03006 r->rk =co.H_H2O_OH_H2;
03007
03008
03009 if(r->next == NULL)
03010 {
03011 int in[]={ipMHm},out[]={ipMH2g};
03012 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03013 }
03014 r = r->next;
03015 rindex++;
03016 r->rk =co.Hminus_HCOP_CO_H2;
03017
03018 if(r->next == NULL)
03019 {
03020 int in[]={ipMHm},out[]={ipMH2g};
03021 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03022 }
03023 r = r->next;
03024 rindex++;
03025 r->rk =co.Hminus_H3OP_H2O_H2;
03026
03027 if(r->next == NULL)
03028 {
03029 int in[]={ipMHm},out[]={ipMH2g, ipMH};
03030 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03031 }
03032 r = r->next;
03033 rindex++;
03034 r->rk =co.Hminus_H3OP_OH_H2_H;
03035
03036 if(r->next == NULL)
03037 {
03038 int in[]={ipMHp},out[]={ipMH};
03039 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03040 }
03041 r = r->next;
03042 rindex++;
03043 r->rk =co.HP_CH_CHP_H;
03044
03045 if(r->next == NULL)
03046 {
03047 int in[]={ipMHp},out[]={ipMH};
03048 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03049 }
03050 r = r->next;
03051 rindex++;
03052 r->rk =co.HP_CH2_CH2P_H;
03053
03054 if(r->next == NULL)
03055 {
03056 int in[]={ipMHp},out[]={ipMH};
03057 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03058 }
03059 r = r->next;
03060 rindex++;
03061 r->rk =co.HP_H2O_H2OP_H;
03062
03063 if(r->next == NULL)
03064 {
03065 int in[]={ipMHp},out[]={ipMH};
03066 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03067 }
03068 r = r->next;
03069 rindex++;
03070 r->rk =co.HP_O2_O2P_H;
03071
03072 if(r->next == NULL)
03073 {
03074 int in[]={ipMHp},out[]={ipMH};
03075 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03076 }
03077 r = r->next;
03078 rindex++;
03079 r->rk =co.HP_OH_OHP_H;
03080
03081 if(r->next == NULL)
03082 {
03083 int in[]={ipMHp},out[]={ipMH};
03084 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03085 }
03086 r = r->next;
03087 rindex++;
03088 r->rk =co.HP_SiO_SiOP_H;
03089
03090 if(r->next == NULL)
03091 {
03092 int in[]={ipMHp},out[]={ipMH2g};
03093 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03094 }
03095 r = r->next;
03096 rindex++;
03097 r->rk =co.HP_CH2_CHP_H2;
03098
03099 if(r->next == NULL)
03100 {
03101 int in[]={ipMHp},out[]={ipMH2g};
03102 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03103 }
03104 r = r->next;
03105 rindex++;
03106 r->rk =co.HP_SiH_SiP_H2;
03107
03108 if(r->next == NULL)
03109 {
03110 int in[]={ipMH2g},out[]={ipMH};
03111 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03112 }
03113 r = r->next;
03114 rindex++;
03115 r->rk =co.H2_CHP_CH2P_H;
03116
03117 if(r->next == NULL)
03118 {
03119 int in[]={ipMH2g},out[]={ipMH};
03120 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03121 }
03122 r = r->next;
03123 rindex++;
03124 r->rk =co.H2_CH2P_CH3P_H;
03125
03126 if(r->next == NULL)
03127 {
03128 int in[]={ipMH2g},out[]={ipMH};
03129 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03130 }
03131 r = r->next;
03132 rindex++;
03133 r->rk =co.H2_OHP_H2OP_H;
03134
03135 if(r->next == NULL)
03136 {
03137 int in[]={ipMH2g},out[]={ipMH};
03138 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03139 }
03140 r = r->next;
03141 rindex++;
03142 r->rk =co.H2_H2OP_H3OP_H;
03143
03144 if(r->next == NULL)
03145 {
03146 int in[]={ipMH2g},out[]={ipMH};
03147 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03148 }
03149 r = r->next;
03150 rindex++;
03151 r->rk =co.H2_COP_HCOP_H;
03152
03153 if(r->next == NULL)
03154 {
03155 int in[]={ipMH2g},out[]={ipMH};
03156 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03157 }
03158 r = r->next;
03159 rindex++;
03160 r->rk =co.H2_OP_OHP_H;
03161
03162 if(r->next == NULL)
03163 {
03164 int in[]={ipMH2g},out[]={ipMH};
03165 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03166 }
03167 r = r->next;
03168 rindex++;
03169 r->rk =co.H2_SiOP_SiOHP_H;
03170
03171 if(r->next == NULL)
03172 {
03173 int in[]={ipMH2g},out[]={ipMH};
03174 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03175 }
03176 r = r->next;
03177 rindex++;
03178 r->rk =co.H2_C_CH_H;
03179
03180 if(r->next == NULL)
03181 {
03182 int in[]={ipMH2g},out[]={ipMH};
03183 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03184 }
03185 r = r->next;
03186 rindex++;
03187 r->rk =co.H2_CP_CHP_H;
03188
03189 if(r->next == NULL)
03190 {
03191 int in[]={ipMH2g},out[]={ipMH};
03192 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03193 }
03194 r = r->next;
03195 rindex++;
03196 r->rk =co.H2_CH_CH2_H;
03197
03198 if(r->next == NULL)
03199 {
03200 int in[]={ipMH2g},out[]={ipMH};
03201 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03202 }
03203 r = r->next;
03204 rindex++;
03205 r->rk =co.H2_OH_H2O_H;
03206
03207 if(r->next == NULL)
03208 {
03209 int in[]={ipMH2g},out[]={ipMH};
03210 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03211 }
03212 r = r->next;
03213 rindex++;
03214 r->rk =co.H2_O_OH_H;
03215
03216 if(r->next == NULL)
03217 {
03218 int in[]={ipMH2g},out[]={ipMH2g,ipMH};
03219 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03220 }
03221 r = r->next;
03222 rindex++;
03223 r->rk =co.H2_CH_C_H2_H;
03224
03225 if(r->next == NULL)
03226 {
03227 int in[]={ipMH2g},out[]={ipMH2g, ipMH};
03228 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03229 }
03230 r = r->next;
03231 rindex++;
03232 r->rk =co.H2_OH_O_H2_H;
03233
03234 if(r->next == NULL)
03235 {
03236 int in[]={ipMH2g},out[]={ipMH2g, ipMH};
03237 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03238 }
03239 r = r->next;
03240 rindex++;
03241 r->rk =co.H2_H2O_OH_H2_H;
03242
03243 if(r->next == NULL)
03244 {
03245 int in[]={ipMH2g},out[]={ipMH2g};
03246 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03247 }
03248 r = r->next;
03249 rindex++;
03250 r->rk =co.H2_O2_O_O_H2;
03251
03252
03253
03254
03255
03256
03257
03258
03259
03260
03261
03262
03263 if(r->next == NULL)
03264 {
03265 int in[]={ipMH2s},out[]={ipMH2g,ipMH};
03266 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03267 }
03268 r = r->next;
03269 rindex++;
03270 r->rk =co.H2s_CH_C_H2_H;
03271
03272 if(r->next == NULL)
03273 {
03274 int in[]={ipMH2s},out[]={ipMH2g,ipMH};
03275 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03276 }
03277 r = r->next;
03278 rindex++;
03279 r->rk =co.H2s_OH_O_H2_H;
03280
03281 if(r->next == NULL)
03282 {
03283 int in[]={ipMH2s},out[]={ipMH2g,ipMH};
03284 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03285 }
03286 r = r->next;
03287 rindex++;
03288 r->rk =co.H2s_H2O_OH_H2_H;
03289
03290 if(r->next == NULL)
03291 {
03292 int in[]={ipMH2s},out[]={ipMH2g};
03293 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03294 }
03295 r = r->next;
03296 rindex++;
03297 r->rk =co.H2s_O2_O_O_H2;
03298
03299 if(r->next == NULL)
03300 {
03301 int in[]={ipMH2p},out[]={ipMH};
03302 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03303 }
03304 r = r->next;
03305 rindex++;
03306 r->rk =co.H2P_C_CHP_H;
03307
03308 if(r->next == NULL)
03309 {
03310 int in[]={ipMH2p},out[]={ipMH};
03311 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03312 }
03313 r = r->next;
03314 rindex++;
03315 r->rk =co.H2P_CH_CH2P_H;
03316
03317 if(r->next == NULL) {
03318 int in[]={ipMH2p},out[]={ipMH};
03319 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03320 }
03321 r = r->next;
03322 rindex++;
03323 r->rk =co.H2P_CH2_CH3P_H;
03324
03325 if(r->next == NULL) {
03326 int in[]={ipMH2p},out[]={ipMH};
03327 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03328 }
03329 r = r->next;
03330 rindex++;
03331 r->rk =co.H2P_OH_H2OP_H;
03332
03333 if(r->next == NULL) {
03334 int in[]={ipMH2p},out[]={ipMH};
03335 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03336 }
03337 r = r->next;
03338 rindex++;
03339 r->rk =co.H2P_H2O_H3OP_H;
03340
03341 if(r->next == NULL) {
03342 int in[]={ipMH2p},out[]={ipMH};
03343 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03344 }
03345 r = r->next;
03346 rindex++;
03347 r->rk =co.H2P_CO_HCOP_H;
03348
03349 if(r->next == NULL) {
03350 int in[]={ipMH2p},out[]={ipMH};
03351 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03352 }
03353 r = r->next;
03354 rindex++;
03355 r->rk =co.H2P_O_OHP_H;
03356
03357 if(r->next == NULL) {
03358 int in[]={ipMH2p},out[]={ipMH2g};
03359 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03360 }
03361 r = r->next;
03362 rindex++;
03363 r->rk =co.H2P_CH_CHP_H2;
03364
03365 if(r->next == NULL) {
03366 int in[]={ipMH2p},out[]={ipMH2g};
03367 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03368 }
03369 r = r->next;
03370 rindex++;
03371 r->rk =co.H2P_CH2_CH2P_H2;
03372
03373 if(r->next == NULL) {
03374 int in[]={ipMH2p},out[]={ipMH2g};
03375 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03376 }
03377 r = r->next;
03378 rindex++;
03379 r->rk =co.H2P_CO_COP_H2;
03380
03381 if(r->next == NULL) {
03382 int in[]={ipMH2p},out[]={ipMH2g};
03383 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03384 }
03385 r = r->next;
03386 rindex++;
03387 r->rk =co.H2P_H2O_H2OP_H2;
03388
03389 if(r->next == NULL) {
03390 int in[]={ipMH2p},out[]={ipMH2g};
03391 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03392 }
03393 r = r->next;
03394 rindex++;
03395 r->rk =co.H2P_O2_O2P_H2;
03396
03397 if(r->next == NULL) {
03398 int in[]={ipMH2p},out[]={ipMH2g};
03399 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03400 }
03401 r = r->next;
03402 rindex++;
03403 r->rk =co.H2P_OH_OHP_H2;
03404
03405 if(r->next == NULL) {
03406 int in[]={ipMH3p},out[]={ipMH2g};
03407 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03408 }
03409 r = r->next;
03410 rindex++;
03411 r->rk =co.H3P_C_CHP_H2;
03412
03413 if(r->next == NULL) {
03414 int in[]={ipMH3p},out[]={ipMH2g};
03415 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03416 }
03417 r = r->next;
03418 rindex++;
03419 r->rk =co.H3P_CH_CH2P_H2;
03420
03421 if(r->next == NULL) {
03422 int in[]={ipMH3p},out[]={ipMH2g};
03423 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03424 }
03425 r = r->next;
03426 rindex++;
03427 r->rk =co.H3P_CH2_CH3P_H2;
03428
03429 if(r->next == NULL) {
03430 int in[]={ipMH3p},out[]={ipMH2g};
03431 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03432 }
03433 r = r->next;
03434 rindex++;
03435 r->rk =co.H3P_OH_H2OP_H2;
03436
03437 if(r->next == NULL) {
03438 int in[]={ipMH3p},out[]={ipMH2g};
03439 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03440 }
03441 r = r->next;
03442 rindex++;
03443 r->rk =co.H3P_H2O_H3OP_H2;
03444
03445 if(r->next == NULL) {
03446 int in[]={ipMH3p},out[]={ipMH2g};
03447 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03448 }
03449 r = r->next;
03450 rindex++;
03451 r->rk =co.H3P_CO_HCOP_H2;
03452
03453 if(r->next == NULL)
03454 {
03455 int in[]={ipMH3p},out[]={ipMH2g};
03456 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03457 }
03458 r = r->next;
03459 rindex++;
03460 r->rk =co.H3P_O_OHP_H2;
03461
03462 if(r->next == NULL)
03463 {
03464 int in[]={ipMH3p},out[]={ipMH2g};
03465 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03466 }
03467 r = r->next;
03468 rindex++;
03469 r->rk =co.H3P_SiH_SiH2P_H2;
03470
03471 if(r->next == NULL) {
03472 int in[]={ipMH3p},out[]={ipMH2g};
03473 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03474 }
03475 r = r->next;
03476 rindex++;
03477 r->rk =co.H3P_SiO_SiOHP_H2;
03478
03479 if(r->next == NULL)
03480 {
03481 int in[]={ipMH2s},out[]={ipMH};
03482 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03483 }
03484 r = r->next;
03485 rindex++;
03486 r->rk =co.H2s_CH_CH2_H;
03487
03488 if(r->next == NULL)
03489 {
03490 int in[]={ipMH2s},out[]={ipMH};
03491 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03492 }
03493 r = r->next;
03494 rindex++;
03495 r->rk =co.H2s_O_OH_H;
03496
03497 if(r->next == NULL)
03498 {
03499 int in[]={ipMH2s},out[]={ipMH};
03500 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03501 }
03502 r = r->next;
03503 rindex++;
03504 r->rk =co.H2s_OH_H2O_H;
03505
03506
03507 if(r->next == NULL)
03508 {
03509 int in[]={ipMH2s},out[]={ipMH};
03510 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03511 }
03512 r = r->next;
03513 rindex++;
03514 r->rk =co.H2s_C_CH_H;
03515
03516 if(r->next == NULL)
03517 {
03518 int in[]={ipMH2s},out[]={ipMH};
03519 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03520 }
03521 r = r->next;
03522 rindex++;
03523 r->rk =co.H2s_CP_CHP_H;
03524
03525 if(r->next == NULL)
03526 {
03527 int in[]={ipMH},out[]={ipMH2g};
03528 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03529 }
03530 r = r->next;
03531 rindex++;
03532 r->rk =co.H_CH3_CH2_H2;
03533
03534 if(r->next == NULL)
03535 {
03536 int in[]={ipMH},out[]={ipMH2g};
03537 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03538 }
03539 r = r->next;
03540 rindex++;
03541 r->rk =co.H_CH4P_CH3P_H2;
03542
03543 if(r->next == NULL)
03544 {
03545 int in[]={ipMH},out[]={ipMH2g};
03546 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03547 }
03548 r = r->next;
03549 rindex++;
03550 r->rk =co.H_CH5P_CH4P_H2;
03551
03552 if(r->next == NULL)
03553 {
03554 int in[]={ipMH2g},out[]={ipMH};
03555 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03556 }
03557 r = r->next;
03558 rindex++;
03559 r->rk =co.H2_CH2_CH3_H;
03560
03561
03562 if(r->next == NULL)
03563 {
03564 int in[]={ipMH2g},out[]={ipMH};
03565 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03566 }
03567 r = r->next;
03568 rindex++;
03569 r->rk = co.H2_CH3_CH4_H;
03570
03571 if(r->next == NULL)
03572 {
03573 int in[]={ipMH2g},out[]={ipMH};
03574 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03575 }
03576 r = r->next;
03577 rindex++;
03578 r->rk = co.H2_CH4P_CH5P_H;
03579
03580 if(r->next == NULL)
03581 {
03582 int in[]={ipMH2s},out[]={ipMH};
03583 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03584 }
03585 r = r->next;
03586 rindex++;
03587 r->rk =co.H2s_CH2_CH3_H;
03588
03589 if(r->next == NULL)
03590 {
03591 int in[]={ipMH2s},out[]={ipMH};
03592 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03593 }
03594 r = r->next;
03595 rindex++;
03596 r->rk = co.H2s_CH3_CH4_H;
03597
03598 if(r->next == NULL)
03599 {
03600 int in[]={ipMH2p},out[]={ipMH2g};
03601 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03602 }
03603 r = r->next;
03604 rindex++;
03605 r->rk = co.H2P_CH4_CH3P_H2;
03606
03607 if(r->next == NULL)
03608 {
03609 int in[]={ipMH2p},out[]={ipMH2g};
03610 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03611 }
03612 r = r->next;
03613 rindex++;
03614 r->rk = co.H2P_CH4_CH4P_H2;
03615
03616 if(r->next == NULL)
03617 {
03618 int in[]={ipMH2p},out[]={ipMH};
03619 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03620 }
03621 r = r->next;
03622 rindex++;
03623 r->rk = co.H2P_CH4_CH5P_H;
03624
03625 if(r->next == NULL)
03626 {
03627 int in[]={ipMH3p},out[]={ipMH2g};
03628 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03629 }
03630 r = r->next;
03631 rindex++;
03632 r->rk = co.H3P_CH3_CH4P_H2;
03633
03634 if(r->next == NULL)
03635 {
03636 int in[]={ipMH3p},out[]={ipMH2g};
03637 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03638 }
03639 r = r->next;
03640 rindex++;
03641 r->rk = co.H3P_CH4_CH5P_H2;
03642
03643 if(r->next == NULL) {
03644 int in[]={ipMHp},out[]={ipMH2g};
03645 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03646 }
03647 r = r->next;
03648 rindex++;
03649 r->rk = co.HP_CH4_CH3P_H2;
03650
03651 if(r->next == NULL) {
03652 int in[]={ipMHp},out[]={ipMH};
03653 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03654 }
03655 r = r->next;
03656 rindex++;
03657 r->rk = co.HP_CH4_CH4P_H;
03658
03659 if(r->next == NULL) {
03660 int in[]={ipMH2g},out[]={ipMH};
03661 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03662 }
03663 r = r->next;
03664 rindex++;
03665 r->rk = co.H2_ClP_HClP_H;
03666
03667 if(r->next == NULL) {
03668 int in[]={ipMH2g},out[]={ipMH};
03669 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03670 }
03671 r = r->next;
03672 rindex++;
03673 r->rk = co.H2_HClP_H2ClP_H;
03674
03675 if(r->next == NULL) {
03676 int in[]={ipMH3p},out[]={ipMH2g};
03677 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03678 }
03679 r = r->next;
03680 rindex++;
03681 r->rk = co.H3P_Cl_HClP_H2;
03682
03683 if(r->next == NULL) {
03684 int in[]={ipMH3p},out[]={ipMH2g};
03685 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03686 }
03687 r = r->next;
03688 rindex++;
03689 r->rk = co.H3P_HCl_H2ClP_H2;
03690
03691 if(r->next == NULL) {
03692 int in[]={ipMHp},out[]={ipMH};
03693 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03694 }
03695 r = r->next;
03696 rindex++;
03697 r->rk = co.HP_HCl_HClP_H;
03698
03699
03700 if(r->next == NULL) {
03701 int in[]={ipMH2g},out[]={ipMH};
03702 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03703 }
03704 r = r->next;
03705 rindex++;
03706 r->rk = co.H2_S_HS_H;
03707
03708
03709
03710
03711
03712 if(r->next == NULL)
03713 {
03714 int in[]={ipMH2g},out[]={ipMH};
03715 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03716 }
03717 r = r->next;
03718 rindex++;
03719 r->rk = co.H2_N_NH_H;
03720
03721 if(r->next == NULL)
03722 {
03723 int in[]={ipMH2g},out[]={ipMH};
03724 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03725 }
03726 r = r->next;
03727 rindex++;
03728 r->rk = co.H2_NH_NH2_H;
03729
03730
03731 if(r->next == NULL)
03732 {
03733 int in[]={ipMH2g},out[]={ipMH};
03734 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03735 }
03736 r = r->next;
03737 rindex++;
03738 r->rk = co.H2_NH2_NH3_H;
03739
03740 if(r->next == NULL)
03741 {
03742 int in[]={ipMH2g},out[]={ipMH};
03743 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03744 }
03745 r = r->next;
03746 rindex++;
03747 r->rk = co.H2_CN_HCN_H;
03748
03749 if(r->next == NULL)
03750 {
03751 int in[]={ipMHp},out[]={ipMH2g};
03752 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03753 }
03754 r = r->next;
03755 rindex++;
03756 r->rk = co.HP_HNO_NOP_H2;
03757
03758 if(r->next == NULL)
03759 {
03760 int in[]={ipMHp},out[]={ipMH2g};
03761 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03762 }
03763 r = r->next;
03764 rindex++;
03765 r->rk = co.HP_HS_SP_H2;
03766
03767 if(r->next == NULL)
03768 {
03769 int in[]={ipMH},out[]={ipMH2g};
03770 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03771 }
03772 r = r->next;
03773 rindex++;
03774 r->rk = co.H_HSP_SP_H2;
03775
03776 if(r->next == NULL)
03777 {
03778 int in[]={ipMH2p},out[]={ipMH};
03779 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03780 }
03781 r = r->next;
03782 rindex++;
03783 r->rk = co.H2P_N_NHP_H;
03784
03785 if(r->next == NULL)
03786 {
03787 int in[]={ipMH2g},out[]={ipMH};
03788 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03789 }
03790 r = r->next;
03791 rindex++;
03792 r->rk = co.H2_NP_NHP_H;
03793
03794 if(r->next == NULL)
03795 {
03796 int in[]={ipMH2g},out[]={ipMH3p};
03797 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03798 }
03799 r = r->next;
03800 rindex++;
03801 r->rk = co.H2_NHP_N_H3P;
03802
03803 if(r->next == NULL)
03804 {
03805 int in[]={ipMH2p},out[]={ipMH};
03806 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03807 }
03808 r = r->next;
03809 rindex++;
03810 r->rk = co.H2P_NH_NH2P_H;
03811
03812
03813 if(r->next == NULL)
03814 {
03815 int in[]={ipMH2g},out[]={ipMH};
03816 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03817 }
03818 r = r->next;
03819 rindex++;
03820 r->rk = co.H2_NHP_NH2P_H;
03821
03822 if(r->next == NULL)
03823 {
03824 int in[]={ipMH2g},out[]={ipMH};
03825 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03826 }
03827 r = r->next;
03828 rindex++;
03829 r->rk = co.H2_NH2P_NH3P_H;
03830
03831 if(r->next == NULL)
03832 {
03833 int in[]={ipMH2g},out[]={ipMH};
03834 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03835 }
03836 r = r->next;
03837 rindex++;
03838 r->rk = co.H2_NH3P_NH4P_H;
03839
03840 if(r->next == NULL)
03841 {
03842 int in[]={ipMH2p},out[]={ipMH};
03843 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03844 }
03845 r = r->next;
03846 rindex++;
03847 r->rk = co.H2P_CN_HCNP_H;
03848
03849 if(r->next == NULL)
03850 {
03851 int in[]={ipMH2g},out[]={ipMH};
03852 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03853 }
03854 r = r->next;
03855 rindex++;
03856 r->rk = co.H2_CNP_HCNP_H;
03857
03858 if(r->next == NULL)
03859 {
03860 int in[]={ipMH2p},out[]={ipMH};
03861 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03862 }
03863 r = r->next;
03864 rindex++;
03865 r->rk = co.H2P_NO_HNOP_H;
03866
03867 if(r->next == NULL)
03868 {
03869 int in[]={ipMH2g},out[]={ipMH};
03870 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03871 }
03872 r = r->next;
03873 rindex++;
03874 r->rk = co.H2_SP_HSP_H;
03875
03876 if(r->next == NULL)
03877 {
03878 int in[]={ipMH2g},out[]={ipMH};
03879 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03880 }
03881 r = r->next;
03882 rindex++;
03883 r->rk = co.H2_CSP_HCSP_H;
03884
03885
03886 if(r->next == NULL)
03887 {
03888 int in[]={ipMH3p},out[]={ipMH2g};
03889 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03890 }
03891 r = r->next;
03892 rindex++;
03893 r->rk = co.H3P_NH_NH2P_H2;
03894
03895 if(r->next == NULL)
03896 {
03897 int in[]={ipMH3p},out[]={ipMH2g};
03898 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03899 }
03900 r = r->next;
03901 rindex++;
03902 r->rk = co.H3P_NH2_NH3P_H2;
03903
03904 if(r->next == NULL)
03905 {
03906 int in[]={ipMH3p},out[]={ipMH2g};
03907 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03908 }
03909 r = r->next;
03910 rindex++;
03911 r->rk = co.H3P_NH3_NH4P_H2;
03912
03913 if(r->next == NULL)
03914 {
03915 int in[]={ipMH3p},out[]={ipMH2g};
03916 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03917 }
03918 r = r->next;
03919 rindex++;
03920 r->rk = co.H3P_CN_HCNP_H2;
03921
03922 if(r->next == NULL)
03923 {
03924 int in[]={ipMH3p},out[]={ipMH2g};
03925 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03926 }
03927 r = r->next;
03928 rindex++;
03929 r->rk = co.H3P_NO_HNOP_H2;
03930
03931 if(r->next == NULL)
03932 {
03933 int in[]={ipMH3p},out[]={ipMH2g};
03934 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03935 }
03936 r = r->next;
03937 rindex++;
03938 r->rk = co.H3P_S_HSP_H2;
03939
03940 if(r->next == NULL)
03941 {
03942 int in[]={ipMH3p},out[]={ipMH2g};
03943 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03944 }
03945 r = r->next;
03946 rindex++;
03947 r->rk = co.H3P_CS_HCSP_H2;
03948
03949 if(r->next == NULL)
03950 {
03951 int in[]={ipMH3p},out[]={ipMH2g};
03952 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03953 }
03954 r = r->next;
03955 rindex++;
03956 r->rk = co.H3P_NO2_NOP_OH_H2;
03957
03958 if(r->next == NULL)
03959 {
03960 int in[]={ipMHp},out[]={ipMH};
03961 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03962 }
03963 r = r->next;
03964 rindex++;
03965 r->rk = co.HP_NH_NHP_H;
03966
03967 if(r->next == NULL)
03968 {
03969 int in[]={ipMHp},out[]={ipMH};
03970 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03971 }
03972 r = r->next;
03973 rindex++;
03974 r->rk = co.HP_NH2_NH2P_H;
03975
03976 if(r->next == NULL)
03977 {
03978 int in[]={ipMHp},out[]={ipMH};
03979 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03980 }
03981 r = r->next;
03982 rindex++;
03983 r->rk = co.HP_NH3_NH3P_H;
03984
03985
03986 if(r->next == NULL)
03987 {
03988 int in[]={ipMH},out[]={ipMHp};
03989 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03990 }
03991 r = r->next;
03992 rindex++;
03993 r->rk = co.H_CNP_CN_HP;
03994
03995 if(r->next == NULL)
03996 {
03997 int in[]={ipMHp},out[]={ipMH};
03998 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
03999 }
04000 r = r->next;
04001 rindex++;
04002 r->rk = co.HP_HCN_HCNP_H;
04003
04004 if(r->next == NULL)
04005 {
04006 int in[]={ipMH},out[]={ipMHp};
04007 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04008 }
04009 r = r->next;
04010 rindex++;
04011 r->rk = co.H_HCNP_HCN_HP;
04012
04013 if(r->next == NULL)
04014 {
04015 int in[]={ipMH},out[]={ipMHp};
04016 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04017 }
04018 r = r->next;
04019 rindex++;
04020 r->rk = co.H_N2P_N2_HP;
04021
04022 if(r->next == NULL)
04023 {
04024 int in[]={ipMHp},out[]={ipMH};
04025 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04026 }
04027 r = r->next;
04028 rindex++;
04029 r->rk = co.HP_NO_NOP_H;
04030
04031 if(r->next == NULL)
04032 {
04033 int in[]={ipMHp},out[]={ipMH};
04034 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04035 }
04036 r = r->next;
04037 rindex++;
04038 r->rk = co.HP_HS_HSP_H;
04039
04040 if(r->next == NULL)
04041 {
04042 int in[]={ipMHp},out[]={ipMH};
04043 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04044 }
04045 r = r->next;
04046 rindex++;
04047 r->rk = co.HP_SiN_SiNP_H;
04048
04049 if(r->next == NULL)
04050 {
04051 int in[]={ipMHp},out[]={ipMH};
04052 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04053 }
04054 r = r->next;
04055 rindex++;
04056 r->rk = co.HP_CS_CSP_H;
04057
04058 if(r->next == NULL)
04059 {
04060 int in[]={ipMHp},out[]={ipMH};
04061 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04062 }
04063 r = r->next;
04064 rindex++;
04065 r->rk = co.HP_NS_NSP_H;
04066
04067 if(r->next == NULL)
04068 {
04069 int in[]={ipMHp},out[]={ipMH};
04070 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04071 }
04072 r = r->next;
04073 rindex++;
04074 r->rk = co.HP_SO_SOP_H;
04075
04076
04077 if(r->next == NULL)
04078 {
04079 int in[]={ipMHp},out[]={ipMH};
04080 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04081 }
04082 r = r->next;
04083 rindex++;
04084 r->rk = co.HP_OCS_OCSP_H;
04085
04086 if(r->next == NULL)
04087 {
04088 int in[]={ipMHp},out[]={ipMH};
04089 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04090 }
04091 r = r->next;
04092 rindex++;
04093 r->rk = co.HP_S2_S2P_H;
04094
04095 if(r->next == NULL)
04096 {
04097 int in[]={ipMH2p},out[]={ipMH2g};
04098 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04099 }
04100 r = r->next;
04101 rindex++;
04102 r->rk = co.H2P_NH_NHP_H2;
04103
04104 if(r->next == NULL)
04105 {
04106 int in[]={ipMH2p},out[]={ipMH2g};
04107 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04108 }
04109 r = r->next;
04110 rindex++;
04111 r->rk = co.H2P_NH2_NH2P_H2;
04112
04113 if(r->next == NULL)
04114 {
04115 int in[]={ipMH2p},out[]={ipMH2g};
04116 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04117 }
04118 r = r->next;
04119 rindex++;
04120 r->rk = co.H2P_NH3_NH3P_H2;
04121
04122 if(r->next == NULL)
04123 {
04124 int in[]={ipMH2p},out[]={ipMH2g};
04125 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04126 }
04127 r = r->next;
04128 rindex++;
04129 r->rk = co.H2P_CN_CNP_H2;
04130
04131 if(r->next == NULL)
04132 {
04133 int in[]={ipMH2p},out[]={ipMH2g};
04134 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04135 }
04136 r = r->next;
04137 rindex++;
04138 r->rk = co.H2P_HCN_HCNP_H2;
04139
04140 if(r->next == NULL)
04141 {
04142 int in[]={ipMH2p},out[]={ipMH2g};
04143 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04144 }
04145 r = r->next;
04146 rindex++;
04147 r->rk = co.H2P_NO_NOP_H2;
04148
04149 if(r->next == NULL)
04150 {
04151 int in[]={ipMHp},out[]={ipMH};
04152 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04153 }
04154 r = r->next;
04155 rindex++;
04156 r->rk = co.HP_C2_C2P_H;
04157
04158
04159 if(r->next == NULL)
04160 {
04161 int in[]={ipMH2p},out[]={ipMH2g};
04162 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04163 }
04164 r = r->next;
04165 rindex++;
04166 r->rk = co.H2P_C2_C2P_H2;
04167
04168 if(r->next == NULL)
04169 {
04170 int in[]={ipMHm},out[]={ipMH2g};
04171 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04172 }
04173 r = r->next;
04174 rindex++;
04175 r->rk = co.Hminus_NH4P_NH3_H2;
04176
04177 if(r->next == NULL)
04178 {
04179 int in[]={ipMHm},out[]={ipMH};
04180 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04181 }
04182 r = r->next;
04183 rindex++;
04184 r->rk = co.Hminus_NP_N_H;
04185
04186
04187
04188 if(r->next == NULL)
04189 {
04190 int in[]={ipMHp},out[]={ipMHp};
04191 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04192 }
04193 r = r->next;
04194 rindex++;
04195 r->rk = co.HP_HNC_HCN_HP;
04196
04197 if(r->next == NULL)
04198 {
04199 int in[]={ipMH},out[]={ipMH};
04200 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04201 }
04202 r = r->next;
04203 rindex++;
04204 r->rk = co.H_HNC_HCN_H;
04205
04206
04207
04208
04209
04210
04211
04212
04213
04214
04215
04216 if(r->next == NULL)
04217 {
04218 int in[]={ipMH2g},out[]={ipMH};
04219 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04220 }
04221 r = r->next;
04222 rindex++;
04223 r->rk = co.H2_HCNP_HCNHP_H;
04224
04225 if(r->next == NULL)
04226 {
04227 int in[]={ipMH3p},out[]={ipMH2g};
04228 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04229 }
04230 r = r->next;
04231 rindex++;
04232 r->rk = co.H3P_HCN_HCNHP_H2;
04233
04234
04235 if(r->next == NULL)
04236 {
04237 int in[]={ipMH2s},out[]={ipMH};
04238 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04239 }
04240 r = r->next;
04241 rindex++;
04242 r->rk =co.H2s_OP_OHP_H;
04243
04244
04245
04246 if(r->next == NULL)
04247 {
04248 int in[]={ipMHp},out[]={ipMH};
04249 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04250 }
04251 r = r->next;
04252 rindex++;
04253 r->rk = co.HP_C2H2_C2H2P_H;
04254
04255 if(r->next == NULL)
04256 {
04257 int in[]={ipMHp},out[]={ipMH2g};
04258 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04259 }
04260 r = r->next;
04261 rindex++;
04262 r->rk = co.HP_C2H2_C2HP_H2;
04263
04264 if(r->next == NULL)
04265 {
04266 int in[]={ipMHp},out[]={ipMH};
04267 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04268 }
04269 r = r->next;
04270 rindex++;
04271 r->rk = co.HP_C3H_C3HP_H;
04272
04273 if(r->next == NULL)
04274 {
04275 int in[]={ipMHp},out[]={ipMH2g};
04276 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04277 }
04278 r = r->next;
04279 rindex++;
04280 r->rk =co.HP_C3H_C3P_H2;
04281
04282 if(r->next == NULL)
04283 {
04284 int in[]={ipMH2p},out[]={ipMH};
04285 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04286 }
04287 r = r->next;
04288 rindex++;
04289 r->rk = co.H2P_C2H_C2H2P_H;
04290
04291 if(r->next == NULL)
04292 {
04293 int in[]={ipMH2p},out[]={ipMH2g};
04294 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04295 }
04296 r = r->next;
04297 rindex++;
04298 r->rk = co.H2P_C2H2_C2H2P_H2;
04299
04300 if(r->next == NULL)
04301 {
04302 int in[]={ipMH3p},out[]={ipMH2g};
04303 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04304 }
04305 r = r->next;
04306 rindex++;
04307 r->rk = co.H3P_C2H_C2H2P_H2;
04308
04309 if(r->next == NULL)
04310 {
04311 int in[]={ipMH3p},out[]={ipMH2g};
04312 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04313 }
04314 r = r->next;
04315 rindex++;
04316 r->rk = co.H3P_C3_C3HP_H2;
04317
04318 if(r->next == NULL)
04319 {
04320 int in[]={ipMH2g},out[]={ipMH};
04321 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04322 }
04323 r = r->next;
04324 rindex++;
04325 r->rk = co.H2_C2HP_C2H2P_H;
04326
04327 if(r->next == NULL)
04328 {
04329 int in[]={ipMH2g},out[]={ipMH};
04330 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04331 }
04332 r = r->next;
04333 rindex++;
04334 r->rk = co.H2_C3P_C3HP_H;
04335
04336 if(r->next == NULL)
04337 {
04338 int in[]={ipMH},out[]={ipMH2g};
04339 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04340 }
04341 r = r->next;
04342 rindex++;
04343 r->rk = co.H_C2H3P_C2H2P_H2;
04344
04345 if(r->next == NULL)
04346 {
04347 int in[]={ipMH3p},out[]={ipMH2g};
04348 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04349 }
04350 r = r->next;
04351 rindex++;
04352 r->rk = co.H3P_C2H2_C2H3P_H2;
04353
04354 if(r->next == NULL)
04355 {
04356 int in[]={ipMH2p},out[]={ipMH};
04357 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04358 }
04359 r = r->next;
04360 rindex++;
04361 r->rk = co.H2P_C2H2_C2H3P_H;
04362
04363 if(r->next == NULL)
04364 {
04365 int in[]={ipMHp},out[]={ipMH};
04366 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04367 }
04368 r = r->next;
04369 rindex++;
04370 r->rk = co.HP_C3_C3P_H;
04371
04372 if(r->next == NULL)
04373 {
04374 int in[]={ipMH2p},out[]={ipMH};
04375 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04376 }
04377 r = r->next;
04378 rindex++;
04379 r->rk = co.H2P_C2_C2HP_H;
04380
04381 if(r->next == NULL)
04382 {
04383 int in[]={ipMH2p},out[]={ipMH2g};
04384 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04385 }
04386 r = r->next;
04387 rindex++;
04388 r->rk = co.H2P_C2H_C2HP_H2;
04389
04390 if(r->next == NULL)
04391 {
04392 int in[]={ipMH3p},out[]={ipMH2g};
04393 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04394 }
04395 r = r->next;
04396 rindex++;
04397 r->rk = co.H3P_C2_C2HP_H2;
04398
04399 if(r->next == NULL)
04400 {
04401 int in[]={ipMH2g},out[]={ipMH};
04402 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04403 }
04404 r = r->next;
04405 rindex++;
04406 r->rk = co.H2_C2P_C2HP_H;
04407
04408 if(r->next == NULL)
04409 {
04410 int in[]={ipMHp},out[]={ipMH2g};
04411 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04412 }
04413 r = r->next;
04414 rindex++;
04415 r->rk = co.HP_C2H_C2P_H2;
04416
04417
04418
04419
04420
04421
04422 if(r->next == NULL)
04423 {
04424 int in[]={ipMH3p},out[]={ipMH2g};
04425 r->next = newreaction(rindex,in,INTSZ(in),out,INTSZ(out),NULL,0);
04426 }
04427 r = r->next;
04428 rindex++;
04429 r->rk = co.N2_H3P_N2HP_H2;
04430 }
04431 else
04432 {
04433 fixit();
04434 }
04435
04436
04437
04438
04439 # if 0
04440
04441
04442
04443
04444
04445
04446
04447
04448 for(i=0; i<N_H_MOLEC; ++i)
04449 {
04450 co.hydro_source[i] = 0;
04451 co.hydro_sink[i] = 0;
04452 }
04453
04454
04455
04456 co.C_H3OP_HCOP_H2_1 = HMRATE(1.0e-11,0,0)*findspecies("H3O+")->hevmol*findspecies("C")->hevmol;
04457 co.hydro_source[ipMH2g] += (realnum)co.C_H3OP_HCOP_H2_1;
04458
04459
04460
04461
04462 co.C_OH_CO_H_1 = HMRATE(1.1e-10,0.5,0)*findspecies("OH")->hevmol*findspecies("C")->hevmol;
04463 co.hydro_source[ipMH] += (realnum)co.C_OH_CO_H_1;
04464
04465
04466
04467 co.CP_OH_CO_HP_1 = HMRATE(7.7e-10,0,0)*findspecies("OH")->hevmol*findspecies("C+")->hevmol;
04468 co.hydro_source[ipMHp] += (realnum)co.CP_OH_CO_HP_1;
04469
04470
04471
04472 co.CP_H2O_HCOP_H_1 = HMRATE(9.0e-10,0,0)*findspecies("H2O")->hevmol*findspecies("C+")->hevmol;
04473 co.hydro_source[ipMH] += (realnum)co.CP_H2O_HCOP_H_1;
04474
04475
04476
04477 co.CP_OH_COP_H_1 = HMRATE(7.7e-10,0,0)*findspecies("OH")->hevmol*findspecies("C+")->hevmol;
04478 co.hydro_source[ipMH] += (realnum)co.CP_OH_COP_H_1;
04479
04480
04481
04482 co.O_CH_CO_H_1 = HMRATE(6.6e-11,0,0)*findspecies("CH")->hevmol*findspecies("O")->hevmol;
04483 co.hydro_source[ipMH] += (realnum)co.O_CH_CO_H_1;
04484
04485
04486
04487 co.O_CHP_COP_H_1 = HMRATE(3.5e-10,0,0)*findspecies("CH+")->hevmol*findspecies("O")->hevmol;
04488 co.hydro_source[ipMH] += (realnum)co.O_CHP_COP_H_1;
04489
04490
04491
04492 co.O_CH2_CO_H_H_1 = HMRATE(1.33e-10,0,0)*findspecies("CH2")->hevmol*findspecies("O")->hevmol;
04493 co.hydro_source[ipMH] += (realnum)co.O_CH2_CO_H_H_1;
04494
04495
04496
04497 co.O_CH2_CO_H2_1 = HMRATE(8.0e-11,0,0)*findspecies("CH2")->hevmol*findspecies("O")->hevmol;
04498 co.hydro_source[ipMH2g] += (realnum)co.O_CH2_CO_H2_1;
04499
04500
04501
04502 co.O_CH2P_HCOP_H_1 = HMRATE(7.5e-10,0,0)*findspecies("CH2+")->hevmol*findspecies("O")->hevmol;
04503 co.hydro_source[ipMH] += (realnum)co.O_CH2P_HCOP_H_1;
04504
04505
04506
04507 co.O_CH3P_HCOP_H2_1 = HMRATE(4.0e-10,0,0)*findspecies("CH3+")->hevmol*findspecies("O")->hevmol;
04508 co.hydro_source[ipMH2g] += (realnum)co.O_CH3P_HCOP_H2_1;
04509
04510
04511
04512 co.O_H2OP_O2P_H2_1 = HMRATE(4.0e-11,0,0)*findspecies("H2O+")->hevmol*findspecies("O")->hevmol;
04513 co.hydro_source[ipMH2g] += (realnum)co.O_H2OP_O2P_H2_1;
04514
04515
04516
04517 co.O_OH_O2_H_1 = HMRATE(4.34e-11,-0.5,30)*findspecies("OH")->hevmol*findspecies("O")->hevmol;
04518 co.hydro_source[ipMH] += (realnum)co.O_OH_O2_H_1;
04519
04520
04521
04522 co.O_OHP_O2P_H_1 = HMRATE(7.1e-10,0,0)*findspecies("OH+")->hevmol*findspecies("O")->hevmol;
04523 co.hydro_source[ipMH] += (realnum)co.O_OHP_O2P_H_1;
04524
04525
04526
04527 co.O_SiH_SiO_H_1 = HMRATE(4.0e-11,0.5,0)*findspecies("SiH")->hevmol*findspecies("O")->hevmol;
04528 co.hydro_source[ipMH] += (realnum)co.O_SiH_SiO_H_1;
04529
04530
04531
04532 co.O_SiH2P_SiOHP_H_1 = HMRATE(6.3e-10,0,0)*findspecies("SiH2+")->hevmol*findspecies("O")->hevmol;
04533 co.hydro_source[ipMH] += (realnum)co.O_SiH2P_SiOHP_H_1;
04534
04535
04536
04537 co.OP_CH_COP_H_1 = HMRATE(3.5e-10,0,0)*findspecies("CH")->hevmol*findspecies("O+")->hevmol;
04538 co.hydro_source[ipMH] += (realnum)co.OP_CH_COP_H_1;
04539
04540
04541
04542 co.OP_OH_O2P_H_1 = HMRATE(3.6e-10,0,0)*findspecies("OH")->hevmol*findspecies("O+")->hevmol;
04543 co.hydro_source[ipMH] += (realnum)co.OP_OH_O2P_H_1;
04544
04545
04546
04547 co.Si_OH_SiO_H_1 = HMRATE(2.0e-10,0.5,0)*findspecies("OH")->hevmol*findspecies("Si")->hevmol;
04548 co.hydro_source[ipMH] += (realnum)co.Si_OH_SiO_H_1;
04549
04550
04551
04552 co.SiP_H2O_SiOHP_H_1 = HMRATE(2.3e-10,0,0)*findspecies("H2O")->hevmol*findspecies("Si+")->hevmol;
04553 co.hydro_source[ipMH] += (realnum)co.SiP_H2O_SiOHP_H_1;
04554
04555
04556
04557 co.SiP_OH_SiOP_H_1 = HMRATE(6.3e-10,0,0)*findspecies("OH")->hevmol*findspecies("Si+")->hevmol;
04558 co.hydro_source[ipMH] += (realnum)co.SiP_OH_SiOP_H_1;
04559
04560
04561
04562 co.CHP_H2O_HCOP_H2_1 = HMRATE(2.9e-9,0,0)*findspecies("H2O")->hevmol*findspecies("CH+")->hevmol;
04563 co.hydro_source[ipMH2g] += (realnum)co.CHP_H2O_HCOP_H2_1;
04564
04565
04566
04567 co.CHP_OH_COP_H2_1 = HMRATE(7.5e-10,0,0)*findspecies("OH")->hevmol*findspecies("CH+")->hevmol;
04568 co.hydro_source[ipMH2g] += (realnum)co.CHP_OH_COP_H2_1;
04569
04570
04571
04572 co.H_C_CH_nu = HMRATE(0.00000000000000001,0,0)*dense.xIonDense[ipHYDROGEN][0]*findspecies("C")->hevmol;
04573 co.hydro_sink[ipMH] += (realnum)co.H_C_CH_nu;
04574
04575
04576
04577 co.H_CP_CHP_nu = HMRATE(1.7e-17,0,0)*findspecies("C+")->hevmol*dense.xIonDense[ipHYDROGEN][0];
04578 co.hydro_sink[ipMH] += (realnum)co.H_CP_CHP_nu;
04579
04580
04581
04582 co.H_OH_H2O_nu = HMRATE(5.26E-18,-5.22,90)*findspecies("OH")->hevmol*dense.xIonDense[ipHYDROGEN][0];
04583 co.hydro_sink[ipMH] += (realnum)co.H_OH_H2O_nu;
04584
04585
04586
04587 co.Hminus_CH_CH2_e = HMRATE(1.0e-10,0,0)*findspecies("CH")->hevmol*hmi.Hmolec[ipMHm];
04588 co.hydro_sink[ipMHm] += (realnum)co.Hminus_CH_CH2_e;
04589
04590
04591
04592 co.Hminus_C_CH_e = HMRATE(1.0e-9,0,0)*findspecies("C")->hevmol*hmi.Hmolec[ipMHm];
04593 co.hydro_sink[ipMHm] += (realnum)co.Hminus_C_CH_e;
04594
04595
04596
04597 co.Hminus_OH_H2O_e = HMRATE(1.0e-10,0,0)*findspecies("OH")->hevmol*hmi.Hmolec[ipMHm];
04598 co.hydro_sink[ipMHm] += (realnum)co.Hminus_OH_H2O_e;
04599
04600
04601
04602 co.Hminus_O_OH_e = HMRATE(1.0e-9,0,0)*findspecies("O")->hevmol*hmi.Hmolec[ipMHm];
04603 co.hydro_sink[ipMHm] += (realnum)co.Hminus_O_OH_e;
04604
04605
04606
04607 co.H2_C_CH2_nu = HMRATE(1.0e-17,0,0)*findspecies("C")->hevmol*hmi.Hmolec[ipMH2g];
04608 co.hydro_sink[ipMH2g] += (realnum)co.H2_C_CH2_nu;
04609
04610
04611
04612 co.H2_CP_CH2P_nu = HMRATE(4.0e-16,-0.2,0)*findspecies("C+")->hevmol*hmi.Hmolec[ipMH2g];
04613 co.hydro_sink[ipMH2g] += (realnum)co.H2_CP_CH2P_nu;
04614
04615
04616
04617 co.H2_SiP_SiH2P_nu = HMRATE(3.0e-18,0,0)*findspecies("Si+")->hevmol*hmi.Hmolec[ipMH2g];
04618 co.hydro_sink[ipMH2g] += (realnum)co.H2_SiP_SiH2P_nu;
04619
04620
04621
04622 co.H2_O2_OH_OH = HMRATE(3.16e-10,0,21890)*findspecies("O2")->hevmol*hmi.Hmolec[ipMH2g];
04623 co.hydro_sink[ipMH2g] += (realnum)co.H2_O2_OH_OH;
04624
04625
04626
04627 co.HeP_CH_CP_He_H = HMRATE(1.1e-9,0,0)*findspecies("CH")->hevmol*dense.xIonDense[ipHELIUM][1];
04628 co.hydro_source[ipMH] += (realnum)co.HeP_CH_CP_He_H;
04629
04630
04631
04632 co.HeP_CH2_CHP_He_H = HMRATE(7.5e-10,0,0)*findspecies("CH2")->hevmol*dense.xIonDense[ipHELIUM][1];
04633 co.hydro_source[ipMH] += (realnum)co.HeP_CH2_CHP_He_H;
04634
04635
04636
04637 co.HeP_OH_OP_He_H = HMRATE(1.1e-9,0,0)*findspecies("OH")->hevmol*dense.xIonDense[ipHELIUM][1];
04638 co.hydro_source[ipMH] += (realnum)co.HeP_OH_OP_He_H;
04639
04640
04641
04642 co.HeP_H2O_OHP_He_H = HMRATE(2.86e-10,0,0)*findspecies("H2O")->hevmol*dense.xIonDense[ipHELIUM][1];
04643 co.hydro_source[ipMH] += (realnum)co.HeP_H2O_OHP_He_H;
04644
04645
04646
04647 co.HeP_SiH_SiP_He_H = HMRATE(1.8e-9,0,0)*findspecies("SiH")->hevmol*dense.xIonDense[ipHELIUM][1];
04648 co.hydro_source[ipMH] += (realnum)co.HeP_SiH_SiP_He_H;
04649
04650
04651
04652 co.HeP_H2O_OH_He_HP = HMRATE(2.04e-10,0,0)*findspecies("H2O")->hevmol*dense.xIonDense[ipHELIUM][1];
04653 co.hydro_source[ipMHp] += (realnum)co.HeP_H2O_OH_He_HP;
04654
04655
04656
04657 co.HeP_CH2_CP_He_H2 = HMRATE(7.5e-10,0,0)*findspecies("CH2")->hevmol*dense.xIonDense[ipHELIUM][1];
04658 co.hydro_source[ipMH2g] += (realnum)co.HeP_CH2_CP_He_H2;
04659
04660
04661
04662 co.crnu_CH_C_H = findspecies("CH")->hevmol*secondaries.csupra[ipHYDROGEN ][0] * 2. * 756;
04663 co.hydro_source[ipMH] += (realnum)co.crnu_CH_C_H;
04664
04665
04666
04667 co.crnu_CHP_CP_H = findspecies("CH+")->hevmol*secondaries.csupra[ipHYDROGEN][0] * 2. * 183;
04668 co.hydro_source[ipMH] += (realnum)co.crnu_CHP_CP_H;
04669
04670
04671
04672 co.crnu_H2O_OH_H = findspecies("H2O")->hevmol*secondaries.csupra[ipHYDROGEN][0] * 2. * 979;
04673 co.hydro_source[ipMH] += (realnum)co.crnu_H2O_OH_H;
04674
04675
04676
04677 co.crnu_OH_O_H = findspecies("OH")->hevmol*secondaries.csupra[ipHYDROGEN][0] * 2. *522;
04678 co.hydro_source[ipMH] += (realnum)co.crnu_OH_O_H;
04679
04680
04681
04682 co.crnu_SiH_Si_H = findspecies("SiH")->hevmol*secondaries.csupra[ipHYDROGEN][0] * 2. *500;
04683 co.hydro_source[ipMH] += (realnum)co.crnu_SiH_Si_H;
04684
04685
04686
04687 co.nu_CH_C_H = HMRATE(8.6e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH")->hevmol;
04688 co.hydro_source[ipMH] += (realnum)co.nu_CH_C_H;
04689
04690
04691
04692 co.nu_CHP_CP_H = HMRATE(2.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH+")->hevmol;
04693 co.hydro_source[ipMH] += (realnum)co.nu_CHP_CP_H;
04694
04695
04696
04697 co.nu_CH2_CH_H = HMRATE(7.2e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH2")->hevmol;
04698 co.hydro_source[ipMH] += (realnum)co.nu_CH2_CH_H;
04699
04700
04701
04702 co.nu_CH2P_CHP_H = HMRATE(1.7e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH2+")->hevmol;
04703 co.hydro_source[ipMH] += (realnum)co.nu_CH2P_CHP_H;
04704
04705
04706
04707 co.nu_CH3P_CH2P_H = HMRATE(1.0e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH3+")->hevmol;
04708 co.hydro_source[ipMH] += (realnum)co.nu_CH3P_CH2P_H;
04709
04710
04711
04712 co.nu_CH3P_CHP_H2 = HMRATE(1.0e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH3+")->hevmol;
04713 co.hydro_source[ipMH2g] += (realnum)co.nu_CH3P_CHP_H2;
04714
04715
04716
04717 co.nu_H2O_OH_H = HMRATE(5.9e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("H2O")->hevmol;
04718 co.hydro_source[ipMH] += (realnum)co.nu_H2O_OH_H;
04719
04720
04721
04722 co.nu_OH_O_H = HMRATE(3.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("OH")->hevmol;
04723 co.hydro_source[ipMH] += (realnum)co.nu_OH_O_H;
04724
04725
04726
04727 co.nu_OHP_O_HP = HMRATE(1.0e-12,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("OH+")->hevmol;
04728 co.hydro_source[ipMHp] += (realnum)co.nu_OHP_O_HP;
04729
04730
04731
04732 co.nu_SiH_Si_H = HMRATE(2.8e-9,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("SiH")->hevmol;
04733 co.hydro_source[ipMH] += (realnum)co.nu_SiH_Si_H;
04734
04735
04736
04737 co.e_CHP_C_H = HMRATE(1.5e-7,-0.42,0)*dense.eden*findspecies("CH+")->hevmol;
04738 co.hydro_source[ipMH] += (realnum)co.e_CHP_C_H;
04739
04740
04741
04742 co.e_CH2P_CH_H = HMRATE(1.6e-7,-0.6,0)*dense.eden*findspecies("CH2+")->hevmol;
04743 co.hydro_source[ipMH] += (realnum)co.e_CH2P_CH_H;
04744
04745
04746
04747 co.e_CH2P_C_H_H = HMRATE(4.03e-7,-0.6,0)*dense.eden*findspecies("CH2+")->hevmol;
04748 co.hydro_source[ipMH] += (realnum)co.e_CH2P_C_H_H;
04749
04750
04751
04752 co.e_CH2P_C_H2 = HMRATE(7.68e-8,-0.6,0)*dense.eden*findspecies("CH2+")->hevmol;
04753 co.hydro_source[ipMH2g] += (realnum)co.e_CH2P_C_H2;
04754
04755
04756
04757 co.e_CH3P_C_H2_H = HMRATE(1.05e-7,-0.5,0)*dense.eden*findspecies("CH3+")->hevmol;
04758 co.hydro_source[ipMH] += (realnum)co.e_CH3P_C_H2_H;
04759 co.hydro_source[ipMH2g] += (realnum)co.e_CH3P_C_H2_H;
04760
04761
04762
04763 co.e_CH3P_CH2_H = HMRATE(1.4e-7,-0.5,0)*dense.eden*findspecies("CH3+")->hevmol;
04764 co.hydro_source[ipMH] += (realnum)co.e_CH3P_CH2_H;
04765
04766
04767
04768 co.e_CH3P_CH_H_H = HMRATE(5.6e-8,-0.5,0)*dense.eden*findspecies("CH3+")->hevmol;
04769 co.hydro_source[ipMH] += (realnum)co.e_CH3P_CH_H_H;
04770
04771
04772
04773 co.e_CH3P_CH_H2 = HMRATE(4.9e-8,-0.5,0)*dense.eden*findspecies("CH3+")->hevmol;
04774 co.hydro_source[ipMH2g] += (realnum)co.e_CH3P_CH_H2;
04775
04776
04777
04778 co.e_H2OP_OH_H = HMRATE(7.92e-8,-0.5,0)*dense.eden*findspecies("H2O+")->hevmol;
04779 co.hydro_source[ipMH] += (realnum)co.e_H2OP_OH_H;
04780
04781
04782
04783 co.e_H2OP_O_H_H = HMRATE(2.45e-7,-0.5,0)*dense.eden*findspecies("H2O+")->hevmol;
04784 co.hydro_source[ipMH] += (realnum)co.e_H2OP_O_H_H;
04785
04786
04787
04788 co.e_H2OP_O_H2 = HMRATE(3.6e-8,-0.5,0)*dense.eden*findspecies("H2O+")->hevmol;
04789 co.hydro_source[ipMH2g] += (realnum)co.e_H2OP_O_H2;
04790
04791
04792
04793 co.e_H3OP_H2O_H = HMRATE(1.08e-7,-0.5,0)*dense.eden*findspecies("H3O+")->hevmol;
04794 co.hydro_source[ipMH] += (realnum)co.e_H3OP_H2O_H;
04795
04796
04797
04798 co.e_H3OP_OH_H_H = HMRATE(2.58e-7,-0.5,0)*dense.eden*findspecies("H3O+")->hevmol;
04799 co.hydro_source[ipMH] += (realnum)co.e_H3OP_OH_H_H;
04800
04801
04802
04803 co.e_H3OP_OH_H2 = HMRATE(6.45e-8,-0.5,0)*dense.eden*findspecies("H3O+")->hevmol;
04804 co.hydro_source[ipMH2g] += (realnum)co.e_H3OP_OH_H2;
04805
04806
04807
04808 co.e_H3OP_O_H2_H = HMRATE(5.59e-9,-0.5,0)*dense.eden*findspecies("H3O+")->hevmol;
04809 co.hydro_source[ipMH] += (realnum)co.e_H3OP_O_H2_H;
04810
04811
04812
04813 co.e_HCOP_CO_H = HMRATE(1.1e-7,-1,0)*dense.eden*findspecies("HCO+")->hevmol;
04814 co.hydro_source[ipMH] += (realnum)co.e_HCOP_CO_H;
04815
04816
04817
04818 co.e_OHP_O_H = HMRATE(3.75e-8,-0.5,0)*dense.eden*findspecies("OH+")->hevmol;
04819 co.hydro_source[ipMH] += (realnum)co.e_OHP_O_H;
04820
04821
04822
04823 co.e_SiH2P_SiH_H = HMRATE(1.5e-7,-0.5,0)*dense.eden*findspecies("SiH2+")->hevmol;
04824 co.hydro_source[ipMH] += (realnum)co.e_SiH2P_SiH_H;
04825
04826
04827
04828 co.e_SiH2P_Si_H_H = HMRATE(2.0e-7,-0.5,0)*dense.eden*findspecies("SiH2+")->hevmol;
04829 co.hydro_source[ipMH] += (realnum)co.e_SiH2P_Si_H_H;
04830
04831
04832
04833 co.e_SiH2P_Si_H2 = HMRATE(1.5e-7,-0.5,0)*dense.eden*findspecies("SiH2+")->hevmol;
04834 co.hydro_source[ipMH2g] += (realnum)co.e_SiH2P_Si_H2;
04835
04836
04837
04838 co.e_SiOHP_SiO_H = HMRATE(1.5e-7,-0.5,0)*dense.eden*findspecies("SiOH+")->hevmol;
04839 co.hydro_source[ipMH] += (realnum)co.e_SiOHP_SiO_H;
04840
04841
04842
04843 co.H2_CH_CH3_nu = HMRATE(5.09E-18,-0.71,11.6)*findspecies("H2")->hevmol*findspecies("CH")->hevmol;
04844 co.hydro_sink[ipMH2g] += (realnum)co.H2_CH_CH3_nu;
04845
04846
04847
04848 co.H2_CH3P_CH5P_nu = HMRATE(1.3e-15,-1,0)*findspecies("CH3+")->hevmol*hmi.Hmolec[ipMH2g];
04849 co.hydro_sink[ipMH2g] += (realnum)co.H2_CH3P_CH5P_nu;
04850
04851
04852
04853 co.H2s_CH_CH3_nu = HMRATE(5.09E-18,-0.71,0)*findspecies("CH")->hevmol*hmi.Hmolec[ipMH2s]*hmi.lgLeiden_Keep_ipMH2s;
04854 co.hydro_sink[ipMH2s] += (realnum)co.H2s_CH_CH3_nu;
04855
04856
04857
04858 co.Hminus_CH2_CH3_e = HMRATE(1.0e-9,0,0)*findspecies("CH2")->hevmol*hmi.Hmolec[ipMHm];
04859 co.hydro_sink[ipMHm] += (realnum)co.Hminus_CH2_CH3_e;
04860
04861
04862
04863 co.Hminus_CH3_CH4_e = HMRATE(1.0e-9,0,0)*findspecies("CH3")->hevmol*hmi.Hmolec[ipMHm];
04864 co.hydro_sink[ipMHm] += (realnum)co.Hminus_CH3_CH4_e;
04865
04866
04867
04868 co.nu_CH3_CH2_H = HMRATE(2.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH3")->hevmol;
04869 co.hydro_source[ipMH] += (realnum)co.nu_CH3_CH2_H;
04870
04871
04872
04873 co.nu_CH3_CH_H2 = HMRATE(2.5e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH3")->hevmol;
04874 co.hydro_source[ipMH2g] += (realnum)co.nu_CH3_CH_H2;
04875
04876
04877
04878 co.nu_CH4_CH3_H = HMRATE(2.2e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH4")->hevmol;
04879 co.hydro_source[ipMH] += (realnum)co.nu_CH4_CH3_H;
04880
04881
04882
04883 co.nu_CH4_CH2_H2 = HMRATE(9.8e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH4")->hevmol;
04884 co.hydro_source[ipMH2g] += (realnum)co.nu_CH4_CH2_H2;
04885
04886
04887
04888 co.nu_CH4_CH_H2 = HMRATE(2.2e-10,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_depth/1.66)*findspecies("CH4")->hevmol;
04889 co.hydro_source[ipMH2g] += (realnum)co.nu_CH4_CH_H2;
04890
04891
04892
04893 co.crnu_CH3_CH2_H = findspecies("CH3")->hevmol*secondaries.csupra[ipHYDROGEN][0] * 2.*500;
04894 co.hydro_source[ipMH] += (realnum)co.crnu_CH3_CH2_H;
04895
04896
04897
04898 co.crnu_CH3_CH_H2 = findspecies("CH3")->hevmol*secondaries.csupra[ipHYDROGEN][0] * 2.*500;
04899 co.hydro_source[ipMH2g] += (realnum)co.crnu_CH3_CH_H2;
04900
04901
04902
04903 co.crnu_CH4_CH2_H2 = findspecies("CH4")->hevmol*secondaries.csupra[ipHYDROGEN][0] * 2.*2272;
04904 co.hydro_source[ipMH2g] += (realnum)co.crnu_CH4_CH2_H2;
04905
04906
04907
04908 co.e_CH5P_CH3_H2 = HMRATE(5.5e-7,-0.3,0)*dense.eden*findspecies("CH5+")->hevmol;
04909 co.hydro_source[ipMH2g] += (realnum)co.e_CH5P_CH3_H2;
04910
04911
04912
04913 co.e_CH5P_CH4_H = HMRATE(5.5e-7,-0.3,0)*dense.eden*findspecies("CH5+")->hevmol;
04914 co.hydro_source[ipMH] += (realnum)co.e_CH5P_CH4_H;
04915
04916
04917
04918 co.e_CH4P_CH3_H = HMRATE(1.75e-7,-0.5,0)*dense.eden*findspecies("CH4+")->hevmol;
04919 co.hydro_source[ipMH] += (realnum)co.e_CH4P_CH3_H;
04920
04921
04922
04923 co.e_CH4P_CH2_H_H = HMRATE(1.75e-7,-0.5,0)*dense.eden*findspecies("CH4+")->hevmol;
04924 co.hydro_source[ipMH] += (realnum)co.e_CH4P_CH2_H_H;
04925 # endif
04926
04927 # if 0
04928 if( iteration>1 )
04929 {
04930 r = rlist;
04931 i = 0;
04932 while (r->next != NULL)
04933 {
04934 ++i;
04935 r = r->next;
04936 fprintf(ioQQQ,"DEBUG r\t%li\t%.3e\n", i, r->rk);
04937 }
04938 }
04939 # endif
04940
04941
04942 for( i=0; i < N_H_MOLEC; i++ )
04943 {
04944 bvec[i] = 0.;
04945 for( j=0; j < N_H_MOLEC; j++ )
04946 {
04947 c[j][i] = 0.;
04948 }
04949 }
04950
04951 for(i=0;i<2;i++)
04952 {
04953 mole.source[ipHYDROGEN][i] = mole.sink[ipHYDROGEN][i] = 0.;
04954 }
04955
04956
04957
04958 r = rlist;
04959
04960
04961
04962
04963
04964
04965
04966
04967
04968
04969
04970
04971
04972
04973
04974
04975
04976
04977
04978
04979
04980
04981
04982
04983
04984
04985
04986
04987
04988
04989 while (r->next != NULL)
04990 {
04991 r = r->next;
04992
04993
04994
04995 rk = r->rk;
04996
04997
04998
04999 ASSERT( !isnan( rk ) );
05000
05001
05002
05003
05004 for(i=0;i<r->nrates;i++)
05005 {
05006 rate_deriv[i] = rk;
05007 for(j=0;j<r->nrates;j++)
05008 {
05009
05010
05011 if(i!=j)
05012 {
05013 rate_deriv[i] *= Hmolec_old[r->rate_species[j]];
05014
05015
05016 ASSERT( !isnan( rate_deriv[i] ) );
05017 }
05018 }
05019 }
05020
05021
05022 rate = rate_deriv[0]*Hmolec_old[r->rate_species[0]];
05023
05024
05025 ASSERT( !isnan( rate ) );
05026
05027
05028 for(i=0;i < r->nreactants;i++)
05029 {
05030 int ok = 0;
05031 for(j=0;j < r->nrates && !ok;j++)
05032 {
05033 if(r->rate_species[j] == r->reactants[i])
05034 {
05035 sinkrate[i] = rate_deriv[j];
05036 ok = true;
05037 }
05038 }
05039 if(!ok)
05040 {
05041
05042
05043
05044
05045
05046
05047
05048 fprintf(ioQQQ,"A chemical rate in hmole was independent of the species it used\n");
05049 fprintf(ioQQQ,"This probably shouldn't happen (so you shouldn't see this message).\n");
05050 cdEXIT(EXIT_FAILURE);
05051 }
05052 }
05053
05054
05055
05056
05057
05058
05059 for(i=0;i<r->nreactants;i++)
05060 {
05061 ratei = r->reactants[i];
05062 bvec[ratei] -= rate;
05063
05064
05065
05066
05067 if(ratei == ipMH || ratei == ipMHp)
05068 mole.sink[ipHYDROGEN][ratei] += sinkrate[i];
05069 }
05070
05071
05072
05073 for(i=0;i<r->nproducts;i++)
05074 {
05075 ratei = r->products[i];
05076 bvec[ratei] += rate;
05077
05078
05079 if(ratei == ipMH || ratei == ipMHp)
05080 {
05081 mole.source[ipHYDROGEN][ratei] += rate;
05082
05083
05084 ASSERT( !isnan( mole.source[ipHYDROGEN][ratei] ) );
05085 }
05086 }
05087
05088
05089
05090
05091
05092
05093
05094
05095
05096
05097
05098
05099
05100
05101
05102
05103
05104
05105
05106
05107
05108
05109
05110
05111
05112
05113
05114
05115
05116
05117
05118
05119
05120
05121
05122
05123
05124
05125
05126
05127
05128
05129
05130
05131
05132
05133
05134
05135
05136
05137
05138
05139
05140
05141
05142
05143
05144
05145
05146
05147
05148
05149 for(j=0;j<r->nrates;j++)
05150 {
05151 ratej = r->rate_species[j];
05152 rated = rate_deriv[j];
05153 for(i=0;i<r->nreactants;i++)
05154 {
05155 c[ratej][r->reactants[i]] -= rated;
05156 }
05157 for(i=0;i<r->nproducts;i++)
05158 {
05159 c[ratej][r->products[i]] += rated;
05160 }
05161 }
05162 }
05163
05164
05165
05166
05167
05168
05169
05170 hmi.H2_rate_destroy = (hmi.Hmolec[ipMH2g] * (-c[ipMH2g][ipMH2g]-c[ipMH2g][ipMH2s]) +
05171 hmi.Hmolec[ipMH2s] * (-c[ipMH2s][ipMH2s]-c[ipMH2s][ipMH2g])) / SDIV(hmi.H2_total);
05172
05173 {
05174
05175 enum {DEBUG_LOC=false};
05176 if( DEBUG_LOC )
05177 {
05178 if( DEBUG_LOC && (nzone > 570) )
05179 {
05180 printsol = 1;
05181 fprintf(ioQQQ,"Temperature %g\n",phycon.te);
05182 fprintf(ioQQQ," Net mol ion rate [%g %g] %g\n",mole.source[ipHYDROGEN][1],mole.sink[ipHYDROGEN][1],
05183 mole.source[ipHYDROGEN][1]-mole.sink[ipHYDROGEN][1]*Hmolec_old[ipMHp]);
05184 }
05185 }
05186 }
05187
05188
05189
05190 desh2p = -c[ipMH2p][ipMH2p];
05191
05192
05193
05194
05197 # if 0
05198
05199 # ifndef NDEBUG
05200 {
05201 double total, mtotal;
05202 for(i=0;i<N_H_MOLEC;i++)
05203 {
05204 total = 0.;
05205 for( j=0;j<N_H_MOLEC;j++)
05206 {
05207 total += c[i][j]*hmi.nProton[j];
05208 }
05209 if( fabs(total) > 1e-5*fabs(c[i][i]*hmi.nProton[i]))
05210 {
05211 fprintf(ioQQQ,"PROBLEM Subtotal1 %.2e\n",fabs(total)/fabs(c[i][i]*hmi.nProton[i]));
05212 fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*hmi.nProton[i]);
05213 }
05214 else if( fabs(total) > 1e-6*fabs(c[i][i]*hmi.nProton[i]) && phycon.te< 1e6 )
05215 {
05216 fprintf(ioQQQ,"NOTE Subtotal1 %.2e Te=%.4e\n",
05217 fabs(total)/fabs(c[i][i]*hmi.nProton[i]),phycon.te);
05218 fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*hmi.nProton[i]);
05219 }
05220 }
05221 total = mtotal = 0.;
05222 for(j=0;j<N_H_MOLEC;j++)
05223 {
05224 total += bvec[j]*hmi.nProton[j];
05225 mtotal += fabs(bvec[j]*hmi.nProton[j]);
05226 }
05227 if(fabs(total) > 1e-30 && fabs(total) > 1e-10*rtot)
05228 {
05229 fprintf(ioQQQ,"PROBLEM Subtotal2 %.2e\n",fabs(total)/mtotal);
05230 fprintf(ioQQQ,"RHS Total %g cf %g\n",total,mtotal);
05231 }
05232 else if(fabs(total) > 1e-7*mtotal)
05233 {
05234 fprintf(ioQQQ,"WARNING zone %li Hmole RHS conservation error %.2e of %.2e\n",nzone,total,mtotal);
05235 fprintf(ioQQQ,"(may be due to high rate equilibrium reactions)\n");
05236 }
05237 }
05238 # endif
05239 # endif
05240
05241
05242 #define MOLMIN 1
05243 #define N_H_MAT (N_H_MOLEC-MOLMIN)
05244
05245
05246
05247
05248
05249 if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection
05250 )
05251 {
05252
05253 ipConserve = -1;
05254
05255 for(i=0;i<N_H_MOLEC;i++)
05256 {
05257 c[i][i] -= dynamics.Rate;
05258 bvec[i] -= (Hmolec_old[i]*dynamics.Rate-dynamics.H2_molec[i]);
05259 }
05260
05261
05262 proton_sum_old = 0.;
05263 for(i=0; i<N_H_MOLEC;i++)
05264 {
05265 proton_sum_old += hmi.nProton[i]*dynamics.H2_molec[i]/dynamics.Rate;
05266 }
05267
05268
05269
05270 for(i=0;i<N_H_MOLEC;i++)
05271 {
05272 c[ipMHp][i] = (Hmolec_old[ipMH]*c[ipMH][i]+Hmolec_old[ipMHp]*c[ipMHp][i])/SDIV(sum_H0_Hp);
05273 c[ipMH][i] = 0.;
05274 }
05275 for(i=1;i<N_H_MOLEC;i++)
05276 {
05277 c[i][ipMHp] += c[i][ipMH];
05278 c[i][ipMH] = 0.;
05279 }
05280 bvec[ipMHp] += bvec[ipMH];
05281 bvec[ipMH] = 0.;
05282 Hmolec_old[ipMHp] += Hmolec_old[ipMH];
05283 Hmolec_old[ipMH] = 0.;
05284 }
05285 else
05286 {
05287
05288
05289
05290 for(i=0;i<N_H_MOLEC;i++)
05291 {
05292
05293
05294 if( sum_H0_Hp > SMALLFLOAT )
05295 c[ipMHp][i] = (Hmolec_old[ipMH]*c[ipMH][i]+Hmolec_old[ipMHp]*c[ipMHp][i])/sum_H0_Hp;
05296 c[ipMH][i] = 0.;
05297 }
05298 Hmolec_old[ipMHp] += Hmolec_old[ipMH];
05299 bvec[ipMH] = Hmolec_old[ipMH] = 0.;
05300 ipConserve = ipMHp;
05301
05302
05303 bvec[ipConserve] = 0.;
05304
05305
05306 proton_sum_old = 0.;
05307 for(i=MOLMIN;i<N_H_MOLEC;i++)
05308 {
05309 c[i][ipConserve] = hmi.nProton[i];
05310 proton_sum_old += hmi.nProton[i]*Hmolec_old[i];
05311 }
05312 }
05313
05314 {
05315
05316 enum {DEBUG_LOC=false};
05317 if( DEBUG_LOC )
05318 {
05319
05320 fprintf( ioQQQ, " HMOLE h2 %.2e h2* %.2e\n" , Hmolec_old[ipMH2g] ,Hmolec_old[ipMH2s] );
05321 }
05322 }
05323
05324
05325 if(printsol || (trace.lgTrace && trace.lgTr_H2_Mole ))
05326 {
05327
05328
05329
05330
05331
05332
05333
05334
05335
05336
05337
05338
05339
05340
05341
05342
05343
05344
05345
05346 fprintf(ioQQQ," MOLE old abundances\t%.2f",fnzone);
05347 for( i=0; i<N_H_MOLEC; i++ )
05348 fprintf(ioQQQ,"\t%.2e", Hmolec_old[i] );
05349 fprintf(ioQQQ,"\n" );
05350
05351
05352 fprintf( ioQQQ, " ");
05353 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05354 {
05355 fprintf( ioQQQ, " %s", hmi.chLab[i] );
05356 }
05357 fprintf( ioQQQ, " bvec \n" );
05358
05359 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05360 {
05361 fprintf( ioQQQ, " MOLE%2ld %s", i-MOLMIN ,hmi.chLab[i] );
05362 for( j=MOLMIN; j < N_H_MOLEC; j++ )
05363 {
05364 fprintf( ioQQQ, "%10.2e", c[j][i] );
05365 }
05366 fprintf( ioQQQ, "%10.2e", bvec[i] );
05367 fprintf( ioQQQ, "\n" );
05368 }
05369 }
05370
05371
05372
05373
05374
05375
05376
05377
05378 if( -c[ipMH2g][ipMH2g] > SMALLFLOAT )
05379 {
05380
05381 timesc.time_H2_Dest_here = -1./c[ipMH2g][ipMH2g];
05382 }
05383 else
05384 {
05385 timesc.time_H2_Dest_here = 0.;
05386 }
05387
05388
05389
05390 timesc.time_H2_Form_here = gv.rate_h2_form_grains_used_total *
05391
05392
05393
05394 dense.gas_phase[ipHYDROGEN]/SDIV(dense.xIonDense[ipHYDROGEN][0]) +
05395 hmi.hminus_rad_attach;
05396
05397 if( timesc.time_H2_Form_here > SMALLFLOAT )
05398 {
05399
05400 timesc.time_H2_Form_here = 1./timesc.time_H2_Form_here;
05401 }
05402 else
05403 {
05404 timesc.time_H2_Form_here = 0.;
05405 }
05406
05407 # ifdef MAT
05408 # undef MAT
05409 # endif
05410 # define MAT(a,I_,J_) (*((a)+(I_)*(N_H_MAT)+(J_)))
05411
05412
05413 for( j=0; j < N_H_MAT; j++ )
05414 {
05415 for( i=0; i < N_H_MAT; i++ )
05416 {
05417 MAT(amat,i,j) = c[i+MOLMIN][j+MOLMIN];
05418 }
05419 }
05420
05421 if(printsol)
05422 {
05423 double total=0;
05424 fprintf(ioQQQ,"Zone %.2f input:\n",fnzone);
05425 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05426 {
05427 fprintf(ioQQQ,"%15.7g\t",Hmolec_old[i]);
05428 total += hmi.nProton[i]*Hmolec_old[i];
05429 }
05430 fprintf(ioQQQ,"sum = %15.7g\n",total);
05431 }
05432
05433 int32 merror1 = 0;
05434 int32 merror2 = 0;
05435
05436
05437 getrf_wrapper(N_H_MAT,N_H_MAT,(double*)amat,N_H_MAT,ipiv,&merror1);
05438 getrs_wrapper('N',N_H_MAT,1,(double*)amat,N_H_MAT,ipiv,bvec+MOLMIN,N_H_MAT,&merror2);
05439
05440 if( merror1 != 0 || merror2 != 0 )
05441 {
05442 fprintf( ioQQQ, "PROBLEM hmole_step: dgetrs finds singular or ill-conditioned matrix\n" );
05443 cdEXIT(EXIT_FAILURE);
05444 }
05445
05446 if(printsol)
05447 {
05448 double total=0;
05449 fprintf(ioQQQ,"solution:\n");
05450 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05451 {
05452 fprintf(ioQQQ,"%15.7g\t",bvec[i]);
05453 total += hmi.nProton[i]*bvec[i];
05454 }
05455 fprintf(ioQQQ,"sum = %15.7g\n",total);
05456 }
05457
05458 *error = 0;
05459
05460
05461
05462
05463
05464
05465 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05466 {
05467
05468 etmp = bvec[i]/(ABSLIM+Hmolec_old[i]);
05469
05470 if(printsol)
05471 fprintf(ioQQQ,"%15.7g\t",etmp);
05472
05473
05474 *error += etmp*etmp;
05475
05476
05477 bvec[i] = Hmolec_old[i]-bvec[i];
05478 }
05479
05480 *error = sqrt(*error)/N_H_MAT;
05481
05482 if(printsol)
05483 {
05484 double total=0;
05485 fprintf(ioQQQ,"err = %15.7g\n",*error);
05486
05487 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05488 {
05489 fprintf(ioQQQ,"%15.7g\t",bvec[i]);
05490 total += hmi.nProton[i]*bvec[i];
05491 }
05492 fprintf(ioQQQ,"sum = %15.7g\n",total);
05493 }
05494
05495 proton_sum_new = 0.;
05496
05497 lgNegPop = false;
05498 fracneg = 0.;
05499 fracnegfac = 0.;
05500 iworst = -1;
05501 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05502 {
05503 if( bvec[i] < 0. )
05504 {
05505 lgNegPop = true;
05506 }
05507
05508 fracnegtmp = -bvec[i]/SDIV(Hmolec_old[i]);
05509
05510
05511 if(fracnegtmp > fracneg)
05512 {
05513 fracneg = fracnegtmp;
05514 fracnegfac = (0.5*Hmolec_old[i]-bvec[i])/(Hmolec_old[i]-bvec[i]);
05515 iworst = i;
05516 }
05517
05518 proton_sum_new += hmi.nProton[i] * bvec[i];
05519 }
05520
05521
05522
05523 conserve = (proton_sum_old - proton_sum_new)/SDIV(proton_sum_old);
05524
05525
05526
05527
05528
05529
05530 if( fabs(conserve) > 10.*FLT_EPSILON )
05531 fprintf(ioQQQ,"PROBLEM hmoleee zn %li proton_sum_old %.8e, proton_sum_new %.8e n(H) %.8e (old-new)/old %.3e nH-old %.3e nH-new %.3e\n",
05532 nzone ,
05533 proton_sum_old ,
05534 proton_sum_new ,
05535 dense.gas_phase[ipHYDROGEN] ,
05536 conserve ,
05537 dense.gas_phase[ipHYDROGEN]-proton_sum_old,
05538 dense.gas_phase[ipHYDROGEN]-proton_sum_new);
05539
05540 # if 0
05541
05542
05543 # ifndef NDEBUG
05544
05545 {
05546 fprintf( ioQQQ, "Zone %li\n",nzone);
05547 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05548 {
05549 fprintf(ioQQQ," %s %.2e", hmi.chLab[i], Hmolec_old[i]);
05550 }
05551 fprintf( ioQQQ, " =>\n" );
05552 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05553 {
05554 fprintf(ioQQQ," %s %.2e", hmi.chLab[i], bvec[i]);
05555 }
05556 fprintf( ioQQQ, "\n" );
05557 }
05558 # endif
05559 # endif
05560
05561 if(lgNegPop)
05562 {
05563 # ifndef NDEBUG
05564
05565
05566 if(*nFixup )
05567 {
05568 fprintf( ioQQQ, " PROBLEM hmole_step finds negative H molecule, in zone %.2f.\n",fnzone );
05569 fprintf( ioQQQ, " Worst is species %d -ve by fraction %g.\n",iworst,fracneg );
05570 fprintf( ioQQQ, " The populations follow:\n");
05571 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05572 {
05573 fprintf(ioQQQ," %s %.2e", hmi.chLab[i], bvec[i]);
05574 }
05575 fprintf( ioQQQ, "\n" );
05576 }
05577 # endif
05578
05579
05580 {
05581 double total=0., ntotal=0., ratio;
05582 enum {FIXUPTYPE = 1};
05583
05584 ++*nFixup;
05585
05586 if(FIXUPTYPE == 1) {
05587 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05588 {
05589 total += hmi.nProton[i]*bvec[i];
05590 if(bvec[i] < 0)
05591 {
05592 ntotal += hmi.nProton[i]*bvec[i];
05593 bvec[i] = 0.;
05594 }
05595 }
05596 ratio = total/(total-ntotal);
05597 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05598 {
05599 bvec[i] *= ratio;
05600 }
05601 }
05602 else if(FIXUPTYPE == 2)
05603 {
05604 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05605 {
05606 bvec[i] = fracnegfac*Hmolec_old[i]+(1-fracnegfac)*bvec[i];
05607 }
05608 }
05609
05610 # ifndef NDEBUG
05611
05612
05613
05614 if( *nFixup>1 )
05615 {
05616 fprintf(ioQQQ," FIXUP taken %i time%s.\n\n", *nFixup, (*nFixup == 1)?"":"s");
05617 fprintf( ioQQQ, " Initially:\n");
05618 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05619 {
05620 fprintf(ioQQQ," %s %.2e", hmi.chLab[i], Hmolec_old[i]);
05621 }
05622 fprintf( ioQQQ, "\n" );
05623 fprintf( ioQQQ, " Changed to:\n");
05624 for( i=MOLMIN; i < N_H_MOLEC; i++ )
05625 {
05626 fprintf(ioQQQ," %s %.2e", hmi.chLab[i], bvec[i]);
05627 }
05628 fprintf( ioQQQ, "\n" );
05629 }
05630 # endif
05631 }
05632 }
05633
05634
05635
05636
05637 h1fnd = bvec[ipMHp];
05638
05639 h1rat = h1fnd/SDIV(sum_H0_Hp);
05640
05641
05642 bvec[ipMH] = dense.xIonDense[ipHYDROGEN][0] * h1rat;
05643 bvec[ipMHp] = dense.xIonDense[ipHYDROGEN][1] * h1rat;
05644
05645
05646 if(fabs(bvec[ipMH]+bvec[ipMHp]-h1fnd) >= 1e-12 * h1fnd)
05647 {
05648 static bool lgPrint=true;
05649 fprintf(ioQQQ,"PROBLEM h1fnd residual error, bvec[ipMH]=%g [ipMHp}=%g"
05650 " h1fnd=%g h1rat=%g bvec[ipMH]+bvec[ipMHp]-h1fnd=%g\n",
05651 bvec[ipMH],bvec[ipMHp],h1fnd,h1rat,bvec[ipMH]+bvec[ipMHp]-h1fnd);
05652
05653
05654 if( lgPrint )
05655 {
05656 fprintf(ioQQQ," This problem is usually caused by little or no sources of ionizaiton.\n");
05657 fprintf(ioQQQ," Is the simulation physically motivated?\n");
05658 if( hextra.cryden==0. && lgPrint )
05659 {
05660 fprintf(ioQQQ,"PROBLEM h1fnd - no cosmic rays are present - is this physical?\n");
05661 fprintf(ioQQQ,"PROBLEM h1fnd - Consider including the COSMIC RAY BACKGROUND command.\n");
05662 }
05663 lgPrint = false;
05664 }
05665 }
05666
05667
05668 for(mol=0;mol<N_H_MOLEC;mol++)
05669 {
05670 hmi.Hmolec[mol] = (realnum) bvec[mol];
05671 }
05672
05673 dense.xIonDense[ipHYDROGEN][0] = (realnum) bvec[ipMH];
05674 dense.xIonDense[ipHYDROGEN][1] = (realnum) bvec[ipMHp];
05675
05676
05677 hmi.H2_total = hmi.Hmolec[ipMH2s] + hmi.Hmolec[ipMH2g];
05678
05679 h2.ortho_density = 0.75*hmi.H2_total;
05680 h2.para_density = 0.25*hmi.H2_total;
05681
05682
05683
05684
05685
05686 dense.xIonDense[LIMELM+2][0] = hmi.H2_total;
05687
05688
05689 {
05690
05691 enum {DEBUG_LOC=false};
05692 if( DEBUG_LOC && (nzone>50) )
05693 {
05694 double createsum ,create_from_Hn2 , create_3body_Ho, create_h2p,
05695 create_h3p, create_h3pe, create_grains, create_hminus;
05696 double destroysum, destroy_hm ,destroy_soloman ,destroy_2h ,destroy_hp,
05697 destroy_h,destroy_hp2,destroy_h3p;
05698
05699
05700 create_from_Hn2 = hmi.radasc*dense.xIonDense[ipHYDROGEN][0];
05701
05702 create_3body_Ho = hmi.bh2dis*dense.xIonDense[ipHYDROGEN][0];
05703
05704 create_h2p = hmi.bh2h2p*dense.xIonDense[ipHYDROGEN][0]*Hmolec_old[ipMH2p];
05705
05706 create_h3p = hmi.h3ph2p*dense.xIonDense[ipHYDROGEN][0]*hmi.Hmolec[ipMH3p];
05707
05708 create_h3pe = hmi.eh3_h2h*dense.eden * hmi.Hmolec[ipMH3p];
05709
05710 create_grains = gv.rate_h2_form_grains_used_total;
05711
05712 create_hminus = Hmolec_old[ipMH]*hmi.assoc_detach*hmi.Hmolec[ipMHm];
05713
05714 createsum = create_from_Hn2 + create_3body_Ho + create_h2p +
05715 create_h3p + create_h3pe + create_grains + create_hminus;
05716
05717 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",
05718 fnzone,
05719 createsum,
05720 create_hminus / createsum,
05721 create_from_Hn2 / createsum,
05722 create_3body_Ho / createsum,
05723 create_h2p / createsum,
05724 create_h3p / createsum,
05725 create_h3pe / createsum,
05726 create_grains / createsum );
05727
05728
05729
05730
05731 destroy_hm = hmi.assoc_detach_backwards_grnd+hmi.assoc_detach_backwards_exct;
05732
05733 destroy_soloman = hmi.H2_Solomon_dissoc_rate_used_H2g;
05734 destroy_2h = hmi.eh2hh;
05735 destroy_hp = hmi.h2hph3p*dense.xIonDense[ipHYDROGEN][1];
05736 destroy_h = hmi.rh2dis*dense.xIonDense[ipHYDROGEN][0];
05737 destroy_hp2 = hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1];
05738 destroy_h3p = hmi.h3petc * hmi.Hmolec[ipMH3p];
05739 destroysum = destroy_hm + destroy_soloman + destroy_2h +
05740 destroy_hp+ destroy_h+ destroy_hp2+ destroy_h3p;
05741
05742 fprintf(ioQQQ,"H2 destroy\t%.3f \t%.2e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
05743 fnzone,
05744 destroysum,
05745 destroy_hm / destroysum ,
05746 destroy_soloman / destroysum ,
05747 destroy_2h / destroysum ,
05748 destroy_hp / destroysum ,
05749 destroy_h / destroysum ,
05750 destroy_hp2 / destroysum ,
05751 destroy_h3p / destroysum );
05752
05753 }
05754 }
05755
05756 {
05757
05758 enum {DEBUG_LOC=false};
05759 if( DEBUG_LOC && (nzone>140) )
05760 {
05761 double create_from_Ho,create_3body_Ho,create_batach,destroy_photo,
05762 destroy_coll_heavies,destroy_coll_electrons,destroy_Hattach,destroy_fhneut,
05763 destsum , createsum;
05764
05765 create_from_Ho = (hmi.hminus_rad_attach + hmi.HMinus_induc_rec_rate);
05766 create_3body_Ho = c3bod;
05767
05768 create_batach = hmi.assoc_detach_backwards_grnd + hmi.assoc_detach_backwards_exct;
05769 destroy_photo = hmi.HMinus_photo_rate;
05770 destroy_coll_heavies = hmi.hmin_ct_firstions*sum_first_ions;
05771 destroy_coll_electrons = cionhm;
05772 destroy_Hattach = Hmolec_old[ipMH]*hmi.assoc_detach;
05773 destroy_fhneut = fhneut;
05774
05775 destsum = destroy_photo + destroy_coll_heavies + destroy_coll_electrons +
05776 destroy_Hattach + destroy_fhneut;
05777 fprintf(ioQQQ,"H- destroy zone\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n",
05778 fnzone,
05779 phycon.te,
05780 destsum,
05781 destroy_photo/destsum ,
05782 destroy_coll_heavies/destsum,
05783 destroy_coll_electrons/destsum,
05784 destroy_Hattach/destsum,
05785 destroy_fhneut/destsum );
05786
05787 createsum = create_from_Ho+create_3body_Ho+create_batach;
05788 fprintf(ioQQQ,"H- create\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
05789 fnzone,
05790 phycon.te,
05791 createsum,
05792 dense.eden,
05793 create_from_Ho/createsum,
05794 create_3body_Ho/createsum,
05795 create_batach/createsum);
05796 }
05797 }
05798
05799
05800 hmi.HalphaHmin = (realnum)(fhneut*hmi.Hmolec[ipMHm]);
05801
05802
05803 if( hmi.lgNoH2Mole )
05804 {
05805 hmi.HeatH2Dish_TH85 = 0.;
05806 hmi.HeatH2Dexc_TH85 = 0.;
05807 hmi.deriv_HeatH2Dexc_TH85 = 0.;
05808 }
05809 else
05810 {
05811
05812
05813
05814
05815
05816
05817
05818
05819
05820
05821
05822
05823
05824
05825
05826
05827
05828
05829
05830
05831
05832
05833
05834 hmi.HeatH2Dish_TH85 = 0.25 * EN1EV *
05835 (hmi.H2_Solomon_dissoc_rate_used_H2g * hmi.Hmolec[ipMH2g] +
05836 hmi.H2_Solomon_dissoc_rate_used_H2s * hmi.Hmolec[ipMH2s]);
05837
05838
05839 hmi.HeatH2Dexc_TH85 = (hmi.Hmolec[ipMH2s]*H2star_deexcit - hmi.Hmolec[ipMH2g]*H2star_excit) * 4.17e-12;
05840
05841 hmi.deriv_HeatH2Dexc_TH85 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_TH85)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
05842
05843 if( hmi.chH2_small_model_type == 'H' )
05844 {
05845
05846 hmi.HeatH2Dish_BHT90 = 0.25 * EN1EV *
05847 (hmi.H2_Solomon_dissoc_rate_used_H2g * hmi.Hmolec[ipMH2g] +
05848 hmi.H2_Solomon_dissoc_rate_used_H2s * hmi.Hmolec[ipMH2s]);
05849
05850
05851 hmi.HeatH2Dexc_BHT90 = (hmi.Hmolec[ipMH2s]*H2star_deexcit - hmi.Hmolec[ipMH2g]*H2star_excit) * 4.17e-12;
05852
05853 hmi.deriv_HeatH2Dexc_BHT90 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BHT90)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
05854 }
05855 else if( hmi.chH2_small_model_type == 'B')
05856 {
05857
05858 hmi.HeatH2Dish_BD96 = 0.25 * EN1EV *
05859 (hmi.H2_Solomon_dissoc_rate_used_H2g * hmi.Hmolec[ipMH2g] +
05860 hmi.H2_Solomon_dissoc_rate_used_H2s * hmi.Hmolec[ipMH2s]);
05861
05862 hmi.HeatH2Dexc_BD96 = (hmi.Hmolec[ipMH2s]*H2star_deexcit - hmi.Hmolec[ipMH2g]*H2star_excit) * 4.17e-12;
05863
05864 hmi.deriv_HeatH2Dexc_BD96 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BD96)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
05865 }
05866 else if(hmi.chH2_small_model_type == 'E')
05867 {
05868
05869
05870
05871
05872
05873 double log_density,
05874 f1, f2,f3, f4, f5;
05875 static double log_G0_face = -1;
05876 static double k_f4;
05877
05878
05879
05880
05881 if( !nzone )
05882 {
05883 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
05884 {
05885 log_G0_face = 0.;
05886 }
05887 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
05888 {
05889 log_G0_face = 7.;
05890 }
05891 else
05892 {
05893 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
05894 }
05895
05896 log_G0_face /= radius.r1r0sq;
05897 }
05898
05899 if(dense.gas_phase[ipHYDROGEN] <= 1.)
05900 {
05901 log_density = 0.;
05902 }
05903 else if(dense.gas_phase[ipHYDROGEN] >= 1e7)
05904 {
05905 log_density = 7.;
05906 }
05907 else
05908 {
05909 log_density = log10(dense.gas_phase[ipHYDROGEN]);
05910 }
05911
05912 f1 = 0.15 * log_density + 0.75;
05913 f2 = -0.5 * log_density + 10.;
05914
05915 hmi.HeatH2Dish_ELWERT = 0.25 * EN1EV * f1 *
05916 (hmi.H2_Solomon_dissoc_rate_used_H2g * hmi.Hmolec[ipMH2g] +
05917 hmi.H2_Solomon_dissoc_rate_used_H2s * hmi.Hmolec[ipMH2s] ) +
05918 f2 * secondaries.x12tot * EN1EV * hmi.H2_total;
05919
05920
05921
05922
05923
05924
05925
05926
05927
05928
05929 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
05930 {
05931 log_G0_face = 0.;
05932 }
05933 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
05934 {
05935 log_G0_face = 7.;
05936 }
05937 else
05938 {
05939 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
05940 }
05941
05942 log_G0_face /= radius.r1r0sq;
05943
05944
05945 k_f4 = -0.25 * log_G0_face + 1.25;
05946
05947
05948 if(dense.gas_phase[ipHYDROGEN] <= 1.)
05949 {
05950 log_density = 0.;
05951 f4 = 0.;
05952 }
05953 else if(dense.gas_phase[ipHYDROGEN] >= 1e7)
05954 {
05955 log_density = 7.;
05956 f4 = pow(k_f4,2) * pow( 10. , 2.2211 * log_density - 29.8506);
05957 }
05958 else
05959 {
05960 log_density = log10(dense.gas_phase[ipHYDROGEN]);
05961 f4 = pow(k_f4,2) * pow( 10., 2.2211 * log_density - 29.8506);
05962 }
05963
05964 f3 = MAX2(0.1, -4.75 * log_density + 24.25);
05965 f5 = MAX2(1.,0.95 * log_density - 1.45) * 0.2 * log_G0_face;
05966
05967 hmi.HeatH2Dexc_ELWERT = (hmi.Hmolec[ipMH2s]*H2star_deexcit - hmi.Hmolec[ipMH2g]*H2star_excit) * 4.17e-12 * f3 +
05968 f4 * (hmi.Hmolec[ipMH]/dense.gas_phase[ipHYDROGEN]) +
05969 f5 * secondaries.x12tot * EN1EV * hmi.H2_total;
05970
05971 if(log_G0_face == 0.&& dense.gas_phase[ipHYDROGEN] > 1.)
05972 hmi.HeatH2Dexc_ELWERT *= 0.001 / dense.gas_phase[ipHYDROGEN];
05973
05974
05975
05976 hmi.HeatH2Dexc_ELWERT /= POW2(radius.r1r0sq);
05977
05978
05979 hmi.deriv_HeatH2Dexc_ELWERT = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_ELWERT)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
05980
05981
05982 }
05983
05984 else
05985 TotalInsanity();
05986
05987 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
05988 {
05989 deexc_htwo = hmi.Average_collH2_deexcit;
05990 deexc_hneut = hmi.Average_collH_deexcit;
05991 }
05992 else
05993 {
05994 deexc_htwo = (1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6.;
05995 deexc_hneut = (1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.;
05996 }
05997
05998 H2star_deexcit = hmi.H2_total*deexc_htwo + hmi.Hmolec[ipMH] * deexc_hneut;
05999
06000 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
06001 {
06002 H2star_excit = hmi.Average_collH2_excit *hmi.H2_total + hmi.Average_collH_excit*hmi.Hmolec[ipMH];
06003 }
06004 else
06005 {
06006 H2star_excit = Boltz_fac_H2_H2star * H2star_deexcit;
06007 }
06008
06009
06010
06011
06012
06013
06014
06015
06016
06017
06018
06019
06020
06021
06022
06023
06024
06025
06026 }
06027
06028 {
06029
06030 enum {DEBUG_LOC=false};
06031 if( DEBUG_LOC )
06032 {
06033
06034 fprintf( ioQQQ, " HMOLE raw; hi\t%.2e" , dense.xIonDense[ipHYDROGEN][0]);
06035 for( i=0; i < N_H_MOLEC; i++ )
06036 {
06037 fprintf( ioQQQ, "\t%s\t%.2e", hmi.chLab[i], bvec[i] );
06038 }
06039 fprintf( ioQQQ, " \n" );
06040 }
06041 }
06042
06043 if( trace.lgTrace && trace.lgTr_H2_Mole )
06044 {
06045
06046 fprintf( ioQQQ, "\n raw; " );
06047 for( i=0; i < N_H_MOLEC; i++ )
06048 {
06049 fprintf( ioQQQ, " %s:%.2e", hmi.chLab[i], bvec[i] );
06050 }
06051 fprintf( ioQQQ, " \n" );
06052 }
06053
06054
06055
06056 hmi.H2_rate_create = gv.rate_h2_form_grains_used_total * dense.xIonDense[ipHYDROGEN][0] +
06057 hmi.assoc_detach * hmi.Hmolec[ipMH] * hmi.Hmolec[ipMHm] +
06058 hmi.bh2dis * dense.xIonDense[ipHYDROGEN][0] +
06059 hmi.bh2h2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2p] +
06060 hmi.radasc * dense.xIonDense[ipHYDROGEN][0] +
06061 hmi.h3ph2p * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH3p] +
06062 hmi.h2phmh2h * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMHm] +
06063 hmi.bh2h22hh2 * 2 * dense.xIonDense[ipHYDROGEN][0] * hmi.Hmolec[ipMH2g] +
06064 hmi.h3phmh2hh * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] +
06065 hmi.h3phm2h2 * hmi.Hmolec[ipMH3p] * hmi.Hmolec[ipMHm] +
06066 hmi.h32h2 * hmi.Hmolec[ipMH2p] * hmi.Hmolec[ipMH3p] +
06067 hmi.eh3_h2h * hmi.Hmolec[ipMH3p] +
06068 hmi.h3ph2hp * hmi.Hmolec[ipMH3p]+
06069 co.H_CH_C_H2 * dense.xIonDense[ipHYDROGEN][0] +
06070 co.H_CHP_CP_H2 * dense.xIonDense[ipHYDROGEN][0] +
06071 co.H_CH2_CH_H2 * dense.xIonDense[ipHYDROGEN][0] +
06072 co.H_CH3P_CH2P_H2 * dense.xIonDense[ipHYDROGEN][0] +
06073 co.H_OH_O_H2 * dense.xIonDense[ipHYDROGEN][0] +
06074 co.Hminus_HCOP_CO_H2 * hmi.Hmolec[ipMHm] +
06075 co.Hminus_H3OP_H2O_H2 * hmi.Hmolec[ipMHm] +
06076 co.Hminus_H3OP_OH_H2_H * hmi.Hmolec[ipMHm] +
06077 co.HP_CH2_CHP_H2 * hmi.Hmolec[ipMHp] +
06078 co.HP_SiH_SiP_H2 * hmi.Hmolec[ipMHp] +
06079 co.H2P_CH_CHP_H2 * hmi.Hmolec[ipMH2p] +
06080 co.H2P_CH2_CH2P_H2 * hmi.Hmolec[ipMH2p] +
06081 co.H2P_CO_COP_H2 * hmi.Hmolec[ipMH2p] +
06082 co.H2P_H2O_H2OP_H2 * hmi.Hmolec[ipMH2p] +
06083 co.H2P_O2_O2P_H2 * hmi.Hmolec[ipMH2p] +
06084 co.H2P_OH_OHP_H2 * hmi.Hmolec[ipMH2p] +
06085 co.H3P_C_CHP_H2 * hmi.Hmolec[ipMH3p] +
06086 co.H3P_CH_CH2P_H2 * hmi.Hmolec[ipMH3p] +
06087 co.H3P_CH2_CH3P_H2 * hmi.Hmolec[ipMH3p] +
06088 co.H3P_OH_H2OP_H2 * hmi.Hmolec[ipMH3p] +
06089 co.H3P_H2O_H3OP_H2 * hmi.Hmolec[ipMH3p] +
06090 co.H3P_CO_HCOP_H2 * hmi.Hmolec[ipMH3p] +
06091 co.H3P_O_OHP_H2 * hmi.Hmolec[ipMH3p] +
06092 co.H3P_SiH_SiH2P_H2 * hmi.Hmolec[ipMH3p] +
06093 co.H3P_SiO_SiOHP_H2 * hmi.Hmolec[ipMH3p] +
06094 co.H_CH3_CH2_H2 * dense.xIonDense[ipHYDROGEN][0] +
06095 co.H_CH4P_CH3P_H2 * dense.xIonDense[ipHYDROGEN][0] +
06096 co.H_CH5P_CH4P_H2 * dense.xIonDense[ipHYDROGEN][0] +
06097 co.H2P_CH4_CH3P_H2 * hmi.Hmolec[ipMH2p] +
06098 co.H2P_CH4_CH4P_H2 * hmi.Hmolec[ipMH2p] +
06099 co.H3P_CH3_CH4P_H2 * hmi.Hmolec[ipMH3p] +
06100 co.H3P_CH4_CH5P_H2 * hmi.Hmolec[ipMH3p] +
06101 co.HP_CH4_CH3P_H2 * hmi.Hmolec[ipMHp] +
06102 co.HP_HNO_NOP_H2 * hmi.Hmolec[ipMHp] +
06103 co.HP_HS_SP_H2 * hmi.Hmolec[ipMHp] +
06104 co.H_HSP_SP_H2 * dense.xIonDense[ipHYDROGEN][0] +
06105 co.H3P_NH_NH2P_H2 * hmi.Hmolec[ipMH3p] +
06106 co.H3P_NH2_NH3P_H2 * hmi.Hmolec[ipMH3p] +
06107 co.H3P_NH3_NH4P_H2 * hmi.Hmolec[ipMH3p] +
06108 co.H3P_CN_HCNP_H2 * hmi.Hmolec[ipMH3p] +
06109 co.H3P_NO_HNOP_H2 * hmi.Hmolec[ipMH3p] +
06110 co.H3P_S_HSP_H2 * hmi.Hmolec[ipMH3p] +
06111 co.H3P_CS_HCSP_H2 * hmi.Hmolec[ipMH3p] +
06112 co.H3P_NO2_NOP_OH_H2 * hmi.Hmolec[ipMH3p] +
06113 co.H2P_NH_NHP_H2 * hmi.Hmolec[ipMH2p] +
06114 co.H2P_NH2_NH2P_H2 * hmi.Hmolec[ipMH2p] +
06115 co.H2P_NH3_NH3P_H2 * hmi.Hmolec[ipMH2p] +
06116 co.H2P_CN_CNP_H2 * hmi.Hmolec[ipMH2p] +
06117 co.H2P_HCN_HCNP_H2 * hmi.Hmolec[ipMH2p] +
06118 co.H2P_NO_NOP_H2 * hmi.Hmolec[ipMH2p] +
06119 co.H3P_Cl_HClP_H2 * hmi.Hmolec[ipMH3p]+
06120 co.H3P_HCl_H2ClP_H2 * hmi.Hmolec[ipMH3p]+
06121 co.H2P_C2_C2P_H2 * hmi.Hmolec[ipMH2p]+
06122 co.Hminus_NH4P_NH3_H2 * hmi.Hmolec[ipMHm]+
06123 co.H3P_HCN_HCNHP_H2 * hmi.Hmolec[ipMH3p];
06124
06125
06126
06127
06128 if( (trace.lgTrace && trace.lgTr_H2_Mole) )
06129 {
06130
06131 if( hmi.H2_rate_create > SMALLFLOAT )
06132 {
06133 fprintf( ioQQQ,
06134 " Create H2, rate=%10.2e grain;%5.3f hmin;%5.3f bhedis;%5.3f h2+;%5.3f hmi.radasc:%5.3f hmi.h3ph2p:%5.3f hmi.h3petc:%5.3f\n",
06135 hmi.H2_rate_create,
06136 gv.rate_h2_form_grains_used_total/hmi.H2_rate_create,
06137 hmi.assoc_detach * hmi.Hmolec[ipMH] * hmi.Hmolec[ipMHm] /hmi.H2_rate_create,
06138 hmi.bh2dis * dense.xIonDense[ipHYDROGEN][0]/hmi.H2_rate_create,
06139 hmi.bh2h2p*dense.xIonDense[ipHYDROGEN][0]*hmi.Hmolec[ipMH2p]/hmi.H2_rate_create,
06140 hmi.radasc*dense.xIonDense[ipHYDROGEN][0]/hmi.H2_rate_create,
06141 hmi.h3ph2p*hmi.Hmolec[ipMH3p]/hmi.H2_rate_create,
06142 hmi.h3petc*hmi.Hmolec[ipMH3p]/hmi.H2_rate_create );
06143 }
06144 else
06145 {
06146 fprintf( ioQQQ, " Create H2, rate=0\n" );
06147 }
06148 }
06149
06150
06151 if( trace.lgTrace && trace.lgTr_H2_Mole )
06152 {
06153 rate = hmi.rh2h2p*hmi.Hmolec[ipMH2g]*dense.xIonDense[ipHYDROGEN][1] + b2pcin*dense.xIonDense[ipHYDROGEN][1]*dense.eden*dense.xIonDense[ipHYDROGEN][0] +
06154 hmi.h3ph2p*hmi.Hmolec[ipMH3p] + hmi.h3petc*hmi.Hmolec[ipMH3p];
06155 if( rate > 1e-25 )
06156 {
06157 fprintf( ioQQQ, " Create H2+, rate=%10.2e hmi.rh2h2p;%5.3f b2pcin;%5.3f hmi.h3ph2p;%5.3f hmi.h3petc+;%5.3f\n",
06158 rate, hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1]*hmi.Hmolec[ipMH2g]/rate,
06159 b2pcin*dense.xIonDense[ipHYDROGEN][1]*dense.xIonDense[ipHYDROGEN][1]*
06160 dense.eden/rate, hmi.h3ph2p*hmi.Hmolec[ipMH3p]/rate, hmi.h3petc*hmi.Hmolec[ipMH3p]/
06161 rate );
06162 }
06163 else
06164 {
06165 fprintf( ioQQQ, " Create H2+, rate=0\n" );
06166 }
06167 }
06168
06169 if( hmi.Hmolec[ipMHm] > 0. && hmi.rel_pop_LTE_Hmin > 0. )
06170 {
06171 hmi.hmidep = (double)hmi.Hmolec[ipMHm]/ SDIV(
06172 ((double)dense.xIonDense[ipHYDROGEN][0]*dense.eden*hmi.rel_pop_LTE_Hmin));
06173 }
06174 else
06175 {
06176 hmi.hmidep = 1.;
06177 }
06178
06179
06180 hmi.hmihet = hmi.HMinus_photo_heat*hmi.Hmolec[ipMHm] - hmi.HMinus_induc_rec_cooling;
06181 hmi.h2plus_heat = h2phet*hmi.Hmolec[ipMH2p];
06182
06183
06184
06185 plte = (double)dense.xIonDense[ipHYDROGEN][0] * hmi.rel_pop_LTE_H2g * (double)dense.xIonDense[ipHYDROGEN][0];
06186 if( plte > 0. )
06187 {
06188 hmi.h2dep = hmi.Hmolec[ipMH2g]/plte;
06189 }
06190 else
06191 {
06192 hmi.h2dep = 1.;
06193 }
06194
06195
06196
06197 plte = (double)dense.xIonDense[ipHYDROGEN][0]*hmi.rel_pop_LTE_H2p*(double)dense.xIonDense[ipHYDROGEN][1];
06198 if( plte > 0. )
06199 {
06200 hmi.h2pdep = hmi.Hmolec[ipMH2p]/plte;
06201 }
06202 else
06203 {
06204 hmi.h2pdep = 1.;
06205 }
06206
06207
06208 if( hmi.rel_pop_LTE_H3p > 0. )
06209 {
06210 hmi.h3pdep = hmi.Hmolec[ipMH3p]/hmi.rel_pop_LTE_H3p;
06211 }
06212 else
06213 {
06214 hmi.h3pdep = 1.;
06215 }
06216
06217
06218 if( trace.lgTrace && trace.lgTr_H2_Mole )
06219 {
06220 fprintf( ioQQQ, " HMOLE, Dep Coef, H-:%10.2e H2:%10.2e H2+:%10.2e\n",
06221 hmi.hmidep, hmi.h2dep, hmi.h2pdep );
06222 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",
06223 hmi.hminus_rad_attach, hmi.HMinus_induc_rec_rate, bhneut, c3bod, hmi.assoc_detach_backwards_grnd, hmi.Hmolec[ipMHm], hmi.hmidep );
06224
06225 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",
06226 hmi.HMinus_photo_rate, hmi.hmin_ct_firstions*sum_first_ions, cionhm,
06227 Hmolec_old[ipMH]*hmi.assoc_detach, iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s],
06228 fhneut );
06229 fprintf( ioQQQ, " H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n",
06230 hmi.hmihet, hmi.HMinus_induc_rec_cooling, hmi.hmicol );
06231 }
06232
06233
06234 if( trace.lgTrace && trace.lgTr_H2_Mole )
06235 {
06236 rate = desh2p;
06237 if( rate != 0. )
06238 {
06239 fprintf( ioQQQ,
06240 " 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",
06241 rate, h2pcin*dense.eden/rate, gamtwo/rate, 2.*iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/
06242 rate, hmi.h2ph3p/rate, h2phhp/rate, h2pion/rate, hmi.bh2h2p*
06243 dense.xIonDense[ipHYDROGEN][0]/rate );
06244
06245 rate *= hmi.Hmolec[ipMH2p];
06246 if( rate > 0. )
06247 {
06248 fprintf( ioQQQ,
06249 " Create H2+: rate=%.2e HII HI;%.3f Col H2;%.3f HII H2;%.3f HI HI;%.3f\n",
06250 rate,
06251 radath*dense.xIonDense[ipHYDROGEN][1]*dense.xIonDense[ipHYDROGEN][0]/rate,
06252 (hmi.H2_photoionize_rate + secondaries.csupra[ipHYDROGEN][0]*2.02)*hmi.Hmolec[ipMH2g]/rate,
06253 hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1]*hmi.Hmolec[ipMH2g]/rate,
06254 b2pcin*dense.xIonDense[ipHYDROGEN][0]*dense.xIonDense[ipHYDROGEN][1]*dense.eden/rate );
06255 }
06256 else
06257 {
06258 fprintf( ioQQQ, " Create H2+: rate= is zero\n" );
06259 }
06260 }
06261 }
06262
06263 {
06264
06265 enum {DEBUG_LOC=false};
06266 if( DEBUG_LOC )
06267 {
06268 fprintf(ioQQQ,"hmole bugg\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
06269 fnzone,
06270 iso.gamnc[ipH_LIKE][ipHYDROGEN][0],
06271 hmi.HMinus_photo_rate,
06272 hmi.Hmolec[ipMH2g] ,
06273 hmi.Hmolec[ipMHm] ,
06274 dense.xIonDense[ipHYDROGEN][1]);
06275 }
06276 }
06277 return;
06278 }
06279 #if defined(__HP_aCC)
06280 #pragma OPTIMIZE OFF
06281 #pragma OPTIMIZE ON
06282 #endif
06283
06284