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