cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*hmole determine populations of hydrogen molecules */ 00004 /*hmole_old determine populations of hydrogen molecules */ 00005 /*hmole_init - initialize some hmole vars */ 00006 /*hmole_reactions update hmole reactions */ 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 /* Define to verify chemistry solution */ 00029 #if 0 00030 #if !defined(NDEBUG) 00031 /* >>chng 02 dec 21, line 14 changed from 00032 if(fabs(total) > 1e-20 && fabs(total) > 1e-14*mtotal) { 00033 to 00034 if(fabs(total) > 1e-20 && fabs(total) > 1e-8*mtotal) { */ 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 /* nothing */ 00058 #endif 00059 00060 #endif 00061 00062 void hmole( void ) 00063 { 00064 int nFixup, i; 00065 double error; 00066 realnum oxy_error=0.; 00067 static realnum abund0=-BIGFLOAT , abund1=-BIGFLOAT; 00068 realnum save1=dense.xIonDense[ipOXYGEN][1], 00069 save0=dense.xIonDense[ipOXYGEN][0]; 00070 00071 DEBUG_ENTRY( "hmole()" ); 00072 00073 /* will be used to keep track of neg solns */ 00074 nFixup = 0; 00075 error = 1.; 00076 00077 /* >>chng 04 apr 29, use PDR solution, scaled to density of H0, as first guess*/ 00078 /* if very first call on this sim, set H mol abundances to scaled TH85 PDR values */ 00079 if( conv.nTotalIoniz==0 && iteration==0 ) 00080 { 00081 realnum 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 /* we should have an H0 soln at this point */ 00091 ASSERT( dense.xIonDense[ipHYDROGEN][0]>SMALLFLOAT ); 00092 /* make sure order of molecules has not changed */ 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 /* update hmole reactions that depend on temperature */ 00110 hmole_reactions(); 00111 00112 /* will test against this error for quitting */ 00113 # define BIGERROR 1e-4 00114 /* >>chng 04 jul 20, upper limit had been 6, why? change to 20 00115 * to get closer to soln */ 00116 # define LIM_LOOP 20 00117 /* loop until run to limit, or no fix ups needed and error is small */ 00118 for(i=0; i < LIM_LOOP && ((error > BIGERROR||nFixup || oxy_error>conv.EdenErrorAllowed));i++) 00119 { 00120 { 00121 /* option to print deduced abundances */ 00122 /*@-redef@*/ 00123 enum {DEBUG_LOC=false}; 00124 /*@+redef@*/ 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 /* nFixup is number of times negative abundances were fixed, should be zero 00137 * at return for valid soln */ 00138 nFixup = 0; 00139 /* >>chng 04 jul 16, bring oxygen into the H solver */ 00140 if( !conv.lgSearch ) 00141 { 00142 IonOxyge(); 00143 } 00144 00145 save1 = dense.xIonDense[ipOXYGEN][1]; 00146 save0 = dense.xIonDense[ipOXYGEN][0]; 00147 /* >>chng 04 jul 29, necessary to damp out small changes in the O^0, O^+ density 00148 * when calling the hydrogen solver, since this can cause changes in H+/H0, which feed 00149 * back into the O+/O0 - small oscillations develop - shown in finely converged 00150 * dynamics models */ 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 /* conserve sum of O0 + O+ */ 00164 dense.xIonDense[ipOXYGEN][1] = abund1; 00165 /* error in this averaging */ 00166 oxy_error = (realnum)fabs(save1-abund1)/SDIV(dense.gas_phase[ipOXYGEN]); 00167 /*fprintf(ioQQQ,"DEBUG hmole\t%.2f\t%.5e\t%.5e\t%.5e\n", fnzone,save0 , save1 , oxy_error);*/ 00168 /*dense.xIonDense[ipOXYGEN][1] = abund1;*/ 00169 /* call the hydrogen solver */ 00170 hmole_step(&nFixup, &error); 00171 dense.xIonDense[ipOXYGEN][0] = save0; 00172 dense.xIonDense[ipOXYGEN][1] = save1; 00173 { 00174 /* option to print deduced abundances */ 00175 /*@-redef@*/ 00176 enum {DEBUG_LOC=false}; 00177 /*@+redef@*/ 00178 if( DEBUG_LOC /*&& (nzone>68)*/ ) 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 /*for( j=0; j<N_H_MOLEC; j++ ) 00198 fprintf(ioQQQ,"\t%.4e", hmi.Hmolec[j] );*/ 00199 fprintf(ioQQQ,"\n" ); 00200 } 00201 } 00202 } 00203 00204 if( (i == LIM_LOOP && error > BIGERROR) || nFixup != 0 ) 00205 { 00206 /* most of these failures occur just one time during search phase - 00207 * that is not a serious problem */ 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 /* check that density is still correct - CDEN is constant density */ 00223 if( strcmp( dense.chDenseLaw, "CDEN" )==0 ) 00224 /* check that we are conserving hydrogen nuclei */ 00225 ASSERT( fabs( dense.gas_phase[ipHYDROGEN] - dense.den0 )/ 00226 dense.gas_phase[ipHYDROGEN]<1e-4 ); 00227 00228 00229 /* total number of H per unit vol in molecules, 00230 * of course not including H0/H+ */ 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 /* remove the atom/ion which was just counted */ 00237 dense.xMolecules[ipHYDROGEN] -= (hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp]); 00238 00239 /* now add on all H in heavy element molecules */ 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 /*fprintf(ioQQQ," hmole return value is %.3e\n",timesc.time_H2_Dest_here);*/ 00246 return; 00247 } 00248 00249 /*hmirat compute radiative association rate for H- */ 00250 double hmirat(double te) 00251 { 00252 double hmirat_v; 00253 00254 DEBUG_ENTRY( "hmirat()" ); 00255 00256 /* fits to radiative association rate coefficients */ 00257 if( te < 31.62 ) 00258 { 00259 hmirat_v = 8.934e-18*phycon.sqrte*phycon.te003*phycon.te001* 00260 phycon.te001; 00261 } 00262 else if( te < 90. ) 00263 { 00264 hmirat_v = 5.159e-18*phycon.sqrte*phycon.te10*phycon.te03* 00265 phycon.te03*phycon.te003*phycon.te001; 00266 } 00267 else if( te < 1200. ) 00268 { 00269 hmirat_v = 2.042e-18*te/phycon.te10/phycon.te03; 00270 } 00271 else if( te < 3800. ) 00272 { 00273 hmirat_v = 8.861e-18*phycon.te70/phycon.te03/phycon.te01* 00274 phycon.te003; 00275 } 00276 else if( te <= 1.4e4 ) 00277 { 00278 /* following really only optimized up to 1e4 */ 00279 hmirat_v = 8.204e-17*phycon.sqrte/phycon.te10/phycon.te01* 00280 phycon.te003; 00281 } 00282 else 00283 { 00284 /* >>chng 00 sep 28, add this branch */ 00285 hmirat_v = 5.44e-16*phycon.te20/phycon.te01; 00286 } 00287 return( hmirat_v ); 00288 } 00289 00290 00291 /* hmole_init - one-time initialize of some hmole vars, called by cdinit */ 00292 void hmole_init(void) 00293 { 00294 00295 DEBUG_ENTRY( "hmole_init()" ); 00296 00297 /* do we need to fill in the molecular labels? */ 00298 strcpy( hmi.chLab[ipMH], "H0 " ); 00299 strcpy( hmi.chLab[ipMHp], "H+ " ); 00300 strcpy( hmi.chLab[ipMHm], "H- " ); 00301 strcpy( hmi.chLab[ipMH2g], "H2g " ); 00302 strcpy( hmi.chLab[ipMH2p], "H2+ " ); 00303 strcpy( hmi.chLab[ipMH3p], "H3+ " ); 00304 strcpy( hmi.chLab[ipMH2s], "H2* " ); 00305 strcpy( hmi.chLab[ipMHeHp], "HeH+" ); 00306 00307 /* rate for He+ + H2 -> */ 00308 hmi.rheph2hpheh = 0.; 00309 return; 00310 00311 } 00312 00313 /*hmole_reactions update hmole reactions */ 00314 void hmole_reactions( void ) 00315 { 00316 static double teused=-1; 00317 double exph2, 00318 exph2s, 00319 exphp, 00320 ex3hp; 00321 long i; 00322 double h2esc, 00323 th2, 00324 cr_H2s , 00325 cr_H2dis, 00326 cr_H2dis_ELWERT_H2g, 00327 cr_H2dis_ELWERT_H2s; 00328 00329 DEBUG_ENTRY( "hmole_reactions()" ); 00330 00331 /* everything here depends only on temperature - don't do anything if we don't 00332 * need to */ 00333 if( ! fp_equal( phycon.te, teused ) ) 00334 { 00335 teused = phycon.te; 00336 00337 /* get LTE populations */ 00338 /* related to population of H- in LTE 00339 * IP is 0.754 eV */ 00340 hmi.exphmi = sexp(8.745e3/phycon.te); 00341 if( hmi.exphmi > 0. ) 00342 { 00343 /* these are ratio n*(H-)/[ n*(ne) n*(Ho) ] */ 00344 hmi.rel_pop_LTE_Hmin = SAHA/(phycon.te32*hmi.exphmi)*(1./(2.*2.)); 00345 } 00346 else 00347 { 00348 hmi.rel_pop_LTE_Hmin = 0.; 00349 } 00350 00351 /* population of H2 in LTE 00352 * dissociation energy of H2g is 4.477eV, for TH85 model */ 00353 /* >>chng 05 oct 17, GS, averaged in big H2 in order to consider correct statistical weight*/ 00354 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 ) 00355 { 00356 /* the terms on the right were computed in the large molecule */ 00357 hmi.rel_pop_LTE_H2g = hmi.H2g_LTE_bigH2; 00358 hmi.rel_pop_LTE_H2s = hmi.H2s_LTE_bigH2; 00359 # if 0 00360 exph2 = sexp((5.195e4-hmi.H2_BigH2_H2g_av*T1CM)/phycon.te); 00361 /* >>chng 05 oct 21, GS, protect against zero Boltzmann factor */ 00362 if( exph2 > 0. ) 00363 hmi.rel_pop_LTE_H2g = SAHA/SDIV(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5; 00364 else 00365 hmi.rel_pop_LTE_H2g = 0.; 00366 /*fprintf(ioQQQ,"test\t%.2e\t%.1e\t%.1e\n", 00367 hmi.rel_pop_LTE_H2g, 00368 exph2, 00369 SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5);*/ 00370 00371 exph2s = sexp((5.195e4-hmi.H2_BigH2_H2s_av * T1CM)/phycon.te); 00372 /* >>chng 05 oct 21, GS, protect against zero Boltzmann factor */ 00373 if( exph2s > 0. ) 00374 hmi.rel_pop_LTE_H2s = SAHA/SDIV(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5; 00375 else 00376 hmi.rel_pop_LTE_H2s = 0.; 00377 /*fprintf(ioQQQ,"testh2s\t%.2e\t%.1e\t%.1e\n", 00378 hmi.rel_pop_LTE_H2s, 00379 exph2s, 00380 SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5);*/ 00381 # endif 00382 } 00383 else 00384 { 00385 00386 /* H2 ground */ 00387 exph2 = sexp((5.195e4)/phycon.te); 00388 /* >>chng 05 oct 17, GS, note that statical weight of H2g is assumed to be 1 if big H2 is not turned on*/ 00389 00390 if( exph2 > 0. ) 00391 { 00392 /* >>chng 05 oct 17, GS, note that statical weight of H2g is assumed to be 1 if big H2 is not turned on*/ 00393 hmi.rel_pop_LTE_H2g = SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5; 00394 } 00395 else 00396 { 00397 hmi.rel_pop_LTE_H2g = 0.; 00398 } 00399 00400 /* H2 star */ 00401 /* population of H2s in LTE 00402 * dissociation energy is 1.877eV, if h2s = 2.6eV, assumed for TH85 model */ 00403 /* >>chng 05 oct 17, GS, averaged in big H2 in order to consider correct statistical weight*/ 00404 exph2s = sexp(2.178e4/phycon.te); 00405 00406 if( exph2s > 0. ) 00407 { 00408 /* >>chng 05 oct 17, GS, note that statical weight of H2s is assumed to be 1 if big H2 is not turned on*/ 00409 hmi.rel_pop_LTE_H2s = SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5; 00410 } 00411 else 00412 { 00413 hmi.rel_pop_LTE_H2s = 0.; 00414 } 00415 } 00416 { 00417 /*@-redef@*/ 00418 /* often the H- route is the most efficient formation mechanism for H2, 00419 * will be through rate called Hmolec_old[ipMH]*hmi.assoc_detach 00420 * this debug print statement is to trace h2 oscillations */ 00421 enum {DEBUG_LOC=false}; 00422 /*@+redef@*/ 00423 if( DEBUG_LOC && nzone>187&& iteration > 1) 00424 { 00425 /* rapid increase in H2 density caused by rapid increase in hmi.rel_pop_LTE_H2g */ 00426 fprintf(ioQQQ,"ph2lteee\t%.2e\t%.1e\t%.1e\n", 00427 hmi.rel_pop_LTE_H2g, 00428 sexp(2.178e4/phycon.te), 00429 phycon.te); 00430 } 00431 } 00432 00433 /* population of H2+ in LTE, hmi.rel_pop_LTE_H2p is H_2+/H / H+ 00434 * dissociation energy is 2.647 */ 00435 exphp = sexp(3.072e4/phycon.te); 00436 if( exphp > 0. ) 00437 { 00438 /* stat weight of H2+ is 4 00439 * last factor was put in ver 85.23, missing before */ 00440 hmi.rel_pop_LTE_H2p = SAHA/(phycon.te32*exphp)*(4./(2.*1.))*3.634e-5; 00441 } 00442 else 00443 { 00444 hmi.rel_pop_LTE_H2p = 0.; 00445 } 00446 00447 /* related to population of H3+ in LTE, hmi.rel_pop_LTE_H3p is H_3+/( H2+ H+ ) 00448 * dissociation energy is 2.647 */ 00449 ex3hp = sexp(1.882e4/phycon.te); 00450 if( ex3hp > 0. ) 00451 { 00452 /* stat weight of H2+ is 4 00453 * last factor was put in ver 85.23, missing before */ 00454 hmi.rel_pop_LTE_H3p = SAHA/(phycon.te32*ex3hp)*(4./(2.*1.))*3.634e-5; 00455 } 00456 else 00457 { 00458 hmi.rel_pop_LTE_H3p = 0.; 00459 } 00460 } 00461 /* end constant temperature - */ 00462 00463 00464 hmi.hminus_rad_attach = hmirat(phycon.te); 00465 /* cooling due to radiative attachment */ 00466 hmi.hmicol = hmi.hminus_rad_attach*EN1RYD*phycon.te*1.15e-5; 00467 00468 /*fprintf(ioQQQ,"%.2e %.2e %.2e %.2e\n", phycon.te, hmi.hminus_rad_attach , hmi.hmicol, 00469 hmi.hmicol/(hmi.hminus_rad_attach*EN1RYD*phycon.te*1.15e-5) );*/ 00470 00471 /* get per unit vol */ 00472 hmi.hminus_rad_attach *= dense.eden; 00473 hmi.hmicol *= dense.eden*hmi.Hmolec[ipMH]; /* was dense.xIonDense[ipHYDROGEN][0]; */ 00474 00475 /* ================================================================= */ 00476 /* evaluate H- photodissociation rate, induced rec and rec cooling rates */ 00477 /* >>chng 00 dec 24, add test so that photo rate only reevaluated two times per zone. 00478 * in grain-free models this was sometimes dominated by Lya and so oscillated. 00479 * especially bad in primal.in - change 2 to 4 and primal.in will stop due to Lya oscil */ 00482 /* >>chng 02 feb 16, add damper on H- photo rate, wild oscillations in Lya photo rate in 00483 * grain free models */ 00484 00485 /*hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , nhe1Com.nhe1[0] , opac.iphmop , 00486 0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling );*/ 00487 hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] , opac.iphmop , 00488 0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling ); 00489 00490 /* save H- photodissociation heating */ 00491 hmi.HMinus_photo_heat = thermal.HeatNet; 00492 00493 { 00494 /* following should be set true to print populations */ 00495 /*@-redef@*/ 00496 enum {DEBUG_LOC=false}; 00497 /*@+redef@*/ 00498 if( DEBUG_LOC) 00499 { 00500 fprintf(ioQQQ,"hminphoto\t%li\t%li\t%.2e\n", nzone, conv.nPres2Ioniz , hmi.HMinus_photo_rate ); 00501 } 00502 } 00503 00504 /* induced recombination */ 00505 hmi.HMinus_induc_rec_rate *= hmi.rel_pop_LTE_Hmin*dense.eden; 00506 00507 /* induced recombination cooling per unit volume */ 00509 hmi.HMinus_induc_rec_cooling *= hmi.rel_pop_LTE_Hmin*dense.eden*hmi.Hmolec[ipMH]; /* dense.xIonDense[ipHYDROGEN][0]; */ 00510 00511 { 00512 /* following should be set true to debug H- photoionization rates */ 00513 /*@-redef@*/ 00514 enum {DEBUG_LOC=false}; 00515 /*@+redef@*/ 00516 if( DEBUG_LOC && nzone>400/*&& iteration > 1*/) 00517 { 00518 fprintf(ioQQQ,"hmoledebugg %.2f ",fnzone); 00519 GammaPrt( 00520 hmi.iphmin-1 , iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] , opac.iphmop , 00521 /* io unit we will write to */ 00522 ioQQQ, 00523 /* total photo rate from previous call to GammaK */ 00524 hmi.HMinus_photo_rate, 00525 /* we will print contributors that are more than this rate */ 00526 hmi.HMinus_photo_rate*0.05); 00527 } 00528 } 00529 /* add on high energy ionization, assume hydrogen cross section 00530 * n.b.; HGAMNC with secondaries */ 00531 /* >>chng 00 dec 24, above goes to HeI edge, no need for this, and was not important*/ 00532 /*hmi.HMinus_photo_rate += iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s];*/ 00533 00534 /* ================================================================= */ 00535 /* photodissociation by Lyman band absorption: esc prob treatment, 00536 * treatment based on 00537 * >>refer HI abs Tielens & Hollenbach 1985 ApJ 291, 722. */ 00538 /* do up to carbon photo edge if carbon is turned on */ 00539 /* >>>chng 00 apr 07, add test for whether element is turned on */ 00540 hmi.UV_Cont_rel2_Habing_TH85_depth = 0.; 00541 /* >>chng 00 apr 07 from explicit ipHeavy to ipLo */ 00542 /* find total intensity over carbon-ionizing continuum */ 00543 /* >>chng 03 jun 09, use exact bounds rather than CI photo threshold for lower bound */ 00544 /*for( i=ipLo; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )*/ 00545 /* the integral is from 6eV to 13.6, or 2060 - 912 Ang */ 00546 for( i=rfield.ipG0_TH85_lo; i < rfield.ipG0_TH85_hi; ++i ) 00547 { 00548 hmi.UV_Cont_rel2_Habing_TH85_depth += (rfield.flux[i-1] + rfield.ConInterOut[i-1]+ 00549 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1]; 00550 } 00551 00552 /* now convert to Habing ISM units 00553 * UV_Cont_rel2_Habing_TH85_face is FUV continuum relative to Habing value 00554 * 1.6e-3 ergs em-2 s-1 is the Habing 1968 field, quoted on page 723, end of first 00555 * full paragraph on left */ 00556 hmi.UV_Cont_rel2_Habing_TH85_depth = (realnum)(hmi.UV_Cont_rel2_Habing_TH85_depth*EN1RYD/1.6e-3); 00557 /* if start of calculation remember G0 at illuminated face */ 00558 if( nzone<=1 ) 00559 { 00560 hmi.UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_depth; 00561 } 00562 00563 00564 /* >>chng 05 jan 09, add special ver of G0 */ 00565 hmi.UV_Cont_rel2_Habing_spec_depth = 0.; 00566 for( i=rfield.ipG0_spec_lo; i < rfield.ipG0_spec_hi; ++i ) 00567 { 00568 hmi.UV_Cont_rel2_Habing_spec_depth += (rfield.flux[i-1] + rfield.ConInterOut[i-1]+ 00569 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1]; 00570 } 00571 hmi.UV_Cont_rel2_Habing_spec_depth = (realnum)(hmi.UV_Cont_rel2_Habing_spec_depth*EN1RYD/1.6e-3); 00572 00573 /* the Draine & Bertoldi version of the same thing, defined over their energy range */ 00574 /* >>chng 04 feb 07, only evaluate at the illuminated face */ 00575 if( hmi.UV_Cont_rel2_Draine_DB96_face ==0 ) 00576 { 00577 /* this is sum of photon number between 912A and 1110, as per BD96 */ 00578 for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; ++i ) 00579 { 00580 hmi.UV_Cont_rel2_Draine_DB96_face += (rfield.flux[i-1] + rfield.ConInterOut[i-1]+ 00581 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1]); 00582 } 00583 /* convert into scaled ISM background field, total number of photons over value for 1 ISM field, 00584 * the coefficient 1.232e7 is the number of photons over this wavelength range for 1H and is 00585 * given in BD96, page 225, 4th line from start of section 4, also page 272, table 1, 2nd line 00586 * from bottom */ 00587 /* >>chng 04 mar 16, introduce the 1.71 */ 00588 /* equation 20 of DB96 gives chi as flux over 1.21e7, to produce one Habing field. 00589 * to get the Draine field we need to further divide by 1.71 as stated on the first 00590 * line after equation 23 */ 00591 hmi.UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face/(1.232e7f * 1.71f); 00592 } 00593 00594 /* escape prob takes into account line shielding, 00595 * next is opacity then optical depth in H2 UV lines, using eqn A7 of TH85, 00596 * LIMELM+2 points to H2 */ 00597 hmi.H2Opacity = (realnum)(1.2e-14*(1e5/DoppVel.doppler[LIMELM+2])); 00598 /* the typical Lyman -Werner H2 line optical depth eq A7 of TH85a */ 00599 th2 = (colden.colden[ipCOL_H2g]+ colden.colden[ipCOL_H2s])*hmi.H2Opacity; 00600 /* the escape probability - chance that continuum photon will penetrate to 00601 * this depth to pump the Lyman Werner bands */ 00602 h2esc = esc_PRD_1side(th2,1e-4); 00603 00604 /* cross section is eqn A8 of 00605 * >>refer H2 dissoc Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722 00606 * branching ratio of 10% is already included, so 10x smaller than their number 00607 * 10% lead to dissociation through H_2 + h nu => 2H */ 00608 /* >>chng 05 mar 10, by TE, break into 2g and 2s 00609 * note use of same shielding column in below - can do far better */ 00610 hmi.H2_Solomon_dissoc_rate_TH85_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc; 00611 hmi.H2_Solomon_dissoc_rate_TH85_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc; 00612 hmi.H2_H2g_to_H2s_rate_TH85 = hmi.H2_Solomon_dissoc_rate_TH85_H2g*9.; 00613 00614 /* these are Burton et al. 1990 rates */ 00615 hmi.H2_Solomon_dissoc_rate_BHT90_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc; 00616 hmi.H2_Solomon_dissoc_rate_BHT90_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc; 00617 hmi.H2_H2g_to_H2s_rate_BHT90 = hmi.H2_Solomon_dissoc_rate_BHT90_H2g*9.; 00618 00619 { 00620 /* the following implements Drain & Bertoldi 1996, equation 37 from 00621 * >>refer H2 destruction Draine, B.T., & Bertoldi, F., 1996, ApJ, 468, 269-289 00622 * but the constant 4.6e-11 comes from Bertoldi & Draine equation 23, 00623 * this is the normalized H2 column density */ 00624 double x = (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) / 5e14; 00625 double sqrtx = sqrt(1.+x); 00626 /* Doppler with of H2 in km/s */ 00627 double b5 = DoppVel.doppler[LIMELM+2]/1e5; 00628 /* the molecular hydrogen line self-shielding factor */ 00629 double fshield = 0.965/POW2(1.+x/b5) + 0.035/sqrtx * 00630 sexp(8.5e-4*sqrtx); 00631 00632 /*double fshield = pow( MAX2(1.,colden.colden[ipCOLH2]/1e14) , -0.75 );*/ 00633 /* this is the Bertoldi & Draine version of the Habing field, 00634 * with dilution of radiation and extinction due to grains */ 00635 /* >>chng 04 apr 18, moved fshield, the line shielding factor, from this defn to 00636 * the following defn of dissociation rate, since following should measure continuum */ 00637 hmi.UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_face * 00638 (realnum)(sexp( opac.TauAbsFace[rfield.ip1000A-1] )/radius.r1r0sq); 00639 00640 /* the following comes from Bertoldi & Draine 1996, equation 23, 00641 * hmi.UV_Cont_rel2_Draine_DB96_depth already includes a factor of 1.71 to correct back from TH85 */ 00642 /* >>chng 05 mar 10, TE, break into 2s and 2s */ 00643 if( co.lgUMISTrates ) 00644 { 00645 /* this is default, when set Leiden hack UMIST rates not entered */ 00646 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield; 00647 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield; 00648 } 00649 else 00650 { 00651 /* done when set Leiden hack UMIST rates command entered */ 00652 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f) 00653 *sexp(3.02*rfield.extin_mag_V_point)* fshield; 00654 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f) 00655 *sexp(3.02*rfield.extin_mag_V_point)* fshield; 00656 } 00657 00658 /* BD do not give an excitation rate, so used 9 times the dissociation 00659 * rate by analogy with 90% going back into X, as per TH85 */ 00660 /*>>chng 04 feb 07, had used 90% relax into X from TH85, 00661 * BD say that 15% dissociate, so 85/15 = 5.67 is factor */ 00662 hmi.H2_H2g_to_H2s_rate_BD96 = 5.67* hmi.H2_Solomon_dissoc_rate_BD96_H2g; 00663 } 00664 00665 /* do Elwert approximation to the dissociation rate */ 00666 if( hmi.UV_Cont_rel2_Draine_DB96_face > SMALLFLOAT ) 00667 { 00668 /* this evaluates the new H2 dissociation rate by Torsten Elwert */ 00669 /* chng 05 jun 23, TE 00670 * >>chng 05 sep 13, update master source with now approximation */ 00671 00672 /* Doppler with of H2 in km/s */ 00673 double b5 = DoppVel.doppler[LIMELM+2]/1e5; 00674 00675 /* split the Solomon rates in H2g and H2s */ 00676 /* >>chng 05 oct 19, 00677 * >>chng 05 dec 05, TE, define new approximation for the heating due to the destruction of H2 00678 * use this approximation for the specified cloud parameters, otherwise 00679 * use limiting cases for 1 <= G0, G0 >= 1e7, n >= 1e7, n <= 1 */ 00680 00681 double x_H2g, x_H2s, 00682 fshield_H2g, fshield_H2s, 00683 f_H2s; 00684 static double a_H2g, a_H2s, 00685 e1_H2g, e1_H2s, 00686 e2_H2g, 00687 b_H2g, 00688 sl_H2g, sl_H2s, 00689 k_f_H2s, 00690 k_H2g_to_H2s, 00691 log_G0_face = -1; 00692 00693 /* define parameter range for the new approximation 00694 * test for G0 00695 *BUGFIX - this tested on lg_G0_face < 0 for initialization needed so did not work 00696 * in grids - change to evaluate in zone 0 */ 00697 /* >>chng 07 feb 24, BUGFIX crash when G0=0 at start and radiation 00698 * field builds up due to diffuse fields - soln is to always reevaluate */ 00699 /*if( !nzone )*/ 00700 { 00701 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.) 00702 { 00703 log_G0_face = 0.; 00704 } 00705 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7) 00706 { 00707 log_G0_face = 7.; 00708 } 00709 else 00710 { 00711 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face); 00712 } 00713 00714 /* terms only dependent on G0_face */ 00715 00716 /* coefficients and exponents */ 00717 a_H2g = 0.06 * log_G0_face + 1.32; 00718 a_H2s = 0.375 * log_G0_face + 2.125; 00719 00720 e1_H2g = -0.05 * log_G0_face + 2.25; 00721 e1_H2s = -0.125 * log_G0_face + 2.625; 00722 00723 e2_H2g = -0.005 * log_G0_face + 0.625; 00724 00725 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2; 00726 00727 /* scalelength for H2g and H2s */ 00728 sl_H2g = 4.0e14; 00729 sl_H2s = 9.0e15; 00730 00731 /* coefficient for 2nd term of Solomon H2s */ 00732 k_f_H2s = MAX2(0.1,2.375 * log_G0_face - 1.875 ); 00733 00734 /* coefficient for branching ratio */ 00735 k_H2g_to_H2s = MAX2(1.,-1.75 * log_G0_face + 11.25); 00736 00737 /*fprintf( ioQQQ, "e1_H2g%.2e, e1_H2s%.2e, e2_H2g%.2e, b_H2g%.2e, a_H2g%.2e, a_H2s%.2e,sl_H2g: %.2e,sl_H2s: %.2e\n", 00738 e1_H2g, e1_H2s, e2_H2g, b_H2g, a_H2g, a_H2s, sl_H2g, sl_H2s); 00739 */ 00740 } 00741 00742 /* Solomon H2s ~G0^0.2 at large depth*/ 00743 f_H2s = k_f_H2s * pow((double)hmi.UV_Cont_rel2_Draine_DB96_depth, 0.2 ); 00744 00745 /* scale length for absorption of UV lines */ 00746 x_H2g = (colden.colden[ipCOL_H2g]) / sl_H2g; 00747 x_H2s = (colden.colden[ipCOL_H2s]) / sl_H2s; 00748 00749 /* the molecular hydrogen line self-shielding factor */ 00750 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g); 00751 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s); 00752 00753 /* the Elwert Solomon rates for H2g and H2s hmi.chH2_small_model_type == 'E' */ 00754 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = a_H2g * 4.6e-11 * fshield_H2g * hmi.UV_Cont_rel2_Draine_DB96_depth; 00755 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = a_H2s * 4.6e-11 * fshield_H2s * (hmi.UV_Cont_rel2_Draine_DB96_depth + f_H2s); 00756 00757 /* assume branching ratio dependent on G0*/ 00758 hmi.H2_H2g_to_H2s_rate_ELWERT = k_H2g_to_H2s * hmi.H2_Solomon_dissoc_rate_ELWERT_H2g; 00759 00760 /* use G0_BD96 as this definition declines faster with depth which is physical as 00761 * the longer wavelengths in the definition of G0_TH85 cannot dissociate 00762 * H2s directly */ 00763 hmi.H2_photodissoc_ELWERT_H2s = hmi.UV_Cont_rel2_Draine_DB96_depth*1e-11; 00764 hmi.H2_photodissoc_ELWERT_H2g = hmi.H2_photodissoc_ELWERT_H2s * 1.0e-10; 00765 } 00766 else 00767 { 00768 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = 0.; 00769 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = 0.; 00770 hmi.H2_photodissoc_ELWERT_H2s = 0.; 00771 hmi.H2_photodissoc_ELWERT_H2g = 0.; 00772 } 00773 00774 /* this is rate of photodissociation of H2*, A12 of TH85 */ 00775 hmi.H2_photodissoc_TH85 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11; 00776 00777 /* dissociation rate from Burton et al. 1990 */ 00778 hmi.H2_photodissoc_BHT90 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11; 00779 00780 /* rates for cosmic ray excitation of singlet bound electronic bound excited states 00781 * only add this to small molecule since automatically included in large 00782 *>>refer H2 cr excit Dalgarno, A., Yan, Min, & Liu, Weihong 1999, ApJS, 125, 237 00783 * this is excitation of H2* */ 00784 /* >>chng 05 sep 13, do not include this process when Leiden hacks are in place */ 00785 cr_H2s = secondaries.x12tot*0.9 / 2. * hmi.lgLeidenCRHack; 00786 /* this is the fraction that dissociate */ 00787 /* >>chng 05 sep 13, do not include this process when Leiden hacks are in place */ 00788 cr_H2dis = secondaries.x12tot*0.1 / 2. * hmi.lgLeidenCRHack; 00789 00790 /* >>chng 05 sep 13, TE, synchronize treatment of CR */ 00791 /* cosmic ray rates for dissociation of ground and H2s 00792 * two factors done to agree with large H2 deep in the cloud where 00793 * cosmic rays are important */ 00794 cr_H2dis_ELWERT_H2g = secondaries.x12tot*5e-8 * hmi.lgLeidenCRHack; 00795 cr_H2dis_ELWERT_H2s = secondaries.x12tot*4e-2 * hmi.lgLeidenCRHack; 00796 00797 /* at this point there are two or three independent estimates of the H2 dissociation rate. 00798 * if the large H2 molecule is on, then H2 Solomon rates has been defined in the last 00799 * call to the large molecule. Just above we have defined hmi.H2_Solomon_dissoc_rate_TH85, 00800 * the dissociation rate from Tielens & Hollenbach 1985, and hmi.H2_Solomon_dissoc_rate_BD96, 00801 * the rate from Bertoldi & Draine 1996. We can use any defined rate. If the big H2 00802 * molecule is on, use its rate. If not, for now use the TH85 rate, since that is the 00803 * rate we always have used in the past. 00804 * The actual rate we will use is given by hmi.H2_Solomon_dissoc_rate_used 00805 */ 00806 /* this is the Solomon process dissociation rate actually used */ 00807 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 ) 00808 { 00809 /* only update after big H2 molecule has been evaluated, 00810 * when very little H2 and big molecule not used, leave at previous (probably TH85) value, 00811 * since that value is always known */ 00812 00813 /* Solomon process rate from X into the X continuum with units s-1 00814 * rates are total rate, and rates from H2g and H2s */ 00815 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_H2g; 00816 if( hmi.H2_Solomon_dissoc_rate_used_H2g <= 0. ) 00817 { 00818 /*fprintf(ioQQQ,"DEBUG his sol dis <0\n");*/ 00819 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g; 00820 } 00821 00822 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_H2s; 00823 if( hmi.H2_Solomon_dissoc_rate_used_H2s <= 0. ) 00824 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2g; 00825 00826 /* photoexcitation from H2g to H2s */ 00827 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BigH2; 00828 if( hmi.H2_H2g_to_H2s_rate_used <= 0. ) 00829 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85; 00830 00831 /* add up H2s + hnu (continuum) => 2H + KE, continuum photodissociation, 00832 * this is not the Solomon process, true continuum, units s-1 */ 00833 /* only include rates from H2s since this is only open channel, this process is well 00834 * shielded against Lyman continuum destruction by atomic hydrogen */ 00835 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BigH2_H2s; 00836 if( hmi.H2_photodissoc_used_H2s <= 0. ) 00837 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85; 00838 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts 00839 * for unfavorable wavelength range of G0*/ 00840 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BigH2_H2g; 00841 if( hmi.H2_photodissoc_used_H2g <= 0. ) 00842 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f; 00843 } 00844 else if( hmi.chH2_small_model_type == 'T' ) 00845 { 00846 /* the TH85 rate */ 00847 /*>>chng 05 jun 23, add cosmic rays */ 00848 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g + cr_H2dis; 00849 /* >>chng 05 sep 13, cr_H2dis was not included */ 00850 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2s + cr_H2dis; 00851 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85 + cr_H2s; 00852 00853 /* continuum photodissociation H2s + hnu => 2H, , 00854 * this is not the Solomon process, true continuum, units s-1 */ 00855 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85; 00856 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts 00857 * for unfavorable wavelength range of G0*/ 00858 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f; 00859 } 00860 00861 else if( hmi.chH2_small_model_type == 'H' ) 00862 { 00863 /* the improved H2 formalism given by 00864 *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M 00865 >>refcon 1990, ApJ, 365, 620 */ 00866 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BHT90_H2g + cr_H2dis; 00867 /* >>chng 05 sep 13, cr_H2dis was not included */ 00868 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BHT90_H2s + cr_H2dis; 00869 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BHT90 + cr_H2s; 00870 00871 /* continuum photodissociation H2s + hnu => 2H, , 00872 * this is not the Solomon process, true continuum, units s-1 */ 00873 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BHT90; 00874 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts 00875 * for unfavorable wavelength range of G0*/ 00876 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BHT90*1.0e-10f; 00877 } 00878 00879 else if( hmi.chH2_small_model_type == 'B' ) 00880 { 00881 /* the Bertoldi & Draine rate - this is the default */ 00882 /*>>chng 05 jun 23, add cosmic rays */ 00883 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BD96_H2g + cr_H2dis; 00884 /* >>chng 05 sep 13, cr_H2dis was not included */ 00885 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BD96_H2s + cr_H2dis; 00886 /* they did not do the excitation or dissoc rate, so use TH85 */ 00887 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BD96 + cr_H2s; 00888 00889 00890 /* continuum photodissociation H2s + hnu => 2H, , 00891 * this is not the Solomon process, true continuum, units s-1 */ 00892 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85; 00893 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts 00894 * for unfavorable wavelength range of G0*/ 00895 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f; 00896 } 00897 else if( hmi.chH2_small_model_type == 'E' ) 00898 { 00899 /* the Elwert et al. rate 00900 *>>chng 05 jun 23, add cosmic rays */ 00901 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + cr_H2dis_ELWERT_H2g; 00902 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_ELWERT_H2s + cr_H2dis_ELWERT_H2s; 00903 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_ELWERT + cr_H2s; 00904 00905 00906 /* continuum photodissociation H2s + hnu => 2H, , 00907 * this is not the Solomon process, true continuum, units s-1 */ 00908 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_ELWERT_H2s; 00909 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_ELWERT_H2g; 00910 } 00911 else 00912 TotalInsanity(); 00913 00914 { 00915 /*@-redef@*/ 00916 enum {DEBUG_LOC=false}; 00917 /*@+redef@*/ 00918 if( DEBUG_LOC && h2.lgH2ON ) 00919 { 00920 fprintf(ioQQQ," Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n", 00921 hmi.H2_Solomon_dissoc_rate_TH85_H2g, 00922 hmi.H2_Solomon_dissoc_rate_BD96_H2g, 00923 hmi.H2_Solomon_dissoc_rate_BigH2_H2g , 00924 hmi.H2_H2g_to_H2s_rate_TH85 , 00925 hmi.H2_H2g_to_H2s_rate_BigH2); 00926 } 00927 } 00928 return; 00929 } 00930