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