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