cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_h_drive.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*hmole determine populations of hydrogen molecules */
4 /*hmole_old determine populations of hydrogen molecules */
5 /*hmole_init - initialize some hmole vars */
6 /*hmole_reactions update hmole reactions */
7 #include "cddefines.h"
8 #include "physconst.h"
9 #include "dense.h"
10 #include "called.h"
11 #include "gammas.h"
12 #include "colden.h"
13 #include "thermal.h"
14 #include "secondaries.h"
15 #include "h2.h"
16 #include "mole.h"
17 #include "radius.h"
18 #include "doppvel.h"
19 #include "rfield.h"
20 #include "ionbal.h"
21 #include "rt.h"
22 #include "opacity.h"
23 #include "iso.h"
24 #include "conv.h"
25 #include "phycon.h"
26 #include "hmi.h"
27 
28 /* Define to verify chemistry solution */
29 #if 0
30 #if !defined(NDEBUG)
31 /* >>chng 02 dec 21, line 14 changed from
32  if(fabs(total) > 1e-20 && fabs(total) > 1e-14*mtotal) {
33  to
34  if(fabs(total) > 1e-20 && fabs(total) > 1e-8*mtotal) { */
35 #define AUDIT(a) { \
36  double total, mtotal; \
37  for( i=0;i<N_H_MOLEC;i++) { \
38  total = 0.; \
39  for( j=0;j<N_H_MOLEC;j++) { \
40  total += c[i][j]*nprot[j]; \
41  } \
42  if( fabs(total) > 1e-6*fabs(c[i][i]*nprot[i])) { \
43  fprintf(ioQQQ,"PROBLEM Subtotal1 %c %.2e\n",a,fabs(total)/fabs(c[i][i]*nprot[i])); \
44  fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*nprot[i]); \
45  } \
46  } \
47  total = mtotal = 0.;for( j=0;j<N_H_MOLEC;j++) { total += bvec[j]*nprot[j]; mtotal += fabs(bvec[j]*nprot[j]); }\
48  if( fabs(total) > 1e-30 && fabs(total) > 1e-10*rtot) { \
49  fprintf(ioQQQ,"PROBLEM Subtotal2 %c %.2e\n",a,fabs(total)/mtotal); \
50  fprintf(ioQQQ,"RHS Total %g cf %g\n",total,mtotal); \
51  } else if( a == '.' && fabs(total) > 1e-7*mtotal) { \
52  fprintf(ioQQQ,"PROBLEM zone %li Hmole RHS conservation error %.2e of %.2e\n",nzone,total,mtotal); \
53  fprintf(ioQQQ,"(may be due to high rate equilibrium reactions)\n"); \
54  } \
55  }
56 #else
57 #define AUDIT /* nothing */
58 #endif
59 
60 #endif
61 
62 void hmole( void )
63 {
64  int nFixup, i;
65  double error;
66  realnum oxy_error=0.;
67  static realnum abund0=-BIGFLOAT , abund1=-BIGFLOAT;
68  realnum save1=dense.xIonDense[ipOXYGEN][1],
69  save0=dense.xIonDense[ipOXYGEN][0];
70 
71  DEBUG_ENTRY( "hmole()" );
72 
73  /* will be used to keep track of neg solns */
74  nFixup = 0;
75  error = 1.;
76 
77  /* >>chng 04 apr 29, use PDR solution, scaled to density of H0, as first guess*/
78  /* if very first call on this sim, set H mol abundances to scaled TH85 PDR values */
79  if( conv.nTotalIoniz==0 && iteration==0 )
80  {
81  realnum pdr_mole_h2[N_H_MOLEC] = {1.00E+00f,
82  3.18E-05f,
83  1.95E-11f,
84  4.00E-08f,
85  1.08E-14f,
86  1.08E-20f,
87  3.85E-07f,
88  8.04E-14f};
89 
90  /* we should have an H0 soln at this point */
92  /* make sure order of molecules has not changed */
93  ASSERT( ipMH==0&&
94  ipMHp == 1&&
95  ipMHm == 2&&
96  ipMH2g == 3&&
97  ipMH2p == 4&&
98  ipMH3p == 5&&
99  ipMH2s == 6&&
100  ipMHeHp== 7&&
101  N_H_MOLEC==8 );
102 
103  for( i=0; i<N_H_MOLEC; ++i )
104  {
105  hmi.Hmolec[i] = dense.xIonDense[ipHYDROGEN][0]*pdr_mole_h2[i];
106  }
107  }
108 
109  /* update hmole reactions that depend on temperature */
110  hmole_reactions();
111 
112  /* will test against this error for quitting */
113 # define BIGERROR 1e-4
114  /* >>chng 04 jul 20, upper limit had been 6, why? change to 20
115  * to get closer to soln */
116 # define LIM_LOOP 20
117  /* loop until run to limit, or no fix ups needed and error is small */
118  for(i=0; i < LIM_LOOP && ((error > BIGERROR||nFixup || oxy_error>conv.EdenErrorAllowed));i++)
119  {
120  {
121  /* option to print deduced abundances */
122  /*@-redef@*/
123  enum {DEBUG_LOC=false};
124  /*@+redef@*/
125  if( DEBUG_LOC && (nzone>140) )
126  {
127  fprintf(ioQQQ,"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
128  fnzone,
129  phycon.te,
130  dense.eden);
131  for( i=0; i<N_H_MOLEC; i++ )
132  fprintf(ioQQQ,"\t%.2e", hmi.Hmolec[i] );
133  fprintf(ioQQQ,"\n" );
134  }
135  }
136  /* nFixup is number of times negative abundances were fixed, should be zero
137  * at return for valid soln */
138  nFixup = 0;
139  /* >>chng 04 jul 16, bring oxygen into the H solver */
140  if( !conv.lgSearch )
141  {
142  IonOxyge();
143  }
144 
145  save1 = dense.xIonDense[ipOXYGEN][1];
146  save0 = dense.xIonDense[ipOXYGEN][0];
147  /* >>chng 04 jul 29, necessary to damp out small changes in the O^0, O^+ density
148  * when calling the hydrogen solver, since this can cause changes in H+/H0, which feed
149  * back into the O+/O0 - small oscillations develop - shown in finely converged
150  * dynamics models */
151  if( nzone )
152  {
153 # define OLD 0.75f
154  abund0 = abund0*OLD+dense.xIonDense[ipOXYGEN][0]*(1.f-OLD);
155  abund1 = abund1*OLD+dense.xIonDense[ipOXYGEN][1]*(1.f-OLD);
156  }
157  else
158  {
159  abund0 = dense.xIonDense[ipOXYGEN][0];
160  abund1 = dense.xIonDense[ipOXYGEN][1];
161  }
162  dense.xIonDense[ipOXYGEN][0] = abund0;
163  /* conserve sum of O0 + O+ */
164  dense.xIonDense[ipOXYGEN][1] = abund1;
165  /* error in this averaging */
166  oxy_error = (realnum)fabs(save1-abund1)/SDIV(dense.gas_phase[ipOXYGEN]);
167  /*fprintf(ioQQQ,"DEBUG hmole\t%.2f\t%.5e\t%.5e\t%.5e\n", fnzone,save0 , save1 , oxy_error);*/
168  /*dense.xIonDense[ipOXYGEN][1] = abund1;*/
169  /* call the hydrogen solver */
170  hmole_step(&nFixup, &error);
171  dense.xIonDense[ipOXYGEN][0] = save0;
172  dense.xIonDense[ipOXYGEN][1] = save1;
173  {
174  /* option to print deduced abundances */
175  /*@-redef@*/
176  enum {DEBUG_LOC=false};
177  /*@+redef@*/
178  if( DEBUG_LOC /*&& (nzone>68)*/ )
179  {
180  fprintf(ioQQQ,"DEBUG hmole out\t%i\t%.2f\t%.5e\t%.5e",
181  i,
182  fnzone,
183  phycon.te,
184  dense.eden);
185  fprintf(ioQQQ,
186  "\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",
187  error,
190  mole.sink[ipHYDROGEN][0],
191  mole.sink[ipHYDROGEN][1],
192  mole.source[ipHYDROGEN][0] ,
193  mole.source[ipHYDROGEN][1] ,
196 
197  /*for( j=0; j<N_H_MOLEC; j++ )
198  fprintf(ioQQQ,"\t%.4e", hmi.Hmolec[j] );*/
199  fprintf(ioQQQ,"\n" );
200  }
201  }
202  }
203 
204  if( (i == LIM_LOOP && error > BIGERROR) || nFixup != 0 )
205  {
206  /* most of these failures occur just one time during search phase -
207  * that is not a serious problem */
208  conv.lgConvPops = false;
209 
210  if( !conv.lgSearch && called.lgTalk )
211  {
212  fprintf(ioQQQ," PROBLEM hmole, zone %li: %d iters, %d bad; final error: %g lgSearch %i\n",
213  nzone,
214  i,
215  nFixup,
216  error,
217  conv.lgSearch);
218  }
219  ConvFail( "pops" , "Hmole");
220  }
221 
222  /* check that density is still correct - CDEN is constant density */
223  if( strcmp( dense.chDenseLaw, "CDEN" )==0 )
224  /* check that we are conserving hydrogen nuclei */
226  dense.gas_phase[ipHYDROGEN]<1e-4 );
227 
228 
229  /* total number of H per unit vol in molecules,
230  * of course not including H0/H+ */
232  for(i=0;i<N_H_MOLEC;i++)
233  {
235  }
236  /* remove the atom/ion which was just counted */
238 
239  /* now add on all H in heavy element molecules */
240  for( i=0; i < mole.num_comole_calc; i++ )
241  {
243  }
244 
245  /*fprintf(ioQQQ," hmole return value is %.3e\n",timesc.time_H2_Dest_here);*/
246  return;
247 }
248 
249 /*hmirat compute radiative association rate for H- */
250 double hmirat(double te)
251 {
252  double hmirat_v;
253 
254  DEBUG_ENTRY( "hmirat()" );
255 
256  /* fits to radiative association rate coefficients */
257  if( te < 31.62 )
258  {
259  hmirat_v = 8.934e-18*phycon.sqrte*phycon.te003*phycon.te001*
260  phycon.te001;
261  }
262  else if( te < 90. )
263  {
264  hmirat_v = 5.159e-18*phycon.sqrte*phycon.te10*phycon.te03*
266  }
267  else if( te < 1200. )
268  {
269  hmirat_v = 2.042e-18*te/phycon.te10/phycon.te03;
270  }
271  else if( te < 3800. )
272  {
273  hmirat_v = 8.861e-18*phycon.te70/phycon.te03/phycon.te01*
274  phycon.te003;
275  }
276  else if( te <= 1.4e4 )
277  {
278  /* following really only optimized up to 1e4 */
279  hmirat_v = 8.204e-17*phycon.sqrte/phycon.te10/phycon.te01*
280  phycon.te003;
281  }
282  else
283  {
284  /* >>chng 00 sep 28, add this branch */
285  hmirat_v = 5.44e-16*phycon.te20/phycon.te01;
286  }
287  return( hmirat_v );
288 }
289 
290 
291 /* hmole_init - one-time initialize of some hmole vars, called by cdinit */
292 void hmole_init(void)
293 {
294 
295  DEBUG_ENTRY( "hmole_init()" );
296 
297  /* do we need to fill in the molecular labels? */
298  strcpy( hmi.chLab[ipMH], "H0 " );
299  strcpy( hmi.chLab[ipMHp], "H+ " );
300  strcpy( hmi.chLab[ipMHm], "H- " );
301  strcpy( hmi.chLab[ipMH2g], "H2g " );
302  strcpy( hmi.chLab[ipMH2p], "H2+ " );
303  strcpy( hmi.chLab[ipMH3p], "H3+ " );
304  strcpy( hmi.chLab[ipMH2s], "H2* " );
305  strcpy( hmi.chLab[ipMHeHp], "HeH+" );
306 
307  /* rate for He+ + H2 -> */
308  hmi.rheph2hpheh = 0.;
309  return;
310 
311 }
312 
313 /*hmole_reactions update hmole reactions */
314 void hmole_reactions( void )
315 {
316  static double teused=-1;
317  double exph2,
318  exph2s,
319  exphp,
320  ex3hp;
321  long i;
322  double h2esc,
323  th2,
324  cr_H2s ,
325  cr_H2dis,
326  cr_H2dis_ELWERT_H2g,
327  cr_H2dis_ELWERT_H2s;
328 
329  DEBUG_ENTRY( "hmole_reactions()" );
330 
331  /* everything here depends only on temperature - don't do anything if we don't
332  * need to */
333  if( ! fp_equal( phycon.te, teused ) )
334  {
335  teused = phycon.te;
336 
337  /* get LTE populations */
338  /* related to population of H- in LTE
339  * IP is 0.754 eV */
340  hmi.exphmi = sexp(8.745e3/phycon.te);
341  if( hmi.exphmi > 0. )
342  {
343  /* these are ratio n*(H-)/[ n*(ne) n*(Ho) ] */
344  hmi.rel_pop_LTE_Hmin = SAHA/(phycon.te32*hmi.exphmi)*(1./(2.*2.));
345  }
346  else
347  {
348  hmi.rel_pop_LTE_Hmin = 0.;
349  }
350 
351  /* population of H2 in LTE
352  * dissociation energy of H2g is 4.477eV, for TH85 model */
353  /* >>chng 05 oct 17, GS, averaged in big H2 in order to consider correct statistical weight*/
355  {
356  /* the terms on the right were computed in the large molecule */
359 # if 0
360  exph2 = sexp((5.195e4-hmi.H2_BigH2_H2g_av*T1CM)/phycon.te);
361  /* >>chng 05 oct 21, GS, protect against zero Boltzmann factor */
362  if( exph2 > 0. )
363  hmi.rel_pop_LTE_H2g = SAHA/SDIV(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
364  else
365  hmi.rel_pop_LTE_H2g = 0.;
366  /*fprintf(ioQQQ,"test\t%.2e\t%.1e\t%.1e\n",
367  hmi.rel_pop_LTE_H2g,
368  exph2,
369  SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5);*/
370 
371  exph2s = sexp((5.195e4-hmi.H2_BigH2_H2s_av * T1CM)/phycon.te);
372  /* >>chng 05 oct 21, GS, protect against zero Boltzmann factor */
373  if( exph2s > 0. )
374  hmi.rel_pop_LTE_H2s = SAHA/SDIV(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
375  else
376  hmi.rel_pop_LTE_H2s = 0.;
377  /*fprintf(ioQQQ,"testh2s\t%.2e\t%.1e\t%.1e\n",
378  hmi.rel_pop_LTE_H2s,
379  exph2s,
380  SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5);*/
381 # endif
382  }
383  else
384  {
385 
386  /* H2 ground */
387  exph2 = sexp((5.195e4)/phycon.te);
388  /* >>chng 05 oct 17, GS, note that statical weight of H2g is assumed to be 1 if big H2 is not turned on*/
389 
390  if( exph2 > 0. )
391  {
392  /* >>chng 05 oct 17, GS, note that statical weight of H2g is assumed to be 1 if big H2 is not turned on*/
393  hmi.rel_pop_LTE_H2g = SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
394  }
395  else
396  {
397  hmi.rel_pop_LTE_H2g = 0.;
398  }
399 
400  /* H2 star */
401  /* population of H2s in LTE
402  * dissociation energy is 1.877eV, if h2s = 2.6eV, assumed for TH85 model */
403  /* >>chng 05 oct 17, GS, averaged in big H2 in order to consider correct statistical weight*/
404  exph2s = sexp(2.178e4/phycon.te);
405 
406  if( exph2s > 0. )
407  {
408  /* >>chng 05 oct 17, GS, note that statical weight of H2s is assumed to be 1 if big H2 is not turned on*/
409  hmi.rel_pop_LTE_H2s = SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
410  }
411  else
412  {
413  hmi.rel_pop_LTE_H2s = 0.;
414  }
415  }
416  {
417  /*@-redef@*/
418  /* often the H- route is the most efficient formation mechanism for H2,
419  * will be through rate called Hmolec_old[ipMH]*hmi.assoc_detach
420  * this debug print statement is to trace h2 oscillations */
421  enum {DEBUG_LOC=false};
422  /*@+redef@*/
423  if( DEBUG_LOC && nzone>187&& iteration > 1)
424  {
425  /* rapid increase in H2 density caused by rapid increase in hmi.rel_pop_LTE_H2g */
426  fprintf(ioQQQ,"ph2lteee\t%.2e\t%.1e\t%.1e\n",
428  sexp(2.178e4/phycon.te),
429  phycon.te);
430  }
431  }
432 
433  /* population of H2+ in LTE, hmi.rel_pop_LTE_H2p is H_2+/H / H+
434  * dissociation energy is 2.647 */
435  exphp = sexp(3.072e4/phycon.te);
436  if( exphp > 0. )
437  {
438  /* stat weight of H2+ is 4
439  * last factor was put in ver 85.23, missing before */
440  hmi.rel_pop_LTE_H2p = SAHA/(phycon.te32*exphp)*(4./(2.*1.))*3.634e-5;
441  }
442  else
443  {
444  hmi.rel_pop_LTE_H2p = 0.;
445  }
446 
447  /* related to population of H3+ in LTE, hmi.rel_pop_LTE_H3p is H_3+/( H2+ H+ )
448  * dissociation energy is 2.647 */
449  ex3hp = sexp(1.882e4/phycon.te);
450  if( ex3hp > 0. )
451  {
452  /* stat weight of H2+ is 4
453  * last factor was put in ver 85.23, missing before */
454  hmi.rel_pop_LTE_H3p = SAHA/(phycon.te32*ex3hp)*(4./(2.*1.))*3.634e-5;
455  }
456  else
457  {
458  hmi.rel_pop_LTE_H3p = 0.;
459  }
460  }
461  /* end constant temperature - */
462 
463 
465  /* cooling due to radiative attachment */
467 
468  /*fprintf(ioQQQ,"%.2e %.2e %.2e %.2e\n", phycon.te, hmi.hminus_rad_attach , hmi.hmicol,
469  hmi.hmicol/(hmi.hminus_rad_attach*EN1RYD*phycon.te*1.15e-5) );*/
470 
471  /* get per unit vol */
473  hmi.hmicol *= dense.eden*hmi.Hmolec[ipMH]; /* was dense.xIonDense[ipHYDROGEN][0]; */
474 
475  /* ================================================================= */
476  /* evaluate H- photodissociation rate, induced rec and rec cooling rates */
477  /* >>chng 00 dec 24, add test so that photo rate only reevaluated two times per zone.
478  * in grain-free models this was sometimes dominated by Lya and so oscillated.
479  * especially bad in primal.in - change 2 to 4 and primal.in will stop due to Lya oscil */
482  /* >>chng 02 feb 16, add damper on H- photo rate, wild oscillations in Lya photo rate in
483  * grain free models */
484 
485  /*hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , nhe1Com.nhe1[0] , opac.iphmop ,
486  0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling );*/
489 
490  /* save H- photodissociation heating */
492 
493  {
494  /* following should be set true to print populations */
495  /*@-redef@*/
496  enum {DEBUG_LOC=false};
497  /*@+redef@*/
498  if( DEBUG_LOC)
499  {
500  fprintf(ioQQQ,"hminphoto\t%li\t%li\t%.2e\n", nzone, conv.nPres2Ioniz , hmi.HMinus_photo_rate );
501  }
502  }
503 
504  /* induced recombination */
506 
507  /* induced recombination cooling per unit volume */
509  hmi.HMinus_induc_rec_cooling *= hmi.rel_pop_LTE_Hmin*dense.eden*hmi.Hmolec[ipMH]; /* dense.xIonDense[ipHYDROGEN][0]; */
510 
511  {
512  /* following should be set true to debug H- photoionization rates */
513  /*@-redef@*/
514  enum {DEBUG_LOC=false};
515  /*@+redef@*/
516  if( DEBUG_LOC && nzone>400/*&& iteration > 1*/)
517  {
518  fprintf(ioQQQ,"hmoledebugg %.2f ",fnzone);
519  GammaPrt(
520  hmi.iphmin-1 , iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] , opac.iphmop ,
521  /* io unit we will write to */
522  ioQQQ,
523  /* total photo rate from previous call to GammaK */
525  /* we will print contributors that are more than this rate */
526  hmi.HMinus_photo_rate*0.05);
527  }
528  }
529  /* add on high energy ionization, assume hydrogen cross section
530  * n.b.; HGAMNC with secondaries */
531  /* >>chng 00 dec 24, above goes to HeI edge, no need for this, and was not important*/
532  /*hmi.HMinus_photo_rate += iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s];*/
533 
534  /* ================================================================= */
535  /* photodissociation by Lyman band absorption: esc prob treatment,
536  * treatment based on
537  * >>refer HI abs Tielens & Hollenbach 1985 ApJ 291, 722. */
538  /* do up to carbon photo edge if carbon is turned on */
539  /* >>>chng 00 apr 07, add test for whether element is turned on */
541  /* >>chng 00 apr 07 from explicit ipHeavy to ipLo */
542  /* find total intensity over carbon-ionizing continuum */
543  /* >>chng 03 jun 09, use exact bounds rather than CI photo threshold for lower bound */
544  /*for( i=ipLo; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )*/
545  /* the integral is from 6eV to 13.6, or 2060 - 912 Ang */
546  for( i=rfield.ipG0_TH85_lo; i < rfield.ipG0_TH85_hi; ++i )
547  {
549  rfield.outlin[i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
550  }
551 
552  /* now convert to Habing ISM units
553  * UV_Cont_rel2_Habing_TH85_face is FUV continuum relative to Habing value
554  * 1.6e-3 ergs em-2 s-1 is the Habing 1968 field, quoted on page 723, end of first
555  * full paragraph on left */
557  /* if start of calculation remember G0 at illuminated face */
558  if( nzone<=1 )
559  {
561  }
562 
563 
564  /* >>chng 05 jan 09, add special ver of G0 */
566  for( i=rfield.ipG0_spec_lo; i < rfield.ipG0_spec_hi; ++i )
567  {
569  rfield.outlin[i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
570  }
572 
573  /* the Draine & Bertoldi version of the same thing, defined over their energy range */
574  /* >>chng 04 feb 07, only evaluate at the illuminated face */
576  {
577  /* this is sum of photon number between 912A and 1110, as per BD96 */
578  for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; ++i )
579  {
581  rfield.outlin[i-1]+ rfield.outlin_noplot[i-1]);
582  }
583  /* convert into scaled ISM background field, total number of photons over value for 1 ISM field,
584  * the coefficient 1.232e7 is the number of photons over this wavelength range for 1H and is
585  * given in BD96, page 225, 4th line from start of section 4, also page 272, table 1, 2nd line
586  * from bottom */
587  /* >>chng 04 mar 16, introduce the 1.71 */
588  /* equation 20 of DB96 gives chi as flux over 1.21e7, to produce one Habing field.
589  * to get the Draine field we need to further divide by 1.71 as stated on the first
590  * line after equation 23 */
592  }
593 
594  /* escape prob takes into account line shielding,
595  * next is opacity then optical depth in H2 UV lines, using eqn A7 of TH85,
596  * LIMELM+2 points to H2 */
597  hmi.H2Opacity = (realnum)(1.2e-14*(1e5/DoppVel.doppler[LIMELM+2]));
598  /* the typical Lyman -Werner H2 line optical depth eq A7 of TH85a */
600  /* the escape probability - chance that continuum photon will penetrate to
601  * this depth to pump the Lyman Werner bands */
602  h2esc = esc_PRD_1side(th2,1e-4);
603 
604  /* cross section is eqn A8 of
605  * >>refer H2 dissoc Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722
606  * branching ratio of 10% is already included, so 10x smaller than their number
607  * 10% lead to dissociation through H_2 + h nu => 2H */
608  /* >>chng 05 mar 10, by TE, break into 2g and 2s
609  * note use of same shielding column in below - can do far better */
613 
614  /* these are Burton et al. 1990 rates */
618 
619  {
620  /* the following implements Drain & Bertoldi 1996, equation 37 from
621  * >>refer H2 destruction Draine, B.T., & Bertoldi, F., 1996, ApJ, 468, 269-289
622  * but the constant 4.6e-11 comes from Bertoldi & Draine equation 23,
623  * this is the normalized H2 column density */
624  double x = (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) / 5e14;
625  double sqrtx = sqrt(1.+x);
626  /* Doppler with of H2 in km/s */
627  double b5 = DoppVel.doppler[LIMELM+2]/1e5;
628  /* the molecular hydrogen line self-shielding factor */
629  double fshield = 0.965/POW2(1.+x/b5) + 0.035/sqrtx *
630  sexp(8.5e-4*sqrtx);
631 
632  /*double fshield = pow( MAX2(1.,colden.colden[ipCOLH2]/1e14) , -0.75 );*/
633  /* this is the Bertoldi & Draine version of the Habing field,
634  * with dilution of radiation and extinction due to grains */
635  /* >>chng 04 apr 18, moved fshield, the line shielding factor, from this defn to
636  * the following defn of dissociation rate, since following should measure continuum */
639 
640  /* the following comes from Bertoldi & Draine 1996, equation 23,
641  * hmi.UV_Cont_rel2_Draine_DB96_depth already includes a factor of 1.71 to correct back from TH85 */
642  /* >>chng 05 mar 10, TE, break into 2s and 2s */
643  if( co.lgUMISTrates )
644  {
645  /* this is default, when set Leiden hack UMIST rates not entered */
648  }
649  else
650  {
651  /* done when set Leiden hack UMIST rates command entered */
653  *sexp(3.02*rfield.extin_mag_V_point)* fshield;
655  *sexp(3.02*rfield.extin_mag_V_point)* fshield;
656  }
657 
658  /* BD do not give an excitation rate, so used 9 times the dissociation
659  * rate by analogy with 90% going back into X, as per TH85 */
660  /*>>chng 04 feb 07, had used 90% relax into X from TH85,
661  * BD say that 15% dissociate, so 85/15 = 5.67 is factor */
663  }
664 
665  /* do Elwert approximation to the dissociation rate */
667  {
668  /* this evaluates the new H2 dissociation rate by Torsten Elwert */
669  /* chng 05 jun 23, TE
670  * >>chng 05 sep 13, update master source with now approximation */
671 
672  /* Doppler with of H2 in km/s */
673  double b5 = DoppVel.doppler[LIMELM+2]/1e5;
674 
675  /* split the Solomon rates in H2g and H2s */
676  /* >>chng 05 oct 19,
677  * >>chng 05 dec 05, TE, define new approximation for the heating due to the destruction of H2
678  * use this approximation for the specified cloud parameters, otherwise
679  * use limiting cases for 1 <= G0, G0 >= 1e7, n >= 1e7, n <= 1 */
680 
681  double x_H2g, x_H2s,
682  fshield_H2g, fshield_H2s,
683  f_H2s;
684  static double a_H2g, a_H2s,
685  e1_H2g, e1_H2s,
686  e2_H2g,
687  b_H2g,
688  sl_H2g, sl_H2s,
689  k_f_H2s,
690  k_H2g_to_H2s,
691  log_G0_face = -1;
692 
693  /* define parameter range for the new approximation
694  * test for G0
695  *BUGFIX - this tested on lg_G0_face < 0 for initialization needed so did not work
696  * in grids - change to evaluate in zone 0 */
697  /* >>chng 07 feb 24, BUGFIX crash when G0=0 at start and radiation
698  * field builds up due to diffuse fields - soln is to always reevaluate */
699  /*if( !nzone )*/
700  {
702  {
703  log_G0_face = 0.;
704  }
705  else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
706  {
707  log_G0_face = 7.;
708  }
709  else
710  {
711  log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
712  }
713 
714  /* terms only dependent on G0_face */
715 
716  /* coefficients and exponents */
717  a_H2g = 0.06 * log_G0_face + 1.32;
718  a_H2s = 0.375 * log_G0_face + 2.125;
719 
720  e1_H2g = -0.05 * log_G0_face + 2.25;
721  e1_H2s = -0.125 * log_G0_face + 2.625;
722 
723  e2_H2g = -0.005 * log_G0_face + 0.625;
724 
725  b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
726 
727  /* scalelength for H2g and H2s */
728  sl_H2g = 4.0e14;
729  sl_H2s = 9.0e15;
730 
731  /* coefficient for 2nd term of Solomon H2s */
732  k_f_H2s = MAX2(0.1,2.375 * log_G0_face - 1.875 );
733 
734  /* coefficient for branching ratio */
735  k_H2g_to_H2s = MAX2(1.,-1.75 * log_G0_face + 11.25);
736 
737  /*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",
738  e1_H2g, e1_H2s, e2_H2g, b_H2g, a_H2g, a_H2s, sl_H2g, sl_H2s);
739  */
740  }
741 
742  /* Solomon H2s ~G0^0.2 at large depth*/
743  f_H2s = k_f_H2s * pow((double)hmi.UV_Cont_rel2_Draine_DB96_depth, 0.2 );
744 
745  /* scale length for absorption of UV lines */
746  x_H2g = (colden.colden[ipCOL_H2g]) / sl_H2g;
747  x_H2s = (colden.colden[ipCOL_H2s]) / sl_H2s;
748 
749  /* the molecular hydrogen line self-shielding factor */
750  fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
751  fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
752 
753  /* the Elwert Solomon rates for H2g and H2s hmi.chH2_small_model_type == 'E' */
755  hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = a_H2s * 4.6e-11 * fshield_H2s * (hmi.UV_Cont_rel2_Draine_DB96_depth + f_H2s);
756 
757  /* assume branching ratio dependent on G0*/
759 
760  /* use G0_BD96 as this definition declines faster with depth which is physical as
761  * the longer wavelengths in the definition of G0_TH85 cannot dissociate
762  * H2s directly */
765  }
766  else
767  {
772  }
773 
774  /* this is rate of photodissociation of H2*, A12 of TH85 */
776 
777  /* dissociation rate from Burton et al. 1990 */
779 
780  /* rates for cosmic ray excitation of singlet bound electronic bound excited states
781  * only add this to small molecule since automatically included in large
782  *>>refer H2 cr excit Dalgarno, A., Yan, Min, & Liu, Weihong 1999, ApJS, 125, 237
783  * this is excitation of H2* */
784  /* >>chng 05 sep 13, do not include this process when Leiden hacks are in place */
785  cr_H2s = secondaries.x12tot*0.9 / 2. * hmi.lgLeidenCRHack;
786  /* this is the fraction that dissociate */
787  /* >>chng 05 sep 13, do not include this process when Leiden hacks are in place */
788  cr_H2dis = secondaries.x12tot*0.1 / 2. * hmi.lgLeidenCRHack;
789 
790  /* >>chng 05 sep 13, TE, synchronize treatment of CR */
791  /* cosmic ray rates for dissociation of ground and H2s
792  * two factors done to agree with large H2 deep in the cloud where
793  * cosmic rays are important */
794  cr_H2dis_ELWERT_H2g = secondaries.x12tot*5e-8 * hmi.lgLeidenCRHack;
795  cr_H2dis_ELWERT_H2s = secondaries.x12tot*4e-2 * hmi.lgLeidenCRHack;
796 
797  /* at this point there are two or three independent estimates of the H2 dissociation rate.
798  * if the large H2 molecule is on, then H2 Solomon rates has been defined in the last
799  * call to the large molecule. Just above we have defined hmi.H2_Solomon_dissoc_rate_TH85,
800  * the dissociation rate from Tielens & Hollenbach 1985, and hmi.H2_Solomon_dissoc_rate_BD96,
801  * the rate from Bertoldi & Draine 1996. We can use any defined rate. If the big H2
802  * molecule is on, use its rate. If not, for now use the TH85 rate, since that is the
803  * rate we always have used in the past.
804  * The actual rate we will use is given by hmi.H2_Solomon_dissoc_rate_used
805  */
806  /* this is the Solomon process dissociation rate actually used */
808  {
809  /* only update after big H2 molecule has been evaluated,
810  * when very little H2 and big molecule not used, leave at previous (probably TH85) value,
811  * since that value is always known */
812 
813  /* Solomon process rate from X into the X continuum with units s-1
814  * rates are total rate, and rates from H2g and H2s */
817  {
818  /*fprintf(ioQQQ,"DEBUG his sol dis <0\n");*/
820  }
821 
825 
826  /* photoexcitation from H2g to H2s */
828  if( hmi.H2_H2g_to_H2s_rate_used <= 0. )
830 
831  /* add up H2s + hnu (continuum) => 2H + KE, continuum photodissociation,
832  * this is not the Solomon process, true continuum, units s-1 */
833  /* only include rates from H2s since this is only open channel, this process is well
834  * shielded against Lyman continuum destruction by atomic hydrogen */
836  if( hmi.H2_photodissoc_used_H2s <= 0. )
838  /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
839  * for unfavorable wavelength range of G0*/
841  if( hmi.H2_photodissoc_used_H2g <= 0. )
843  }
844  else if( hmi.chH2_small_model_type == 'T' )
845  {
846  /* the TH85 rate */
847  /*>>chng 05 jun 23, add cosmic rays */
849  /* >>chng 05 sep 13, cr_H2dis was not included */
852 
853  /* continuum photodissociation H2s + hnu => 2H, ,
854  * this is not the Solomon process, true continuum, units s-1 */
856  /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
857  * for unfavorable wavelength range of G0*/
859  }
860 
861  else if( hmi.chH2_small_model_type == 'H' )
862  {
863  /* the improved H2 formalism given by
864  *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M
865  >>refcon 1990, ApJ, 365, 620 */
867  /* >>chng 05 sep 13, cr_H2dis was not included */
870 
871  /* continuum photodissociation H2s + hnu => 2H, ,
872  * this is not the Solomon process, true continuum, units s-1 */
874  /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
875  * for unfavorable wavelength range of G0*/
877  }
878 
879  else if( hmi.chH2_small_model_type == 'B' )
880  {
881  /* the Bertoldi & Draine rate - this is the default */
882  /*>>chng 05 jun 23, add cosmic rays */
884  /* >>chng 05 sep 13, cr_H2dis was not included */
886  /* they did not do the excitation or dissoc rate, so use TH85 */
888 
889 
890  /* continuum photodissociation H2s + hnu => 2H, ,
891  * this is not the Solomon process, true continuum, units s-1 */
893  /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
894  * for unfavorable wavelength range of G0*/
896  }
897  else if( hmi.chH2_small_model_type == 'E' )
898  {
899  /* the Elwert et al. rate
900  *>>chng 05 jun 23, add cosmic rays */
904 
905 
906  /* continuum photodissociation H2s + hnu => 2H, ,
907  * this is not the Solomon process, true continuum, units s-1 */
910  }
911  else
912  TotalInsanity();
913 
914  {
915  /*@-redef@*/
916  enum {DEBUG_LOC=false};
917  /*@+redef@*/
918  if( DEBUG_LOC && h2.lgH2ON )
919  {
920  fprintf(ioQQQ," Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
926  }
927  }
928  return;
929 }
930 

Generated for cloudy by doxygen 1.8.3.1