00001
00002
00003
00004 #include "cddefines.h"
00005 #include "taulines.h"
00006 #include "dense.h"
00007 #include "hmi.h"
00008 #include "conv.h"
00009 #include "doppvel.h"
00010 #include "thermal.h"
00011 #include "gammas.h"
00012 #include "iso.h"
00013 #include "opacity.h"
00014 #include "phycon.h"
00015 #include "physconst.h"
00016 #include "radius.h"
00017 #include "rfield.h"
00018 #include "secondaries.h"
00019 #include "timesc.h"
00020 #include "trace.h"
00021 #include "thermal.h"
00022 #include "called.h"
00023 #include "hcmap.h"
00024 #include "coolheavy.h"
00025 #include "mole.h"
00026 #include "mole_priv.h"
00027 #include "dynamics.h"
00028 #include "h2.h"
00029 #include "co.h"
00030 #include "ionbal.h"
00031
00032
00033 static const double MOLETOLER = 0.10;
00034
00035 STATIC void mole_effects(void);
00036 STATIC void mole_h_rate_diagnostics(void);
00037 STATIC void mole_ion_trim(void);
00038 STATIC void mole_update_limiting_reactants();
00039
00040
00041 void mole_drive(void)
00042 {
00043 DEBUG_ENTRY( "mole_drive()" );
00044
00045 mole_update_species_cache();
00046
00047 mole_update_limiting_reactants();
00048
00049 mole_update_rks();
00050
00051 double error = mole_solve();
00052
00053 bool lgConverged = (error < MOLETOLER);
00054
00055 if( !lgConverged )
00056 {
00057
00058 }
00059
00060 mole_ion_trim();
00061
00062 mole_effects();
00063
00064 return;
00065 }
00066
00067 STATIC void mole_ion_trim(void)
00068 {
00069 for (ChemAtomList::iterator atom = unresolved_atom_list.begin();
00070 atom != unresolved_atom_list.end(); ++atom)
00071 {
00072 const long int nelem=(*atom)->el->Z-1;
00073 if( !dense.lgElmtOn[nelem] )
00074 continue;
00075
00076 for (long int ion=0;ion<nelem+2;ion++)
00077 {
00078 if ((*atom)->ipMl[ion] != -1)
00079 {
00080 if (dense.IonLow[nelem] > ion )
00081 {
00082 if( dense.xIonDense[nelem][ion] > ( ionbal.trimlo * dense.gas_phase[nelem] ) )
00083 dense.IonLow[nelem] = ion;
00084 else
00085 {
00086 dense.xIonDense[nelem][dense.IonLow[nelem]] += dense.xIonDense[nelem][ion];
00087 dense.xIonDense[nelem][ion] = 0.;
00088 }
00089 }
00090
00091 if (dense.IonHigh[nelem] < ion )
00092 {
00093 if( dense.xIonDense[nelem][ion] > ( ionbal.trimhi * dense.gas_phase[nelem] ) )
00094 dense.IonHigh[nelem] = ion;
00095 else
00096 {
00097 dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion];
00098 dense.xIonDense[nelem][ion] = 0.;
00099 }
00100 }
00101 }
00102 }
00103 }
00104 }
00105
00106 void mole_update_sources(void)
00107 {
00108 DEBUG_ENTRY( "mole_update_sources()" );
00109
00110 mole_update_species_cache();
00111 mole_eval_sources(mole_global.num_total);
00112 }
00113
00114 STATIC void mole_effects()
00115 {
00116 double
00117 plte;
00118
00119 DEBUG_ENTRY( "mole_effects()" );
00120
00121 mole_eval_sources(mole_global.num_total);
00122
00123 co.CODissHeat = (realnum)(mole.findrate("PHOTON,CO=>C,O")*1e-12);
00124
00125 thermal.heating[0][9] = co.CODissHeat;
00126
00127
00128 hmi.H2_total = findspecieslocal("H2*")->den + findspecieslocal("H2")->den;
00129 hmi.HD_total = findspecieslocal("HD")->den;
00130
00131 h2.ortho_density = 0.75*hmi.H2_total;
00132 h2.para_density = 0.25*hmi.H2_total;
00133 {
00134 hmi.H2_total_f = (realnum)hmi.H2_total;
00135 h2.ortho_density_f = (realnum)h2.ortho_density;
00136 h2.para_density_f = (realnum)h2.para_density;
00137 }
00138
00139
00140
00141
00142
00143
00144 hmi.H2_rate_destroy = mole.sink_rate_tot("H2");
00145
00146
00147 if(hmi.H2_rate_destroy > SMALLFLOAT )
00148 {
00149
00150 timesc.time_H2_Dest_here = 1./hmi.H2_rate_destroy;
00151 }
00152 else
00153 {
00154 timesc.time_H2_Dest_here = 0.;
00155 }
00156
00157
00158
00159
00160
00161 timesc.time_H2_Form_here = (mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn"));
00162
00163 if( timesc.time_H2_Form_here > SMALLFLOAT )
00164 {
00165
00166 timesc.time_H2_Form_here = hmi.H2_total/timesc.time_H2_Form_here;
00167 }
00168 else
00169 {
00170 timesc.time_H2_Form_here = 0.;
00171 }
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 hmi.HeatH2Dish_TH85 = 0.25 * EN1EV *
00200 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
00201 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
00202
00203
00204 hmi.HeatH2Dexc_TH85 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
00205 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
00206
00207 hmi.deriv_HeatH2Dexc_TH85 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_TH85)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
00208
00209 if( hmi.chH2_small_model_type == 'H' )
00210 {
00211
00212 hmi.HeatH2Dish_BHT90 = 0.25 * EN1EV *
00213 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
00214 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
00215
00216
00217 hmi.HeatH2Dexc_BHT90 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
00218 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
00219
00220 hmi.deriv_HeatH2Dexc_BHT90 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BHT90)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
00221 }
00222 else if( hmi.chH2_small_model_type == 'B')
00223 {
00224
00225 hmi.HeatH2Dish_BD96 = 0.25 * EN1EV *
00226 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
00227 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
00228
00229 hmi.HeatH2Dexc_BD96 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
00230 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
00231
00232 hmi.deriv_HeatH2Dexc_BD96 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BD96)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
00233 }
00234 else if(hmi.chH2_small_model_type == 'E')
00235 {
00236
00237
00238
00239
00240
00241 double log_density,
00242 f1, f2,f3, f4, f5;
00243 static double log_G0_face = -1;
00244 static double k_f4;
00245
00246
00247
00248
00249 if( !nzone )
00250 {
00251 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
00252 {
00253 log_G0_face = 0.;
00254 }
00255 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
00256 {
00257 log_G0_face = 7.;
00258 }
00259 else
00260 {
00261 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
00262 }
00263
00264 log_G0_face /= radius.r1r0sq;
00265 }
00266
00267
00268 if(dense.gas_phase[ipHYDROGEN] <= 1.)
00269 {
00270 log_density = 0.;
00271 }
00272 else if(dense.gas_phase[ipHYDROGEN] >= 1e7)
00273 {
00274 log_density = 7.;
00275 }
00276 else
00277 {
00278 log_density = log10(dense.gas_phase[ipHYDROGEN]);
00279 }
00280
00281 f1 = 0.15 * log_density + 0.75;
00282 f2 = -0.5 * log_density + 10.;
00283
00284 hmi.HeatH2Dish_ELWERT = 0.25 * EN1EV * f1 *
00285 (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
00286 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den ) +
00287 f2 * secondaries.x12tot * EN1EV * hmi.H2_total;
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
00299 {
00300 log_G0_face = 0.;
00301 }
00302 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
00303 {
00304 log_G0_face = 7.;
00305 }
00306 else
00307 {
00308 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
00309 }
00310
00311 log_G0_face /= radius.r1r0sq;
00312
00313
00314 k_f4 = -0.25 * log_G0_face + 1.25;
00315
00316
00317 if(dense.gas_phase[ipHYDROGEN] <= 1.)
00318 {
00319 log_density = 0.;
00320 f4 = 0.;
00321 }
00322 else if(dense.gas_phase[ipHYDROGEN] >= 1e7)
00323 {
00324 log_density = 7.;
00325 f4 = pow2(k_f4) * pow( 10. , 2.2211 * log_density - 29.8506);
00326 }
00327 else
00328 {
00329 log_density = log10(dense.gas_phase[ipHYDROGEN]);
00330 f4 = pow2(k_f4) * pow( 10., 2.2211 * log_density - 29.8506);
00331 }
00332
00333 f3 = MAX2(0.1, -4.75 * log_density + 24.25);
00334 f5 = MAX2(1.,0.95 * log_density - 1.45) * 0.2 * log_G0_face;
00335
00336 hmi.HeatH2Dexc_ELWERT = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
00337 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12 * f3 +
00338 f4 * (mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[0]].den/dense.gas_phase[ipHYDROGEN]) +
00339 f5 * secondaries.x12tot * EN1EV * hmi.H2_total;
00340
00341 if(log_G0_face == 0.&& dense.gas_phase[ipHYDROGEN] > 1.)
00342 hmi.HeatH2Dexc_ELWERT *= 0.001 / dense.gas_phase[ipHYDROGEN];
00343
00344
00345
00346 hmi.HeatH2Dexc_ELWERT /= POW2(radius.r1r0sq);
00347
00348
00349 hmi.deriv_HeatH2Dexc_ELWERT = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_ELWERT)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
00350
00351
00352 }
00353
00354 else
00355 TotalInsanity();
00356
00357
00358 if( findspecieslocal("H-")->den > 0. && hmi.rel_pop_LTE_Hmin > 0. )
00359 {
00360 hmi.hmidep = (double)findspecieslocal("H-")->den/ SDIV(
00361 ((double)dense.xIonDense[ipHYDROGEN][0]*dense.eden*hmi.rel_pop_LTE_Hmin));
00362 }
00363 else
00364 {
00365 hmi.hmidep = 1.;
00366 }
00367
00368
00369 hmi.hmihet = hmi.HMinus_photo_heat*findspecieslocal("H-")->den - hmi.HMinus_induc_rec_cooling;
00370
00371
00372 t_phoHeat photoHeat;
00373 (void) GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.,&photoHeat);
00374 hmi.h2plus_heat = photoHeat.HeatNet*findspecieslocal("H2+")->den;
00375
00376
00377
00378 plte = (double)dense.xIonDense[ipHYDROGEN][0] * hmi.rel_pop_LTE_H2g * (double)dense.xIonDense[ipHYDROGEN][0];
00379 if( plte > 0. )
00380 {
00381 hmi.h2dep = findspecieslocal("H2")->den/plte;
00382 }
00383 else
00384 {
00385 hmi.h2dep = 1.;
00386 }
00387
00388
00389
00390 plte = (double)dense.xIonDense[ipHYDROGEN][0]*hmi.rel_pop_LTE_H2p*(double)dense.xIonDense[ipHYDROGEN][1];
00391 if( plte > 0. )
00392 {
00393 hmi.h2pdep = findspecieslocal("H2+")->den/plte;
00394 }
00395 else
00396 {
00397 hmi.h2pdep = 1.;
00398 }
00399
00400
00401 if( hmi.rel_pop_LTE_H3p > 0. )
00402 {
00403 hmi.h3pdep = findspecieslocal("H3+")->den/hmi.rel_pop_LTE_H3p;
00404 }
00405 else
00406 {
00407 hmi.h3pdep = 1.;
00408 }
00409
00410 mole_h_rate_diagnostics();
00411
00412 return;
00413 }
00414 #if defined(__HP_aCC)
00415 #pragma OPTIMIZE OFF
00416 #pragma OPTIMIZE ON
00417 #endif
00418 STATIC void mole_h_rate_diagnostics(void)
00419 {
00420 int i;
00421 double
00422 rate;
00423
00424
00425 {
00426
00427
00428 enum {DEBUG_LOC=false};
00429
00430 if( DEBUG_LOC && (nzone>50) )
00431 {
00432 double createsum ,create_from_Hn2 , create_3body_Ho, create_h2p,
00433 create_h3p, create_h3pe, create_grains, create_hminus;
00434 double destroysum, destroy_hm ,destroy_solomon ,destroy_2h ,destroy_hp,
00435 destroy_h,destroy_hp2,destroy_h3p;
00436
00437 create_from_Hn2 = mole.findrate("H,H=>H2");
00438 create_3body_Ho = mole.findrate("H,H,H=>H2,H");
00439 create_h2p = mole.findrate("H,H2+=>H+,H2*");
00440 create_h3p = mole.findrate("H,H3+=>H2,H2+");
00441 create_h3pe = mole.findrate("e-,H3+=>H2,H");
00442 create_grains = mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn");
00443 create_hminus = (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"));
00444
00445 createsum = create_from_Hn2 + create_3body_Ho + create_h2p +
00446 create_h3p + create_h3pe + create_grains + create_hminus;
00447
00448 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",
00449 fnzone,
00450 createsum,
00451 create_hminus / createsum,
00452 create_from_Hn2 / createsum,
00453 create_3body_Ho / createsum,
00454 create_h2p / createsum,
00455 create_h3p / createsum,
00456 create_h3pe / createsum,
00457 create_grains / createsum );
00458
00459
00460
00461
00462 destroy_hm = (mole.findrate("H2,e-=>H,H-")+mole.findrate("H2*,e-=>H,H-"));
00463
00464 destroy_solomon = hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den;
00465 destroy_2h = (mole.findrate("H2,e-=>H,H,e-")+mole.findrate("H2*,e-=>H,H,e-"));
00466 destroy_hp = mole.findrate("H2,H+=>H3+");
00467 destroy_h = mole.findrate("H,H2=>H,H,H");
00468 destroy_hp2 = mole.findrate("H2,H+=>H2+,H");
00469 destroy_h3p = mole.findrate("H2,H3+=>H2,H2+,H");
00470 destroysum = destroy_hm + destroy_solomon + destroy_2h +
00471 destroy_hp+ destroy_h+ destroy_hp2+ destroy_h3p;
00472
00473 fprintf(ioQQQ,"H2 destroy\t%.3f \t%.2e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00474 fnzone,
00475 destroysum,
00476 destroy_hm / destroysum ,
00477 destroy_solomon / destroysum ,
00478 destroy_2h / destroysum ,
00479 destroy_hp / destroysum ,
00480 destroy_h / destroysum ,
00481 destroy_hp2 / destroysum ,
00482 destroy_h3p / destroysum );
00483
00484 }
00485 }
00486
00487 {
00488
00489
00490 enum {DEBUG_LOC=false};
00491
00492 if( DEBUG_LOC && (nzone>140) )
00493 {
00494 double create_from_Ho,create_3body_Ho,create_batach,destroy_photo,
00495 destroy_coll_heavies,destroy_coll_electrons,destroy_Hattach,destroy_fhneut,
00496 destsum , createsum;
00497
00498 create_from_Ho = mole.findrate("H,e-=>H-,PHOTON");
00499 create_3body_Ho = mole.findrate("H,e-,e-=>H-,e-");
00500
00501 create_batach = (mole.findrate("H2,e-=>H,H-") + mole.findrate("H2*,e-=>H,H-"));
00502 destroy_photo = mole.findrate("H-,PHOTON=>H,e-");
00503 destroy_coll_heavies = 0.;
00504 for( i=ipLITHIUM; i < LIMELM; i++ )
00505 {
00506 char react[32];
00507 if (dense.lgElmtOn[i])
00508 {
00509 sprintf(react,"H-,%s=>H,%s",mole_global.list[unresolved_atom_list[i]->ipMl[1]]->label.c_str(),unresolved_atom_list[i]->label().c_str() );
00510 destroy_coll_heavies += mole.findrate(react);
00511 }
00512 }
00513 destroy_coll_electrons = mole.findrate("H-,e-=>H-,e-,e-");
00514 destroy_Hattach = (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"));
00515 destroy_fhneut = mole.findrate("H-,H+=>H,H");
00516
00517 destsum = destroy_photo + destroy_coll_heavies + destroy_coll_electrons +
00518 destroy_Hattach + destroy_fhneut;
00519 fprintf(ioQQQ,"H- destroy zone\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n",
00520 fnzone,
00521 phycon.te,
00522 destsum,
00523 destroy_photo/destsum ,
00524 destroy_coll_heavies/destsum,
00525 destroy_coll_electrons/destsum,
00526 destroy_Hattach/destsum,
00527 destroy_fhneut/destsum );
00528
00529 createsum = create_from_Ho+create_3body_Ho+create_batach;
00530 fprintf(ioQQQ,"H- create\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00531 fnzone,
00532 phycon.te,
00533 createsum,
00534 dense.eden,
00535 create_from_Ho/createsum,
00536 create_3body_Ho/createsum,
00537 create_batach/createsum);
00538 }
00539 }
00540 {
00541
00542
00543 enum {DEBUG_LOC=false};
00544
00545 if( DEBUG_LOC )
00546 {
00547
00548 fprintf( ioQQQ, " MOLE raw; hi\t%.2e" , dense.xIonDense[ipHYDROGEN][0]);
00549 for( i=0; i < mole_global.num_calc; i++ )
00550 {
00551 fprintf( ioQQQ, "\t%-4.4s\t%.2e", mole_global.list[i]->label.c_str(), mole.species[i].den );
00552 }
00553 fprintf( ioQQQ, " \n" );
00554 }
00555 }
00556
00557 if( trace.lgTrace && trace.lgTraceMole )
00558 {
00559
00560 fprintf( ioQQQ, " raw; " );
00561 for( i=0; i < mole_global.num_calc; i++ )
00562 {
00563 fprintf( ioQQQ, " %-4.4s:%.2e", mole_global.list[i]->label.c_str(), mole.species[i].den );
00564 }
00565 fprintf( ioQQQ, " \n" );
00566
00567 }
00568
00569
00570
00571
00572 if( (trace.lgTrace && trace.lgTraceMole) )
00573 {
00575 double H2_rate_create = mole.source_rate_tot("H2") + mole.source_rate_tot("H2*");
00576
00577 if( H2_rate_create > SMALLFLOAT )
00578 {
00579 fprintf( ioQQQ, " 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",
00580 H2_rate_create,
00581 (mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn"))/H2_rate_create,
00582 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"))/H2_rate_create,
00583 mole.findrate("H,H,H=>H2,H")/H2_rate_create,
00584 mole.findrate("H,H2+=>H+,H2*")/H2_rate_create,
00585 mole.findrate("H,H=>H2")/H2_rate_create,
00586 mole.findrate("H,H3+=>H2,H2+")/H2_rate_create,
00587 mole.findrate("H2,H3+=>H2,H2+,H")/H2_rate_create );
00588 }
00589 else
00590 {
00591 fprintf( ioQQQ, " Create H2, rate=0\n" );
00592 }
00593 }
00594
00595
00596 if( trace.lgTrace && trace.lgTraceMole )
00597 {
00598 rate = mole.findrate("H2,H+=>H2+,H") + mole.findrate("H,H+,e-=>H2+,e-") +
00599 mole.findrate("H,H3+=>H2,H2+") + mole.findrate("H2,H3+=>H2,H2+,H");
00600 if( rate > 1e-25 )
00601 {
00602 fprintf( ioQQQ, " Create H2+, rate=%10.2e hmi.rh2h2p;%5.3f b2pcin;%5.3f hmi.h3ph2p;%5.3f hmi.h3petc+;%5.3f\n",
00603 rate, mole.findrate("H2,H+=>H2+,H")/rate,
00604 mole.findrate("H,H+,e-=>H2+,e-")/rate, mole.findrate("H,H3+=>H2,H2+")/rate, mole.findrate("H2,H3+=>H2,H2+,H")/rate );
00605 }
00606 else
00607 {
00608 fprintf( ioQQQ, " Create H2+, rate=0\n" );
00609 }
00610 }
00611
00612 if( trace.lgTrace && trace.lgTraceMole )
00613 {
00614
00615 double destroy_coll_heavies = 0.;
00616
00617 for( i=ipLITHIUM; i < LIMELM; i++ )
00618 {
00619 char react[32];
00620 if (dense.lgElmtOn[i])
00621 {
00622 sprintf(react,"H-,%s=>H,%s",mole_global.list[unresolved_atom_list[i]->ipMl[1]]->label.c_str(),unresolved_atom_list[i]->label().c_str() );
00623 destroy_coll_heavies += mole.findrate(react);
00624 }
00625 }
00626 fprintf( ioQQQ, " MOLE, Dep Coef, H-:%10.2e H2:%10.2e H2+:%10.2e\n",
00627 hmi.hmidep, hmi.h2dep, hmi.h2pdep );
00628 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",
00629 mole.findrate("H,e-=>H-,PHOTON"), hmi.HMinus_induc_rec_rate*findspecieslocal("H-")->den, mole.findrate("H,H=>H-,H+"),
00630 mole.findrate("H,e-,e-=>H-,e-"),
00631 mole.findrate("H2,e-=>H,H-"), findspecieslocal("H-")->den, hmi.hmidep );
00632
00633 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",
00634 mole.findrate("H-,PHOTON=>H,e-"), destroy_coll_heavies,
00635 mole.findrate("H-,e-=>H-,e-,e-"),
00636 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")), iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*findspecieslocal("H-")->den,
00637 mole.findrate("H-,H+=>H,H") );
00638 fprintf( ioQQQ, " H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n",
00639 hmi.hmihet, hmi.HMinus_induc_rec_cooling, hmi.hmicol );
00640 }
00641
00642
00643 if( trace.lgTrace && trace.lgTraceMole )
00644 {
00645 rate = findspecieslocal("H2+")->snk;
00646 if( rate != 0. )
00647 {
00648 rate *= findspecieslocal("H2+")->den;
00649 if( rate > 0. )
00650 {
00651 fprintf( ioQQQ,
00652 " 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",
00653 rate, mole.findrate("H2+,e-=>H,H+,e-")/rate, mole.findrate("H2+,PHOTON=>H,H+")/rate,
00654 mole.findrate("H2+,CRPHOT=>H,H+")/rate, mole.findrate("H2,H2+=>H,H3+")/rate, mole.findrate("H2+,H2=>H,H+,H2")/rate,
00655 mole.findrate("H2+,H+=>H,H+,H+")/rate, mole.findrate("H,H2+=>H+,H2*")/rate );
00656
00657 fprintf( ioQQQ,
00658 " Create H2+: rate=%.2e HII HI;%.3f Col H2;%.3f HII H2;%.3f HI HI;%.3f\n",
00659 rate,
00660 mole.findrate("H+,H=>H2+,PHOTON")/rate,
00661 mole.findrate("H2,CRPHOT=>H2+,e-")/rate,
00662 mole.findrate("H2,H+=>H2+,H")/rate,
00663 mole.findrate("H,H+,e-=>H2+,e-")/rate );
00664 }
00665 else
00666 {
00667 fprintf( ioQQQ, " Create H2+: rate= is zero\n" );
00668 }
00669 }
00670 }
00671
00672 {
00673
00674
00675 enum {DEBUG_LOC=false};
00676
00677 if( DEBUG_LOC )
00678 {
00679 fprintf(ioQQQ,"mole bugg\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00680 fnzone,
00681 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].gamnc,
00682 hmi.HMinus_photo_rate,
00683 findspecieslocal("H2")->den ,
00684 findspecieslocal("H-")->den ,
00685 findspecieslocal("H+")->den);
00686 }
00687 }
00688
00689 {
00690
00691
00692 enum {DEBUG_LOC=false};
00693
00694 if( DEBUG_LOC && nzone>140 )
00695 {
00696 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",
00697 fnzone ,
00698 mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn") ,
00699 hmi.H2_forms_grains+hmi.H2star_forms_grains ,
00700 mole.findrate("H,H,grn=>H2*,grn")/SDIV(mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn")),
00701 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),
00702 hmi.H2star_forms_hminus+hmi.H2_forms_hminus,
00703 frac_H2star_hminus(),
00704 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),findspecieslocal("H")->den,findspecieslocal("H-")->den
00705 );
00706 }
00707 }
00708
00709 {
00710
00711
00712
00713
00714 enum {DEBUG_LOC=false};
00715
00716 if( DEBUG_LOC && nzone>140)
00717 {
00718
00719 fprintf(ioQQQ,
00720 "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",
00721
00722 fnzone,
00723 phycon.te,
00724 dense.eden,
00725 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),
00726 mole.findrate("H2,e-=>H,H-"),
00727 mole.findrate("H2*,e-=>H,H-"),
00728 mole.findrate("H,e-=>H-,PHOTON"),
00729 hmi.HMinus_induc_rec_rate*findspecieslocal("H-")->den,
00730 mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[0]].den,
00731 mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[1]].den,
00732 findspecieslocal("H-")->den,
00733 hmi.H2_total,
00734 hmi.rel_pop_LTE_Hmin,
00735 hmi.rel_pop_LTE_H2g,
00736 hmi.rel_pop_LTE_H2s
00737 );
00738 }
00739 }
00740
00741 if( (trace.lgTrace && trace.lgTraceMole) )
00742 {
00743 if( hmi.H2_rate_destroy > SMALLFLOAT )
00744 {
00745 fprintf( ioQQQ,
00746 " H2 destroy rate=%.2e DIS;%.3f bat;%.3f h2dis;%.3f photoionize_rate;%.3f h2h2p;%.3f E-h;%.3f hmi.h2hph3p;%.3f sec;%.3f\n",
00747 hmi.H2_rate_destroy,
00748 hmi.H2_Solomon_dissoc_rate_used_H2g / hmi.H2_rate_destroy,
00749 mole.findrate("H2,e-=>H,H-") / (hmi.H2_total*hmi.H2_rate_destroy),
00750 mole.findrate("H,H2=>H,H,H") / (hmi.H2_total*hmi.H2_rate_destroy),
00751 h2.photoionize_rate / hmi.H2_rate_destroy,
00752 mole.findrate("H2,H+=>H2+,H") / (hmi.H2_total* hmi.H2_rate_destroy),
00753 (mole.findrate("H2,e-=>H,H,e-")+mole.findrate("H2*,e-=>H,H,e-")) / (hmi.H2_total* hmi.H2_rate_destroy),
00754 mole.findrate("H2,H+=>H3+") / (hmi.H2_total* hmi.H2_rate_destroy) ,
00755 secondaries.csupra[ipHYDROGEN][0]*2.02 / hmi.H2_rate_destroy
00756 );
00757 }
00758 else
00759 {
00760 fprintf( ioQQQ, " Destroy H2: rate=0\n" );
00761 }
00762 }
00763
00764 {
00765
00766
00767 enum {DEBUG_LOC=false};
00768
00769 if( DEBUG_LOC )
00770 {
00771 if( DEBUG_LOC && (nzone > 570) )
00772 {
00773 fprintf(ioQQQ,"Temperature %g\n",phycon.te);
00774 fprintf(ioQQQ," Net mol ion rate [%g %g] %g\n",mole.source[ipHYDROGEN][1],mole.sink[ipHYDROGEN][1],
00775 mole.source[ipHYDROGEN][1]-mole.sink[ipHYDROGEN][1]*mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[1]].den);
00776 }
00777 }
00778 }
00779 }
00780
00781 STATIC void mole_update_limiting_reactants()
00782 {
00783 DEBUG_ENTRY( "mole_update_limiting_reactants()" );
00784
00785
00786 for( long i=0; i<mole_global.num_calc; ++i )
00787 {
00788 mole.species[i].xFracLim = 0.;
00789 mole.species[i].nAtomLim = -1;
00790 for (molecule::nAtomsMap::iterator atom = mole_global.list[i]->nAtom.begin();
00791 atom != mole_global.list[i]->nAtom.end(); ++atom)
00792 {
00793 long nelem = atom->first->el->Z-1;
00794 if( dense.lgElmtOn[nelem])
00795 {
00796 realnum dens_elemsp = (realnum)( mole.species[i].den * atom->second * atom->first->frac );
00797 if (mole.species[i].xFracLim*dense.gas_phase[nelem] < dens_elemsp )
00798 {
00799 mole.species[i].xFracLim = dens_elemsp/dense.gas_phase[nelem];
00800 mole.species[i].nAtomLim = nelem;
00801 }
00802 }
00803 }
00804
00805 }
00806
00807 return;
00808 }
00809