00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "physconst.h"
00009 #include "dense.h"
00010 #include "called.h"
00011 #include "gammas.h"
00012 #include "colden.h"
00013 #include "thermal.h"
00014 #include "secondaries.h"
00015 #include "h2.h"
00016 #include "mole.h"
00017 #include "radius.h"
00018 #include "doppvel.h"
00019 #include "rfield.h"
00020 #include "ionbal.h"
00021 #include "rt.h"
00022 #include "opacity.h"
00023 #include "iso.h"
00024 #include "conv.h"
00025 #include "phycon.h"
00026 #include "hmi.h"
00027
00028
00029 #if 0
00030 #if !defined(NDEBUG)
00031
00032
00033
00034
00035 #define AUDIT(a) { \
00036 double total, mtotal; \
00037 for( i=0;i<N_H_MOLEC;i++) { \
00038 total = 0.; \
00039 for( j=0;j<N_H_MOLEC;j++) { \
00040 total += c[i][j]*nprot[j]; \
00041 } \
00042 if( fabs(total) > 1e-6*fabs(c[i][i]*nprot[i])) { \
00043 fprintf(ioQQQ,"PROBLEM Subtotal1 %c %.2e\n",a,fabs(total)/fabs(c[i][i]*nprot[i])); \
00044 fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*nprot[i]); \
00045 } \
00046 } \
00047 total = mtotal = 0.;for( j=0;j<N_H_MOLEC;j++) { total += bvec[j]*nprot[j]; mtotal += fabs(bvec[j]*nprot[j]); }\
00048 if( fabs(total) > 1e-30 && fabs(total) > 1e-10*rtot) { \
00049 fprintf(ioQQQ,"PROBLEM Subtotal2 %c %.2e\n",a,fabs(total)/mtotal); \
00050 fprintf(ioQQQ,"RHS Total %g cf %g\n",total,mtotal); \
00051 } else if( a == '.' && fabs(total) > 1e-7*mtotal) { \
00052 fprintf(ioQQQ,"PROBLEM zone %li Hmole RHS conservation error %.2e of %.2e\n",nzone,total,mtotal); \
00053 fprintf(ioQQQ,"(may be due to high rate equilibrium reactions)\n"); \
00054 } \
00055 }
00056 #else
00057 #define AUDIT
00058 #endif
00059
00060 #endif
00061
00062 void hmole( void )
00063 {
00064 int nFixup, i;
00065 double error;
00066 float oxy_error=0.;
00067 static float abund0=-BIGFLOAT , abund1=-BIGFLOAT;
00068 float save1=dense.xIonDense[ipOXYGEN][1],
00069 save0=dense.xIonDense[ipOXYGEN][0];
00070
00071 DEBUG_ENTRY( "hmirat()" );
00072
00073
00074 nFixup = 0;
00075 error = 1.;
00076
00077
00078
00079 if( conv.nTotalIoniz==0 && iteration==0 )
00080 {
00081 float pdr_mole_h2[N_H_MOLEC] = {1.00E+00f,
00082 3.18E-05f,
00083 1.95E-11f,
00084 4.00E-08f,
00085 1.08E-14f,
00086 1.08E-20f,
00087 3.85E-07f,
00088 8.04E-14f};
00089
00090
00091 ASSERT( dense.xIonDense[ipHYDROGEN][0]>SMALLFLOAT );
00092
00093 ASSERT( ipMH==0&&
00094 ipMHp == 1&&
00095 ipMHm == 2&&
00096 ipMH2g == 3&&
00097 ipMH2p == 4&&
00098 ipMH3p == 5&&
00099 ipMH2s == 6&&
00100 ipMHeHp== 7&&
00101 N_H_MOLEC==8 );
00102
00103 for( i=0; i<N_H_MOLEC; ++i )
00104 {
00105 hmi.Hmolec[i] = dense.xIonDense[ipHYDROGEN][0]*pdr_mole_h2[i];
00106 }
00107 }
00108
00109
00110 hmole_reactions();
00111
00112
00113 # define BIGERROR 1e-4
00114
00115
00116 # define LIM_LOOP 20
00117
00118 for(i=0; i < LIM_LOOP && ((error > BIGERROR||nFixup || oxy_error>conv.EdenErrorAllowed));i++)
00119 {
00120 {
00121
00122
00123 enum {DEBUG_LOC=false};
00124
00125 if( DEBUG_LOC && (nzone>140) )
00126 {
00127 fprintf(ioQQQ,"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
00128 fnzone,
00129 phycon.te,
00130 dense.eden);
00131 for( i=0; i<N_H_MOLEC; i++ )
00132 fprintf(ioQQQ,"\t%.2e", hmi.Hmolec[i] );
00133 fprintf(ioQQQ,"\n" );
00134 }
00135 }
00136
00137
00138 nFixup = 0;
00139
00140 if( !conv.lgSearch )
00141 {
00142 IonOxyge();
00143 }
00144
00145 save1 = dense.xIonDense[ipOXYGEN][1];
00146 save0 = dense.xIonDense[ipOXYGEN][0];
00147
00148
00149
00150
00151 if( nzone )
00152 {
00153 # define OLD 0.75f
00154 abund0 = abund0*OLD+dense.xIonDense[ipOXYGEN][0]*(1.f-OLD);
00155 abund1 = abund1*OLD+dense.xIonDense[ipOXYGEN][1]*(1.f-OLD);
00156 }
00157 else
00158 {
00159 abund0 = dense.xIonDense[ipOXYGEN][0];
00160 abund1 = dense.xIonDense[ipOXYGEN][1];
00161 }
00162 dense.xIonDense[ipOXYGEN][0] = abund0;
00163
00164 dense.xIonDense[ipOXYGEN][1] = abund1;
00165
00166 oxy_error = (float)fabs(save1-abund1)/SDIV(dense.gas_phase[ipOXYGEN]);
00167
00168
00169
00170 hmole_step(&nFixup, &error);
00171 dense.xIonDense[ipOXYGEN][0] = save0;
00172 dense.xIonDense[ipOXYGEN][1] = save1;
00173 {
00174
00175
00176 enum {DEBUG_LOC=false};
00177
00178 if( DEBUG_LOC )
00179 {
00180 fprintf(ioQQQ,"DEBUG hmole out\t%i\t%.2f\t%.5e\t%.5e",
00181 i,
00182 fnzone,
00183 phycon.te,
00184 dense.eden);
00185 fprintf(ioQQQ,
00186 "\terror\t%.4e\tH0\t%.4e\tH+\t%.4e\tsink\t%.4e\t%.4e\tsour\t%.4e\t%.4e\tion\t%.4e\trec\t%.4e",
00187 error,
00188 dense.xIonDense[ipHYDROGEN][0],
00189 dense.xIonDense[ipHYDROGEN][1],
00190 mole.sink[ipHYDROGEN][0],
00191 mole.sink[ipHYDROGEN][1],
00192 mole.source[ipHYDROGEN][0] ,
00193 mole.source[ipHYDROGEN][1] ,
00194 ionbal.RateIonizTot[ipHYDROGEN][0],
00195 ionbal.RateRecomTot[ipHYDROGEN][0]);
00196
00197
00198
00199 fprintf(ioQQQ,"\n" );
00200 }
00201 }
00202 }
00203
00204 if( (i == LIM_LOOP && error > BIGERROR) || nFixup != 0 )
00205 {
00206
00207
00208 conv.lgConvPops = false;
00209
00210 if( !conv.lgSearch && called.lgTalk )
00211 {
00212 fprintf(ioQQQ," PROBLEM hmole, zone %li: %d iters, %d bad; final error: %g lgSearch %i\n",
00213 nzone,
00214 i,
00215 nFixup,
00216 error,
00217 conv.lgSearch);
00218 }
00219 ConvFail( "pops" , "Hmole");
00220 }
00221
00222
00223 if( strcmp( dense.chDenseLaw, "CDEN" )==0 )
00224
00225 ASSERT( fabs( dense.gas_phase[ipHYDROGEN] - dense.den0 )/
00226 dense.gas_phase[ipHYDROGEN]<1e-4 );
00227
00228
00229
00230
00231 dense.xMolecules[ipHYDROGEN] = 0.;
00232 for(i=0;i<N_H_MOLEC;i++)
00233 {
00234 dense.xMolecules[ipHYDROGEN] += hmi.Hmolec[i]*hmi.nProton[i];
00235 }
00236
00237 dense.xMolecules[ipHYDROGEN] -= (hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp]);
00238
00239
00240 for( i=0; i < mole.num_comole_calc; i++ )
00241 {
00242 dense.xMolecules[ipHYDROGEN] += COmole[i]->hevmol*COmole[i]->nElem[ipHYDROGEN];
00243 }
00244
00245
00246 DEBUG_EXIT( "hmole()" );
00247 return;
00248 }
00249
00250
00251 double hmirat(double te)
00252 {
00253 double hmirat_v;
00254
00255 DEBUG_ENTRY( "hmirat()" );
00256
00257
00258 if( te < 31.62 )
00259 {
00260 hmirat_v = 8.934e-18*phycon.sqrte*phycon.te003*phycon.te001*
00261 phycon.te001;
00262 }
00263 else if( te < 90. )
00264 {
00265 hmirat_v = 5.159e-18*phycon.sqrte*phycon.te10*phycon.te03*
00266 phycon.te03*phycon.te003*phycon.te001;
00267 }
00268 else if( te < 1200. )
00269 {
00270 hmirat_v = 2.042e-18*te/phycon.te10/phycon.te03;
00271 }
00272 else if( te < 3800. )
00273 {
00274 hmirat_v = 8.861e-18*phycon.te70/phycon.te03/phycon.te01*
00275 phycon.te003;
00276 }
00277 else if( te <= 1.4e4 )
00278 {
00279
00280 hmirat_v = 8.204e-17*phycon.sqrte/phycon.te10/phycon.te01*
00281 phycon.te003;
00282 }
00283 else
00284 {
00285
00286 hmirat_v = 5.44e-16*phycon.te20/phycon.te01;
00287 }
00288
00289 DEBUG_EXIT( "hmirat()" );
00290 return( hmirat_v );
00291 }
00292
00293
00294
00295 void hmole_init(void)
00296 {
00297
00298 DEBUG_ENTRY( "hmole_init()" );
00299
00300
00301 strcpy( hmi.chLab[ipMH], "H0 " );
00302 strcpy( hmi.chLab[ipMHp], "H+ " );
00303 strcpy( hmi.chLab[ipMHm], "H- " );
00304 strcpy( hmi.chLab[ipMH2g], "H2g " );
00305 strcpy( hmi.chLab[ipMH2p], "H2+ " );
00306 strcpy( hmi.chLab[ipMH3p], "H3+ " );
00307 strcpy( hmi.chLab[ipMH2s], "H2* " );
00308 strcpy( hmi.chLab[ipMHeHp], "HeH+" );
00309
00310
00311 hmi.rheph2hpheh = 0.;
00312
00313 DEBUG_EXIT( "hmole_init()" );
00314 return;
00315
00316 }
00317
00318
00319 void hmole_reactions( void )
00320 {
00321 static double teused=-1;
00322 double exph2,
00323 exph2s,
00324 exphp,
00325 ex3hp;
00326 long i;
00327 double h2esc,
00328 th2,
00329 cr_H2s ,
00330 cr_H2dis,
00331 cr_H2dis_ELWERT_H2g,
00332 cr_H2dis_ELWERT_H2s;
00333
00334 DEBUG_ENTRY( "hmole_reactions()" );
00335
00336
00337
00338
00339 if( phycon.te!=teused )
00340 {
00341 teused = phycon.te;
00342
00343
00344
00345
00346
00347 hmi.exphmi = sexp(8.745e3/phycon.te);
00348 if( hmi.exphmi > 0. )
00349 {
00350
00351 hmi.rel_pop_LTE_Hmin = SAHA/(phycon.te32*hmi.exphmi)*(1./(2.*2.));
00352 }
00353 else
00354 {
00355 hmi.rel_pop_LTE_Hmin = 0.;
00356 }
00357
00358
00359
00360
00361 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
00362 {
00363
00364 hmi.rel_pop_LTE_H2g = hmi.H2g_LTE_bigH2;
00365 hmi.rel_pop_LTE_H2s = hmi.H2s_LTE_bigH2;
00366 # if 0
00367 exph2 = sexp((5.195e4-hmi.H2_BigH2_H2g_av*T1CM)/phycon.te);
00368
00369 if( exph2 > 0. )
00370 hmi.rel_pop_LTE_H2g = SAHA/SDIV(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
00371 else
00372 hmi.rel_pop_LTE_H2g = 0.;
00373
00374
00375
00376
00377
00378 exph2s = sexp((5.195e4-hmi.H2_BigH2_H2s_av * T1CM)/phycon.te);
00379
00380 if( exph2s > 0. )
00381 hmi.rel_pop_LTE_H2s = SAHA/SDIV(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
00382 else
00383 hmi.rel_pop_LTE_H2s = 0.;
00384
00385
00386
00387
00388 # endif
00389 }
00390 else
00391 {
00392
00393
00394 exph2 = sexp((5.195e4)/phycon.te);
00395
00396
00397 if( exph2 > 0. )
00398 {
00399
00400 hmi.rel_pop_LTE_H2g = SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
00401 }
00402 else
00403 {
00404 hmi.rel_pop_LTE_H2g = 0.;
00405 }
00406
00407
00408
00409
00410
00411 exph2s = sexp(2.178e4/phycon.te);
00412
00413 if( exph2s > 0. )
00414 {
00415
00416 hmi.rel_pop_LTE_H2s = SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
00417 }
00418 else
00419 {
00420 hmi.rel_pop_LTE_H2s = 0.;
00421 }
00422 }
00423 {
00424
00425
00426
00427
00428 enum {DEBUG_LOC=false};
00429
00430 if( DEBUG_LOC && nzone>187&& iteration > 1)
00431 {
00432
00433 fprintf(ioQQQ,"ph2lteee\t%.2e\t%.1e\t%.1e\n",
00434 hmi.rel_pop_LTE_H2g,
00435 sexp(2.178e4/phycon.te),
00436 phycon.te);
00437 }
00438 }
00439
00440
00441
00442 exphp = sexp(3.072e4/phycon.te);
00443 if( exphp > 0. )
00444 {
00445
00446
00447 hmi.rel_pop_LTE_H2p = SAHA/(phycon.te32*exphp)*(4./(2.*1.))*3.634e-5;
00448 }
00449 else
00450 {
00451 hmi.rel_pop_LTE_H2p = 0.;
00452 }
00453
00454
00455
00456 ex3hp = sexp(1.882e4/phycon.te);
00457 if( ex3hp > 0. )
00458 {
00459
00460
00461 hmi.rel_pop_LTE_H3p = SAHA/(phycon.te32*ex3hp)*(4./(2.*1.))*3.634e-5;
00462 }
00463 else
00464 {
00465 hmi.rel_pop_LTE_H3p = 0.;
00466 }
00467 }
00468
00469
00470
00471 hmi.hminus_rad_attach = hmirat(phycon.te);
00472
00473 hmi.hmicol = hmi.hminus_rad_attach*EN1RYD*phycon.te*1.15e-5;
00474
00475
00476
00477
00478
00479 hmi.hminus_rad_attach *= dense.eden;
00480 hmi.hmicol *= dense.eden*hmi.Hmolec[ipMH];
00481
00482
00483
00484
00485
00486
00489
00490
00491
00492
00493
00494 hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] , opac.iphmop ,
00495 0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling );
00496
00497
00498 hmi.HMinus_photo_heat = thermal.HeatNet;
00499
00500 {
00501
00502
00503 enum {DEBUG_LOC=false};
00504
00505 if( DEBUG_LOC)
00506 {
00507 fprintf(ioQQQ,"hminphoto\t%li\t%li\t%.2e\n", nzone, conv.nPres2Ioniz , hmi.HMinus_photo_rate );
00508 }
00509 }
00510
00511
00512 hmi.HMinus_induc_rec_rate *= hmi.rel_pop_LTE_Hmin*dense.eden;
00513
00514
00516 hmi.HMinus_induc_rec_cooling *= hmi.rel_pop_LTE_Hmin*dense.eden*hmi.Hmolec[ipMH];
00517
00518 {
00519
00520
00521 enum {DEBUG_LOC=false};
00522
00523 if( DEBUG_LOC && nzone>400)
00524 {
00525 fprintf(ioQQQ,"hmoledebugg %.2f ",fnzone);
00526 GammaPrt(
00527 hmi.iphmin-1 , iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] , opac.iphmop ,
00528
00529 ioQQQ,
00530
00531 hmi.HMinus_photo_rate,
00532
00533 hmi.HMinus_photo_rate*0.05);
00534 }
00535 }
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547 hmi.UV_Cont_rel2_Habing_TH85_depth = 0.;
00548
00549
00550
00551
00552
00553 for( i=rfield.ipG0_TH85_lo; i < rfield.ipG0_TH85_hi; ++i )
00554 {
00555 hmi.UV_Cont_rel2_Habing_TH85_depth += (rfield.flux[i-1] + rfield.ConInterOut[i-1]+
00556 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
00557 }
00558
00559
00560
00561
00562
00563 hmi.UV_Cont_rel2_Habing_TH85_depth = (float)(hmi.UV_Cont_rel2_Habing_TH85_depth*EN1RYD/1.6e-3);
00564
00565 if( nzone<=1 )
00566 {
00567 hmi.UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_depth;
00568 }
00569
00570
00571
00572 hmi.UV_Cont_rel2_Habing_spec_depth = 0.;
00573 for( i=rfield.ipG0_spec_lo; i < rfield.ipG0_spec_hi; ++i )
00574 {
00575 hmi.UV_Cont_rel2_Habing_spec_depth += (rfield.flux[i-1] + rfield.ConInterOut[i-1]+
00576 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
00577 }
00578 hmi.UV_Cont_rel2_Habing_spec_depth = (float)(hmi.UV_Cont_rel2_Habing_spec_depth*EN1RYD/1.6e-3);
00579
00580
00581
00582 if( hmi.UV_Cont_rel2_Draine_DB96_face ==0 )
00583 {
00584
00585 for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; ++i )
00586 {
00587 hmi.UV_Cont_rel2_Draine_DB96_face += (rfield.flux[i-1] + rfield.ConInterOut[i-1]+
00588 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1]);
00589 }
00590
00591
00592
00593
00594
00595
00596
00597
00598 hmi.UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face/(1.232e7f * 1.71f);
00599 }
00600
00601
00602
00603
00604 hmi.H2Opacity = (float)(1.2e-14*(1e5/DoppVel.doppler[LIMELM+2]));
00605
00606 th2 = (colden.colden[ipCOL_H2g]+ colden.colden[ipCOL_H2s])*hmi.H2Opacity;
00607
00608
00609 h2esc = esc_PRD_1side(th2,1e-4);
00610
00611
00612
00613
00614
00615
00616
00617 hmi.H2_Solomon_dissoc_rate_TH85_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00618 hmi.H2_Solomon_dissoc_rate_TH85_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00619 hmi.H2_H2g_to_H2s_rate_TH85 = hmi.H2_Solomon_dissoc_rate_TH85_H2g*9.;
00620
00621
00622 hmi.H2_Solomon_dissoc_rate_BHT90_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00623 hmi.H2_Solomon_dissoc_rate_BHT90_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
00624 hmi.H2_H2g_to_H2s_rate_BHT90 = hmi.H2_Solomon_dissoc_rate_BHT90_H2g*9.;
00625
00626 {
00627
00628
00629
00630
00631 double x = (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) / 5e14;
00632 double sqrtx = sqrt(1.+x);
00633
00634 double b5 = DoppVel.doppler[LIMELM+2]/1e5;
00635
00636 double fshield = 0.965/POW2(1.+x/b5) + 0.035/sqrtx *
00637 sexp(8.5e-4*sqrtx);
00638
00639
00640
00641
00642
00643
00644 hmi.UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_face *
00645 (float)(sexp( opac.TauAbsFace[rfield.ip1000A-1] )/radius.r1r0sq);
00646
00647
00648
00649
00650 if( co.lgUMISTrates )
00651 {
00652
00653 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
00654 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
00655 }
00656 else
00657 {
00658
00659 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
00660 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
00661 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
00662 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
00663 }
00664
00665
00666
00667
00668
00669 hmi.H2_H2g_to_H2s_rate_BD96 = 5.67* hmi.H2_Solomon_dissoc_rate_BD96_H2g;
00670 }
00671
00672
00673 if( hmi.UV_Cont_rel2_Draine_DB96_face > SMALLFLOAT )
00674 {
00675
00676
00677
00678
00679
00680 double b5 = DoppVel.doppler[LIMELM+2]/1e5;
00681
00682
00683
00684
00685
00686
00687
00688 double x_H2g, x_H2s,
00689 fshield_H2g, fshield_H2s,
00690 f_H2s;
00691 static double a_H2g, a_H2s,
00692 e1_H2g, e1_H2s,
00693 e2_H2g,
00694 b_H2g,
00695 sl_H2g, sl_H2s,
00696 k_f_H2s,
00697 k_H2g_to_H2s,
00698 log_G0_face = -1;
00699
00700
00701
00702
00703
00704
00705
00706
00707 {
00708 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
00709 {
00710 log_G0_face = 0.;
00711 }
00712 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
00713 {
00714 log_G0_face = 7.;
00715 }
00716 else
00717 {
00718 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
00719 }
00720
00721
00722
00723
00724 a_H2g = 0.06 * log_G0_face + 1.32;
00725 a_H2s = 0.375 * log_G0_face + 2.125;
00726
00727 e1_H2g = -0.05 * log_G0_face + 2.25;
00728 e1_H2s = -0.125 * log_G0_face + 2.625;
00729
00730 e2_H2g = -0.005 * log_G0_face + 0.625;
00731
00732 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
00733
00734
00735 sl_H2g = 4.0e14;
00736 sl_H2s = 9.0e15;
00737
00738
00739 k_f_H2s = MAX2(0.1,2.375 * log_G0_face - 1.875 );
00740
00741
00742 k_H2g_to_H2s = MAX2(1.,-1.75 * log_G0_face + 11.25);
00743
00744
00745
00746
00747 }
00748
00749
00750 f_H2s = k_f_H2s * pow((double)hmi.UV_Cont_rel2_Draine_DB96_depth, 0.2 );
00751
00752
00753 x_H2g = (colden.colden[ipCOL_H2g]) / sl_H2g;
00754 x_H2s = (colden.colden[ipCOL_H2s]) / sl_H2s;
00755
00756
00757 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
00758 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
00759
00760
00761 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = a_H2g * 4.6e-11 * fshield_H2g * hmi.UV_Cont_rel2_Draine_DB96_depth;
00762 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = a_H2s * 4.6e-11 * fshield_H2s * (hmi.UV_Cont_rel2_Draine_DB96_depth + f_H2s);
00763
00764
00765 hmi.H2_H2g_to_H2s_rate_ELWERT = k_H2g_to_H2s * hmi.H2_Solomon_dissoc_rate_ELWERT_H2g;
00766
00767
00768
00769
00770 hmi.H2_photodissoc_ELWERT_H2s = hmi.UV_Cont_rel2_Draine_DB96_depth*1e-11;
00771 hmi.H2_photodissoc_ELWERT_H2g = hmi.H2_photodissoc_ELWERT_H2s * 1.0e-10;
00772 }
00773 else
00774 {
00775 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = 0.;
00776 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = 0.;
00777 hmi.H2_photodissoc_ELWERT_H2s = 0.;
00778 hmi.H2_photodissoc_ELWERT_H2g = 0.;
00779 }
00780
00781
00782 hmi.H2_photodissoc_TH85 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
00783
00784
00785 hmi.H2_photodissoc_BHT90 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
00786
00787
00788
00789
00790
00791
00792 cr_H2s = secondaries.x12tot*0.9 / 2. * hmi.lgLeidenCRHack;
00793
00794
00795 cr_H2dis = secondaries.x12tot*0.1 / 2. * hmi.lgLeidenCRHack;
00796
00797
00798
00799
00800
00801 cr_H2dis_ELWERT_H2g = secondaries.x12tot*5e-8 * hmi.lgLeidenCRHack;
00802 cr_H2dis_ELWERT_H2s = secondaries.x12tot*4e-2 * hmi.lgLeidenCRHack;
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
00815 {
00816
00817
00818
00819
00820
00821
00822 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_H2g;
00823 if( hmi.H2_Solomon_dissoc_rate_used_H2g <= 0. )
00824 {
00825
00826 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g;
00827 }
00828
00829 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_H2s;
00830 if( hmi.H2_Solomon_dissoc_rate_used_H2s <= 0. )
00831 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2g;
00832
00833
00834 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BigH2;
00835 if( hmi.H2_H2g_to_H2s_rate_used <= 0. )
00836 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85;
00837
00838
00839
00840
00841
00842 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BigH2_H2s;
00843 if( hmi.H2_photodissoc_used_H2s <= 0. )
00844 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
00845
00846
00847 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BigH2_H2g;
00848 if( hmi.H2_photodissoc_used_H2g <= 0. )
00849 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
00850 }
00851 else if( hmi.chH2_small_model_type == 'T' )
00852 {
00853
00854
00855 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g + cr_H2dis;
00856
00857 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2s + cr_H2dis;
00858 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85 + cr_H2s;
00859
00860
00861
00862 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
00863
00864
00865 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
00866 }
00867
00868 else if( hmi.chH2_small_model_type == 'H' )
00869 {
00870
00871
00872
00873 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BHT90_H2g + cr_H2dis;
00874
00875 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BHT90_H2s + cr_H2dis;
00876 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BHT90 + cr_H2s;
00877
00878
00879
00880 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BHT90;
00881
00882
00883 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BHT90*1.0e-10f;
00884 }
00885
00886 else if( hmi.chH2_small_model_type == 'B' )
00887 {
00888
00889
00890 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BD96_H2g + cr_H2dis;
00891
00892 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BD96_H2s + cr_H2dis;
00893
00894 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BD96 + cr_H2s;
00895
00896
00897
00898
00899 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
00900
00901
00902 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
00903 }
00904 else if( hmi.chH2_small_model_type == 'E' )
00905 {
00906
00907
00908 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + cr_H2dis_ELWERT_H2g;
00909 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_ELWERT_H2s + cr_H2dis_ELWERT_H2s;
00910 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_ELWERT + cr_H2s;
00911
00912
00913
00914
00915 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_ELWERT_H2s;
00916 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_ELWERT_H2g;
00917 }
00918 else
00919 TotalInsanity();
00920
00921 {
00922
00923 enum {DEBUG_LOC=false};
00924
00925 if( DEBUG_LOC && h2.lgH2ON )
00926 {
00927 fprintf(ioQQQ," Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
00928 hmi.H2_Solomon_dissoc_rate_TH85_H2g,
00929 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
00930 hmi.H2_Solomon_dissoc_rate_BigH2_H2g ,
00931 hmi.H2_H2g_to_H2s_rate_TH85 ,
00932 hmi.H2_H2g_to_H2s_rate_BigH2);
00933 }
00934 }
00935
00936 DEBUG_EXIT( "hmole_reactions()" );
00937 return;
00938 }
00939
00940