cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_h2.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 /*H2_ContPoint set the ipCont struc element for the H2 molecule, called by ContCreatePointers */
4 /*H2_Accel radiative acceleration due to H2 */
5 /*H2_RadPress rad pressure due to h2 lines called in PresTotCurrent */
6 /*H2_InterEnergy internal energy of H2 called in PresTotCurrent */
7 /*H2_RT_diffuse do emission from H2 - called from RT_diffuse */
8 /*H2_itrzn - average number of H2 pop evaluations per zone */
9 /*H2_RTMake do RT for H2 - called from RT_line_all */
10 /*H2_RT_tau_inc increment optical depth for the H2 molecule, called from RT_tau_inc */
11 /*H2_LineZero initialize optical depths in H2, called from RT_tau_init */
12 /*H2_RT_tau_reset the large H2 molecule, called from RT_tau_reset */
13 /*H2_Colden maintain H2 column densities within X */
14 /*H2_LevelPops do level H2_populations for H2, called by Hydrogenic */
15 /*H2_Level_low_matrix evaluate CO rotation cooling */
16 /*H2_cooling evaluate cooling and heating due to H2 molecule */
17 /*H2_X_coll_rate_evaluate find collisional rates within X */
18 /*cdH2_colden return column density in H2, negative -1 if cannot find state,
19  * header is cddrive */
20 /*H2_DR choose next zone thickness based on H2 big molecule */
21 /* turn this flag on to do minimal debug print of pops */
22 #define PRT_POPS false
23 /* this is limit to number of loops over H2 pops before giving up */
24 #define LIM_H2_POP_LOOP 100
25 /* this is a typical dissociation cross section (cm2) for H2 + Hnu -> 2H + ke */
26 /* >>chng 05 may 11, had been 2.5e-19 */
27 #define H2_DISS_ALLISON_DALGARNO 6e-19f
28 #include "cddefines.h"
29 #include "cddrive.h"
30 #include "physconst.h"
31 #include "taulines.h"
32 #include "atoms.h"
33 #include "conv.h"
34 #include "secondaries.h"
35 #include "pressure.h"
36 #include "trace.h"
37 #include "hmi.h"
38 #include "hextra.h"
39 #include "rt.h"
40 #include "radius.h"
41 #include "ipoint.h"
42 #include "phycon.h"
43 #include "thermal.h"
44 #include "dense.h"
45 #include "rfield.h"
46 #include "lines_service.h"
47 #include "mole.h"
48 #include "h2.h"
49 #include "h2_priv.h"
50 
51 /* this counts how many times we go through the H2 level H2_populations loop */
52 static long int loop_h2_pops;
53 
54 realnum H2_te_hminus[nTE_HMINUS] = {10.,30.,100.,300.,1000.,3000.,10000.};
55 
56 /* this will contain a vector for collisions within the X ground electronic state,
57  * CollRateFit[vib_up][rot_up][vib_lo][rot_lo][coll_type][3] */
60 
61 /* this integer is added to rotation quantum number J for the test of whether
62  * a particular J state is ortho or para - the state is ortho if J+below is odd,
63  * and para if J+below is even */
64 int H2_nRot_add_ortho_para[N_H2_ELEC] = {0 , 1 , 1 , 0, 1, 1 , 0};
65 
66 /* dissociation energies (cm-1) for each electronic state, from
67  * >>refer H2 energies Sharp, T. E., 1971, Atomic Data, 2, 119
68  * energy for H_2 + hnu -> 2H,
69  * note that energies for excited states are in the Lyman continuum
70  * (1 Ryd = 109670 cm-1) */
71 /* >>chng 02 oct 08, improved energies */
73 { 36118.11, 118375.6, 118375.6, 118375.6, 118375.6,133608.6,133608.6 };
74 /* original values
75  { 36113., 118372., 118372., 118372., 118372.,0.,0. };*/
76 
77 double exp_disoc;
78 
79 /*H2_X_coll_rate_evaluate find collisional rates within X -
80  * this is one time upon entry into H2_LevelPops */
82 {
83  long int nColl,
84  ipHi ,
85  ipLo,
86  ip,
87  iVibHi,
88  iRotHi,
89  iVibLo,
90  iRotLo;
91  realnum colldown;
92  double exph2_big,
93  exph2_big_0,
94  rel_pop_LTE_H2_big;
95 
96  DEBUG_ENTRY( "H2_X_coll_rate_evaluate()" );
97 
98  /* set collider density
99  * the colliders are:
100  * [0] = H
101  * [1], [5] = He (old and new cs data)
102  * [2] = H2 ortho
103  * [3] = H2 para
104  * [4] = H+ + H3+ */
105  /* atomic hydrogen */
107  /* all ortho h2 */
108  /* He - H2 */
110  /* H2 - H2(ortho) */
112  /* all para H2 */
114  /* protons - ionized hydrogen */
116  /* H3+ - assume that H3+ has same rates as proton */
118 
120 
121  /* this is total density of all colliders, is only used for collisional dissociation
122  * rates for H2 are not included here, will be added separately*/
125  (realnum)dense.eden;
126 
128  {
129  fprintf(ioQQQ," Collider densities are:");
130  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
131  {
132  fprintf(ioQQQ,"\t%.3e", collider_density[nColl]);
133  }
134  fprintf(ioQQQ,"\n");
135  }
136 
137  for( ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
138  {
139  H2_X_source[ipHi] = 0.;
140  H2_X_sink[ipHi] = 0.;
141  }
142  H2_X_coll_rate.zero();
143 
144  /*>>chng 05 sep 18, GS, determine LTE populations for H2+ e = H- + H*/
146  exph2_big_0 = exp_disoc/SDIV(H2_Boltzmann[0][0][0]);
147  /*>>chng 05 oct 17, GS, (2*m_e/m_p)^3/2 = 3.634e-5 */
148  rel_pop_LTE_H2_big =SAHA/SDIV((phycon.te32*exph2_big_0))*(H2_stat[0][0][0]/(2.*2.))*3.634e-5;
149  /* now find rates for all collisions within X */
150  /* this is special J=1 to J=0 collision, which is only fast at
151  * very low grain temperatures
152  H2_X_coll_rate[1][0] =
153  (realnum)(hmi.rate_grain_h2_J1_to_J0);*/
154 
155  /* count formation from grains and H- as a collisional formation process */
156  /* cm-3 s-1, evaluated in mole_H2_form */
157  H2_X_source[0] += H2_X_formation[0][0];
158  /*>>chng 05 sep 18, GS, H2+ e = H- + H, unit s-1
159  * H2_X_Hmin_back[iVib][iRot] is resolved formation rate, units cm3 s-1 */
161  SDIV(rel_pop_LTE_H2_big)*dense.eden);
162  /* this represents collisional dissociation into continuum of X,
163  * rates are just guesses */
166  /*>>chng 05 jul 20, GS, collisional dissociation with H2g and H2s are added here*/
167  H2_X_sink[0] += hmi.H2_total*
169 
170  /* this is cosmic ray or secondary photodissociation into mainly H2+ */
172 
173  /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H- */
175 
176  /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H + e */
178 
179  /*>>chng 05 sep 26, GS, H2 + CR -> H+ + H + e */
180  H2_X_sink[0] += (realnum)(secondaries.csupra[ipHYDROGEN][0]*0.0478);
181 
182  /* collisional dissociation by non-thermal electrons
183  * and cosmic rays to the triplet state of H2 (a,b,c ) from
184  *>>>refer Dalgarno, Yan, & Liu 1999 ApJs
185  * rates depend on energy of secondary electrons, this is an estimate
186  * from fig 4b around incident electron energy 15 to 20 eV.
187  * The cross section of state b is divided by the cross-section of Lya
188  * (the ratio is ~ 10)and then multiplied by the cosmic ray excitation
189  * rate of Lya (secondaries.x12tot) the triplet state of H2 (b )
190  *>>chng 07 apr 08, GS from 3 to 10 after study of Dalgarno et al */
191  H2_X_sink[0] += 10.f*secondaries.x12tot;
192 
193  /*>>chng 05 jul 07, GS, ionization by hard photons are added to be consistent with chemistry-network */
195 
196  /* rate (s-1) out of this level, this is dissociation into the continuum of singlet B and C states*/
198 
199  /*>>chng 05 sept 28, GS, H2 + H+ = H2+ + H; note that H2 is destroyed from v=0, J=0,
200  * forward and backward reactions are not in detailed balance, astro-ph/0404288, final unit s-1
201  * ihmi.rh2h2p is a rate coefficient, cm3 s-1, and needs a density to become
202  * the sink inverse lifetime. The process is H2(v=0, J=0) + H+ -> H0 + H2+ */
204 
205  /* now find total coll rates - this loop is over the energy-sorted levels */
206  ipHi = nLevels_per_elec[0];
207  while( (--ipHi) > 0 )
208  {
209  /* array of energy sorted indices within X */
210  ip = H2_ipX_ener_sort[ipHi];
211  iVibHi = ipVib_H2_energy_sort[ip];
212  iRotHi = ipRot_H2_energy_sort[ip];
213 
214  /*>>chng 05 sep 18, GS, determine LTE populations for H2+ e = H- + H, unit s-1*/
215  exph2_big = exp_disoc/SDIV(H2_Boltzmann[0][iVibHi][iRotHi]);
216  /* >>chng 05 oct 13, replace 1 with stat wght of level */
217  /*>>chng 05 oct 17, GS, (2*m_e/m_p)^3/2 = 3.634e-5 */
218  rel_pop_LTE_H2_big = SAHA/SDIV((phycon.te32*exph2_big))*(H2_stat[0][iVibHi][iRotHi]/(2.*2.))*3.634e-5;
219  /* formation */
220  /* count formation from grains and H- as a collisional formation process */
221  /* cm-3 s-1, evaluated in mole_H2_form */
222  H2_X_source[ipHi] += H2_X_formation[iVibHi][iRotHi];
223  /*>>chng 05 sep 18, GS, H2+ e = H- + H*, H2_X_Hmin_back has units cm3 s-1 so final units s-1 */
224  H2_X_sink[ipHi] += (realnum)(H2_X_Hmin_back[iVibHi][iRotHi]*hmi.rel_pop_LTE_Hmin/
225  SDIV(rel_pop_LTE_H2_big)*dense.eden);
226  /* this represents collisional dissociation into continuum of X,
227  * rates are just guesses */
230 
231  /*>>chng 05 jul 20, GS, collisional dissociation with H2g and H2s are added here*/
232  H2_X_sink[ipHi] += hmi.H2_total*
234 
235  /* this is cosmic ray or secondary photodissociation into mainly H2+ */
237 
238  /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H- */
240 
241  /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H + e */
243 
244  /*>>chng 05 sep 26, GS, H2 + CR -> H+ + H + e */
245  H2_X_sink[ipHi] += (realnum)(secondaries.csupra[ipHYDROGEN][0]*0.0478);
246 
247  /* collisional dissociation by non-thermal electrons
248  * and cosmic rays to the triplet state of H2 (a,b,c ) from
249  *>>>refer Dalgarno, Yan, & Liu 1999 ApJs
250  * rates depend on energy of secondary electrons, this is an estimate
251  *from fig 4b around incident electron energy 15 to 20 eV.
252  * The cross section of state b is divided by the cross-section of Lya
253  * (the ratio is ~ 10)and then multiplied by the cosmic ray excitation
254  * rate of Lya (secondaries.x12tot) the triplet state of H2 (b )
255  *>>chng 07 apr 08, GS from 3 to 10 after study of Dalgarno et al */
256  H2_X_sink[ipHi] += 10.f*secondaries.x12tot;
257 
258  /*>>chng 05 jul 07, GS, ionization by hard photons are added to be consistent with chemistry-network */
260 
261  /* rate (s-1) out of this level */
262  H2_X_sink[ipHi] += rfield.flux_accum[H2_ipPhoto[iVibHi][iRotHi]-1]*H2_DISS_ALLISON_DALGARNO;
263 
265  {
266  /* excitation within X due to thermal particles */
267  for( ipLo=0; ipLo<ipHi; ++ipLo )
268  {
269  /* array of energy sorted indices within X */
270  ip = H2_ipX_ener_sort[ipLo];
271  iVibLo = ipVib_H2_energy_sort[ip];
272  iRotLo = ipRot_H2_energy_sort[ip];
273 
274  /* collisional interactions with upper levels within X */
275  colldown = 0.;
276  mr5ci H2CollRate = H2_CollRate.begin(iVibHi,iRotHi,iVibLo,iRotLo);
277  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
278  {
279  /* downward collision rate, units s-1 */
280  colldown += H2CollRate[nColl]*collider_density[nColl];
281  ASSERT( H2CollRate[nColl]*collider_density[nColl] >= 0. );
282  }
283  /* rate in from upper level, units cm-3 s-1 */
284  H2_X_coll_rate[ipHi][ipLo] += colldown;
285  }/* end loop over ipLo */
286  }
287  }/* end loop over ipHi */
288  return;
289 }
290 
291 /*H2_itrzn - average number of H2 pop evaluations per zone */
292 double H2_itrzn( void )
293 {
294  if( h2.lgH2ON && nH2_zone>0 )
295  {
296  return( (double)nH2_pops / (double)nH2_zone );
297  }
298  else
299  {
300  return 0.;
301  }
302 }
303 
304 /* set the ipCont struc element for the H2 molecule, called by ContCreatePointers */
305 void H2_ContPoint( void )
306 {
307  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
308 
309  if( !h2.lgH2ON )
310  return;
311 
312  DEBUG_ENTRY( "H2_ContPoint()" );
313 
314  /* set array index for line energy within continuum array */
315  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
316  {
317  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
318  {
319  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
320  {
321  /* now the lower levels */
322  /* NB - X is the only lower level considered here, since we are only
323  * concerned with excited electronic levels as a photodissociation process
324  * code exists to relax this assumption - simply change following to iElecHi */
325  long int lim_elec_lo = 0;
326  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
327  {
328  /* want to include all vibration states in lower level if different electronic level,
329  * but only lower vibration levels if same electronic level */
330  long int nv = h2.nVib_hi[iElecLo];
331  if( iElecLo==iElecHi )
332  nv = iVibHi;
333  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
334  {
335  long nr = h2.nRot_hi[iElecLo][iVibLo];
336  if( iElecLo==iElecHi && iVibHi==iVibLo )
337  nr = iRotHi-1;
338 
339  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
340  {
341  /* >>chng 05 feb 07, use lgH2_line_exists */
342  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
343  {
344  ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul > 0. );
345  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont =
346  ipLineEnergy(
347  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD ,
348  "H2 " , 0 );
349  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ipFine =
350  ipFineCont(
351  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD );
352  }
353  }
354  }
355  }
356  }
357  }
358  }
359  return;
360 }
361 
362 /* ===================================================================== */
363 /* radiative acceleration due to H2 called in rt_line_driving */
364 double H2_Accel(void)
365 {
366  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
367  double h2_drive;
368 
369  /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */
370  if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ )
371  return 0.;
372 
373  DEBUG_ENTRY( "H2_Accel()" );
374 
375  /* this routine computes the line driven radiative acceleration
376  * due to H2 molecule*/
377 
378  h2_drive = 0.;
379  /* loop over all possible lines */
380  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
381  {
382  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
383  {
384  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
385  {
386  /* now the lower levels */
387  /* NB - X is the only lower level considered here, since we are only
388  * concerned with excited electronic levels as a photodissociation process
389  * code exists to relax this assumption - simply change following to iElecHi */
390  long int lim_elec_lo = 0;
391  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
392  {
393  /* want to include all vibration states in lower level if different electronic level,
394  * but only lower vibration levels if same electronic level */
395  long int nv = h2.nVib_hi[iElecLo];
396  if( iElecLo==iElecHi )
397  nv = iVibHi;
398  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
399  {
400  long nr = h2.nRot_hi[iElecLo][iVibLo];
401  if( iElecLo==iElecHi && iVibHi==iVibLo )
402  nr = iRotHi-1;
403 
404  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
405  mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
406  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
407  {
408  /* >>chng 03 feb 14, from !=0 to >0 */
409  /* >>chng 05 feb 07, use lgH2_line_exists */
410  if( *lgH2le++ )
411  {
412  ASSERT( H2L->ipCont > 0 );
413  h2_drive += H2L->Emis->pump*H2L->EnergyErg*H2L->Emis->PopOpc;
414  }
415  ++H2L;
416  }
417  }
418  }
419  }
420  }
421  }
422  return h2_drive;
423 }
424 
425 /* ===================================================================== */
426 /* rad pressure due to H2 lines called in PresTotCurrent */
427 double H2_RadPress(void)
428 {
429  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
430  double press;
431 
432  /* will be used to check on size of opacity, was capped at this value */
433  realnum smallfloat=SMALLFLOAT*10.f;
434 
435  /* radiation pressure sum is expensive - do not evaluate if we did not
436  * bother evaluating large molecule */
437  if( !h2.lgH2ON || !h2.nCallH2_this_zone )
438  return 0.;
439 
440  DEBUG_ENTRY( "H2_RadPress()" );
441 
442  press = 0.;
443  /* loop over all possible lines */
444  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
445  {
446  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
447  {
448  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
449  {
450  /* now the lower levels - X is the only lower level considered
451  * here, we only do excited electronic levels as a
452  * photodissociation process - code exists to relax this
453  * assumption - simply change following to iElecHi */
454  long int lim_elec_lo = 0;
455  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
456  {
457  /* want to include all vibration states in lower level if different electronic level,
458  * but only lower vibration levels if same electronic level */
459  long int nv = h2.nVib_hi[iElecLo];
460  if( iElecLo==iElecHi )
461  nv = iVibHi;
462  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
463  {
464  long nr = h2.nRot_hi[iElecLo][iVibLo];
465  if( iElecLo==iElecHi && iVibHi==iVibLo )
466  nr = iRotHi-1;
467 
468  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
469  mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
470  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
471  {
472  if( *lgH2le++ )
473  {
474  ASSERT( H2L->ipCont > 0 );
475 
476  if( H2L->Hi->Pop > smallfloat && H2L->Emis->PopOpc > smallfloat )
477  {
478  double RadPres1 = PressureRadiationLine( &(*H2L), 1.0 );
479  press += RadPres1;
480  }
481  }
482  ++H2L;
483  }
484  }
485  }
486  }
487  }
488  }
489 
491  fprintf(ioQQQ,
492  " H2_RadPress returns, radiation pressure is %.2e\n",
493  press );
494  return press;
495 }
496 
497 #if 0
498 /* ===================================================================== */
499 /* internal energy of H2 called in PresTotCurrent */
500 double H2_InterEnergy(void)
501 {
502  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
503  double energy;
504 
505  /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */
506  if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ )
507  return 0.;
508 
509  DEBUG_ENTRY( "H2_InterEnergy()" );
510 
511  broken(); /* the code below contains serious bugs! It is supposed to implement a loop
512  * over all states in order to evaluate the potential energy stored in
513  * excited states of H2. The code below however implements loops over all
514  * combinations of upper and lower levels! */
515 
516  energy = 0.;
517  /* loop over all possible lines */
518  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
519  {
520  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
521  {
522  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
523  {
524  double H2popHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop;
525 
526  /* now the lower levels
527  * NB - X is the only lower level considered here, since we are only
528  * concerned with excited electronic levels as a photodissociation process
529  * code exists to relax this assumption - simply change following to iElecHi */
530  long int lim_elec_lo = 0;
531  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
532  {
533  /* want to include all vibration states in lower level if different electronic level,
534  * but only lower vibration levels if same electronic level */
535  long int nv = h2.nVib_hi[iElecLo];
536  if( iElecLo==iElecHi )
537  nv = iVibHi;
538  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
539  {
540  long nr = h2.nRot_hi[iElecLo][iVibLo];
541  if( iElecLo==iElecHi && iVibHi==iVibLo )
542  nr = iRotHi-1;
543 
544  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
545  mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
546  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
547  {
548  /* >>chng 03 feb 14, from !=0 to >0 */
549  /* >>chng 05 feb 07, use lgH2_line_exists */
550  if( *lgH2le++ )
551  {
552  energy += H2popHi*H2L->EnergyErg;
553  }
554  ++H2L;
555  }
556  }
557  }
558  }
559  }
560  }
561  return energy;
562 }
563 #endif
564 
565 /*H2_RT_diffuse do emission from H2 - called from RT_diffuse */
566 void H2_RT_diffuse(void)
567 {
568  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
569 
570  if( !h2.lgH2ON || !h2.nCallH2_this_zone )
571  return;
572 
573  DEBUG_ENTRY( "H2_RT_diffuse()" );
574 
575  /* loop over all possible lines */
576  /* NB - this loop does not include the electronic lines */
577  for( iElecHi=0; iElecHi<1; ++iElecHi )
578  {
579  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
580  {
581  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
582  {
583  /* now the lower levels */
584  /* NB - X is the only lower level considered here, since we are only
585  * concerned with excited electronic levels as a photodissociation process
586  * code exists to relax this assumption - simply change following to iElecHi */
587  long int lim_elec_lo = 0;
588  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
589  {
590  /* want to include all vibration states in lower level if different electronic level,
591  * but only lower vibration levels if same electronic level */
592  long int nv = h2.nVib_hi[iElecLo];
593  if( iElecLo==iElecHi )
594  nv = iVibHi;
595  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
596  {
597  long nr = h2.nRot_hi[iElecLo][iVibLo];
598  if( iElecLo==iElecHi && iVibHi==iVibLo )
599  nr = iRotHi-1;
600 
601  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
602  {
603  /* >>chng 03 feb 14, from !=0 to > 0 */
604  /* >>chng 05 feb 07, use lgH2_line_exists */
605  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
606  {
607  outline( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
608  }
609  }
610  }
611  }
612  }
613  }
614  }
615  return;
616 }
617 
618 /* do RT for H2 - called from RT_line_all */
619 void H2_RTMake( bool lgDoEsc , bool lgUpdateFineOpac )
620 {
621  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
622 
623  if( !h2.lgH2ON )
624  return;
625 
626  DEBUG_ENTRY( "H2_RTMake()" );
627 
628  /* this routine drives calls to make RT relations for H2 molecule */
629  /* loop over all possible lines */
630  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
631  {
632  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
633  {
634  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
635  {
636  /* now the lower levels */
637  /* NB - X is the only lower level considered here, since we are only
638  * concerned with excited electronic levels as a photodissociation process
639  * code exists to relax this assumption - simply change following to iElecHi */
640  long int lim_elec_lo = 0;
641  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
642  {
643  /* want to include all vibration states in lower level if different electronic level,
644  * but only lower vibration levels if same electronic level */
645  long int nv = h2.nVib_hi[iElecLo];
646  if( iElecLo==iElecHi )
647  nv = iVibHi;
648  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
649  {
650  long nr = h2.nRot_hi[iElecLo][iVibLo];
651  if( iElecLo==iElecHi && iVibHi==iVibLo )
652  nr = iRotHi-1;
653 
654  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
655  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
656  {
657  /* >>chng 03 feb 14, change test from !=0 to >0 */
658  /* >>chng 05 feb 07, use lgH2_line_exists */
659  if( *lgH2le++ )
660  {
661  /* >>chng 03 jun 18, added 4th parameter in call to this routine - says to not
662  * include self-shielding of line across this zone. This introduces a dr dependent
663  * variation in the line pumping rate, which made H2 abundance fluctuate due to
664  * Solomon process having slight dr-caused mole. */
665  RT_line_one( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo], lgDoEsc, lgUpdateFineOpac, false );
666  }
667  }
668  }
669  }
670  }
671  }
672  }
673 
674  /* debug print population weighted transition probability */
675  {
676  enum {DEBUG_LOC=false};
677  if( DEBUG_LOC )
678  {
679  double sumpop = 0.;
680  double sumpopA = 0.;
681  for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
682  {
683  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
684  {
685  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
686  {
687  /* now the lower levels */
688  /* NB - X is the only lower level considered here, since we are only
689  * concerned with excited electronic levels as a photodissociation process
690  * code exists to relax this assumption - simply change following to iElecHi */
691  iElecLo = 0;
692  /* want to include all vibration states in lower level if different electronic level,
693  * but only lower vibration levels if same electronic level */
694  for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo )
695  {
696  long nr = h2.nRot_hi[iElecLo][iVibLo];
697  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
698  {
699  /* >>chng 03 feb 14, change test from !=0 to >0 */
700  /* >>chng 05 feb 07, use lgH2_line_exists */
701  /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )*/
702  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
703  {
704  ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul > 0. );
705 
706  sumpop += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop;
707  sumpopA += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop*
708  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul;
709  }
710  }
711  }
712  }
713  }
714  }
715  fprintf(ioQQQ,"DEBUG sumpop = %.3e sumpopA= %.3e A=%.3e\n", sumpop, sumpopA,
716  sumpopA/SDIV(sumpop) );
717  }
718  }
719  return;
720 }
721 
722 /* increment optical depth for the H2 molecule, called from RT_tau_inc which is called by cloudy,
723  * one time per zone */
724 void H2_RT_tau_inc(void)
725 {
726  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
727 
728  /* >>chng 05 jan 26, now use LTE populations for small H2 abundance case, since electronic
729  * lines become self-shielding surprisingly quickly */
730  if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ )
731  return;
732 
733  DEBUG_ENTRY( "H2_RT_tau_inc()" );
734 
735  /* remember largest and smallest chemistry renorm factor -
736  * if both networks are parallel will be unity,
737  * but only do this after we have stable solution */
738  if( nzone > 0 && nCallH2_this_iteration>2 )
739  {
742  }
743 
744  /* loop over all possible lines */
745  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
746  {
747  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
748  {
749  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
750  {
751  /* now the lower levels */
752  /* NB - X is the only lower level considered here, since we are only
753  * concerned with excited electronic levels as a photodissociation process
754  * code exists to relax this assumption - simply change following to iElecHi */
755  long int lim_elec_lo = 0;
756  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
757  {
758  /* want to include all vibration states in lower level if different electronic level,
759  * but only lower vibration levels if same electronic level */
760  long int nv = h2.nVib_hi[iElecLo];
761  if( iElecLo==iElecHi )
762  nv = iVibHi;
763  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
764  {
765  long nr = h2.nRot_hi[iElecLo][iVibLo];
766  if( iElecLo==iElecHi && iVibHi==iVibLo )
767  nr = iRotHi-1;
768 
769  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
770  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
771  {
772  /* >>chng 03 feb 14, from !=0 to >0 */
773  /* >>chng 05 feb 07, use lgH2_line_exists */
774  if( *lgH2le++ )
775  {
776  ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 );
777  RT_line_one_tauinc( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ,
778  -9 , iRotHi , iVibLo , iRotLo);
779  }
780  }/* iRotLo loop */
781  }/* iVibLo loop */
782  }/* iElecLo loop */
783  }/* iRotHi loop */
784  }/* iVibHi loop */
785  }/* iElecHi loop */
786  /*fprintf(ioQQQ,"\t%.3e\n",H2Lines[1][0][0][0][0][1].TauCon);*/
787  return;
788 }
789 
790 
791 /* initialize optical depths in H2, called from RT_tau_init */
792 void H2_LineZero( void )
793 {
794  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
795 
796  if( !h2.lgH2ON )
797  return;
798 
799  DEBUG_ENTRY( "H2_LineZero()" );
800 
801  /* loop over all possible lines */
802  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
803  {
804  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
805  {
806  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
807  {
808  /* now the lower levels */
809  /* NB - X is the only lower level considered here, since we are only
810  * concerned with excited electronic levels as a photodissociation process
811  * code exists to relax this assumption - simply change following to iElecHi */
812  long int lim_elec_lo = 0;
813  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
814  {
815  /* want to include all vibration states in lower level if different electronic level,
816  * but only lower vibration levels if same electronic level */
817  long int nv = h2.nVib_hi[iElecLo];
818  if( iElecLo==iElecHi )
819  nv = iVibHi;
820  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
821  {
822  long nr = h2.nRot_hi[iElecLo][iVibLo];
823  if( iElecLo==iElecHi && iVibHi==iVibLo )
824  nr = iRotHi-1;
825 
826  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
827  {
828  /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul != 0. )*/
829  /* >>chng 03 feb 14, from !=0 to > 0 */
830  /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )*/
831  /* >>chng 05 feb 07, use lgH2_line_exists */
832  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
833  {
834  /* >>chng 03 feb 14, use TransitionZero rather than explicit sets */
835  TransitionZero( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
836  }
837  }
838  }
839  }
840  }
841  }
842  }
843  return;
844 }
845 
846 /* the large H2 molecule, called from RT_tau_reset */
847 void H2_RT_tau_reset( void )
848 {
849  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
850 
851  if( !h2.lgH2ON )
852  return;
853 
854  DEBUG_ENTRY( "H2_RT_tau_reset()" );
855 
856  /* loop over all possible lines */
857  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
858  {
859  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
860  {
861  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
862  {
863  /* now the lower levels */
864  /* NB - X is the only lower level considered here, since we are only
865  * concerned with excited electronic levels as a photodissociation process
866  * code exists to relax this assumption - simply change following to iElecHi */
867  long int lim_elec_lo = 0;
868  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
869  {
870  /* want to include all vibration states in lower level if different electronic level,
871  * but only lower vibration levels if same electronic level */
872  long int nv = h2.nVib_hi[iElecLo];
873  if( iElecLo==iElecHi )
874  nv = iVibHi;
875  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
876  {
877  long nr = h2.nRot_hi[iElecLo][iVibLo];
878  if( iElecLo==iElecHi && iVibHi==iVibLo )
879  nr = iRotHi-1;
880 
881  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
882  {
883  /* >>chng 03 feb 14, change test from !=0 to >0 */
884  /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )*/
885  /* >>chng 05 feb 07, use lgH2_line_exists */
886  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
887  {
888  /* inward optical depth */
889  RT_line_one_tau_reset( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , 0.5);
890  }
891  }
892  }
893  }
894  }
895  }
896  }
897  return;
898 }
899 
900 /* this is fraction of population that is within levels done with matrix */
901 static double frac_matrix;
902 
903 /*H2_Level_low_matrix evaluate lower populations within X */
905  /* total abundance within matrix */
906  realnum abundance )
907 {
908 
909  /* will need to MALLOC space for these but only on first call */
910  static double **AulEscp ,
911  **col_str ,
912  **AulDest,
913  /* AulPump[low][high] is rate (s^-1) from lower to upper level */
914  **AulPump,
915  **CollRate_levn,
916  *pops,
917  *create,
918  *destroy,
919  *depart,
920  /* statistical weight */
921  *stat_levn ,
922  /* excitation energies in kelvin */
923  *excit;
924  bool lgDoAs;
925  static long int levelAsEval=-1;
926  long int ip;
927  static bool lgFirst=true;
928  long int i,
929  j,
930  ilo ,
931  ihi,
932  iElec,
933  iElecHi,
934  iVib,
935  iRot,
936  iVibHi,
937  iRotHi;
938  int nNegPop;
939  bool lgDeBug,
940  lgZeroPop;
941  double rot_cooling , dCoolDT;
942  static long int ndimMalloced = 0;
943  double rateout , ratein;
944 
945  DEBUG_ENTRY( "H2_Level_low_matrix()" );
946 
947  /* option to not use the matrix */
948  if( nXLevelsMatrix <= 1 )
949  {
950  return;
951  }
952 
953  if( lgFirst )
954  {
955  /* check that not more levels than there are in X */
957  {
958  /* number is greater than number of levels within X */
959  fprintf( ioQQQ,
960  " The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n",
961  nLevels_per_elec[0]);
962  cdEXIT(EXIT_FAILURE);
963  }
964  /* will never do this again */
965  lgFirst = false;
966  /* remember how much space we malloced in case ever called with more needed */
967  /* >>chng 05 jan 19, allocate max number of levels
968  ndimMalloced = nXLevelsMatrix;*/
969  ndimMalloced = nLevels_per_elec[0];
970  /* allocate the 1D arrays*/
971  excit = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
972  stat_levn = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
973  pops = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
974  create = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
975  destroy = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
976  depart = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) );
977  /* create space for the 2D arrays */
978  AulPump = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
979  CollRate_levn = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
980  AulDest = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
981  AulEscp = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
982  col_str = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
983  for( i=0; i<(ndimMalloced); ++i )
984  {
985  AulPump[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
986  CollRate_levn[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
987  AulDest[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
988  AulEscp[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
989  col_str[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
990  }
991 
992  for( j=0; j < ndimMalloced; j++ )
993  {
994  stat_levn[j]=0;
995  excit[j] =0;
996  }
997  /* the statistical weights of the levels
998  * and excitation potentials of each level relative to ground */
999  for( j=0; j < ndimMalloced; j++ )
1000  {
1001  /* obtain the proper indices for the upper level */
1002  ip = H2_ipX_ener_sort[j];
1003  iVib = ipVib_H2_energy_sort[ip];
1004  iRot = ipRot_H2_energy_sort[ip];
1005 
1006  /* statistical weights for each level */
1007  stat_levn[j] = H2_stat[0][iVib][iRot];
1008  /* excitation energy of each level relative to ground, in K */
1009  excit[j] = energy_wn[0][iVib][iRot]*T1CM;
1010  }
1011 
1012  for( j=0; j < ndimMalloced-1; j++ )
1013  {
1014  /* make sure that the energies are ok */
1015  ASSERT( excit[j+1] > excit[j] );
1016  }
1017  }
1018  /* end malloc space and creating constant terms */
1019 
1020  /* this is test for call with too many rotation levels to handle -
1021  * logic needs for largest model atom to be called first */
1022  if( nXLevelsMatrix > ndimMalloced )
1023  {
1024  fprintf(ioQQQ," H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n");
1025  cdEXIT(EXIT_FAILURE);
1026  }
1027 
1028  /* all elements are used, and must be set to zero */
1029  for( i=0; i < nXLevelsMatrix; i++ )
1030  {
1031  pops[i] = 0.;
1032  depart[i] = 0;
1033  for( j=0; j < nXLevelsMatrix; j++ )
1034  {
1035  col_str[j][i] = 0.;
1036  }
1037  }
1038 
1039  /* do we need to reevaluate radiative quantities? only do this one time per zone */
1040  if( nzone!=nzoneAsEval || iteration!=iterationAsEval || nXLevelsMatrix!=levelAsEval)
1041  {
1042  lgDoAs = true;
1043  nzoneAsEval = nzone;
1045  levelAsEval = nXLevelsMatrix;
1046  ASSERT( levelAsEval <= ndimMalloced );
1047  }
1048  else
1049  {
1050  lgDoAs = false;
1051  }
1052 
1053  /* all elements are used, and must be set to zero */
1054  if( lgDoAs )
1055  {
1056  for( i=0; i < nXLevelsMatrix; i++ )
1057  {
1058  pops[i] = 0.;
1059  depart[i] = 0;
1060  for( j=0; j < nXLevelsMatrix; j++ )
1061  {
1062  AulEscp[j][i] = 0.;
1063  AulDest[j][i] = 0.;
1064  AulPump[j][i] = 0.;
1065  CollRate_levn[j][i] = 0.;
1066  }
1067  }
1068  }
1069 
1070  /* find all radiative interactions within matrix, and between
1071  * matrix and upper X and excited electronic states */
1072  iElec = 0;
1073  for( ilo=0; ilo < nXLevelsMatrix; ilo++ )
1074  {
1075  ip = H2_ipX_ener_sort[ilo];
1076  iRot = ipRot_H2_energy_sort[ip];
1077  iVib = ipVib_H2_energy_sort[ip];
1078 
1079  /* H2_X_sink[ilo] includes all processes that destroy H2 in one step,
1080  * these include cosmic ray ionization and dissociation, photodissociation,
1081  * BUT NOT THE SOLOMON process, which, directly, only goes to excited
1082  * electronic states */
1083  destroy[ilo] = H2_X_sink[ilo];
1084 
1085  /* rates H2 is created from grains and H- units cm-3 s-1, evaluated in mole_H2_form */
1086  create[ilo] = H2_X_source[ilo];
1087 
1088  /* this loop does radiative decays from upper states inside matrix,
1089  * and upward pumps within matrix region into this lower level */
1090  if( lgDoAs )
1091  {
1092  for( ihi=ilo+1; ihi<nXLevelsMatrix; ++ihi )
1093  {
1094  ip = H2_ipX_ener_sort[ihi];
1095  iRotHi = ipRot_H2_energy_sort[ip];
1096  iVibHi = ipVib_H2_energy_sort[ip];
1097  ASSERT( H2_energies[ip] <= H2_energies[H2_ipX_ener_sort[nXLevelsMatrix-1]] );
1098  /* general case - but line may not actually exist */
1099  if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
1100  {
1101  /* >>chng 05 feb 07, use lgH2_line_exists */
1102  if( lgH2_line_exists[0][iVibHi][iRotHi][0][iVib][iRot] )
1103  {
1104  ASSERT( H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].ipCont > 0 );
1105 
1106  /* NB - the destruction probability is included in
1107  * the total and the destruction is set to zero
1108  * since we want to only count one ots rate, in
1109  * main calling routine, and do not want matrix
1110  * solver below to include it */
1111  AulEscp[ihi][ilo] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Aul*(
1112  H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pesc +
1113  H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pdest +
1114  H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pelec_esc);
1115  AulDest[ilo][ihi] = 0.;
1116  AulPump[ilo][ihi] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->pump;
1117  }
1118  }
1119  }
1120  }
1121 
1122  iElecHi = 0;
1123  iElec = 0;
1124  rateout = 0.;
1125  ratein = 0.;
1126  /* now do all levels within X, which are above nXLevelsMatrix,
1127  * the highest level inside the matrix */
1128  for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
1129  {
1130  ip = H2_ipX_ener_sort[ihi];
1131  iRotHi = ipRot_H2_energy_sort[ip];
1132  iVibHi = ipVib_H2_energy_sort[ip];
1133  if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
1134  {
1135  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot] )
1136  {
1137  ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].ipCont > 0 );
1138 
1139  /* these will enter as net creation terms in creation vector, with
1140  * units cm-3 s-1
1141  * radiative transitions from above the matrix within X */
1142  ratein +=
1143  H2_populations[iElecHi][iVibHi][iRotHi] *
1144  (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Aul*
1145  (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pesc +
1146  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pelec_esc +
1147  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pdest)+H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->pump *
1148  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Lo->g/
1149  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Hi->g);
1150  /* rate out has units s-1 - destroys current lower level */
1151  rateout +=
1152  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->pump;
1153  }
1154  }
1155  }
1156 
1157  /* all states above the matrix but within X */
1158  create[ilo] += ratein;
1159 
1160  /* rates out of matrix into levels in X but above matrix */
1161  destroy[ilo] += rateout;
1162 
1163  /* Solomon process, this sum dos all pump and decays from all electronic excited states */
1164  /* radiative rates [cm-3 s-1] from electronic excited states into X only vibration and rot */
1165  create[ilo] += H2_X_rate_from_elec_excited[iVib][iRot];
1166 
1167  /* radiative & cosmic ray rates [s-1] to electronic excited states from X only vibration and rot */
1168  destroy[ilo] += H2_X_rate_to_elec_excited[iVib][iRot];
1169  }
1170 
1171  /* this flag set with atom H2 trace matrix */
1173  lgDeBug = true;
1174  else
1175  lgDeBug = false;
1176 
1177  /* now evaluate the rates for all transitions within matrix */
1178  for( ilo=0; ilo < nXLevelsMatrix; ilo++ )
1179  {
1180  ip = H2_ipX_ener_sort[ilo];
1181  iRot = ipRot_H2_energy_sort[ip];
1182  iVib = ipVib_H2_energy_sort[ip];
1183  if( lgDoAs )
1184  {
1185  if(lgDeBug)fprintf(ioQQQ,"DEBUG H2_Level_low_matrix, ilo=%li",ilo);
1186  for( ihi=ilo+1; ihi < nXLevelsMatrix; ihi++ )
1187  {
1188  ip = H2_ipX_ener_sort[ihi];
1189  iRotHi = ipRot_H2_energy_sort[ip];
1190  iVibHi = ipVib_H2_energy_sort[ip];
1191  /* >>chng 05 may 31, replace with simple expresion */
1192  CollRate_levn[ihi][ilo] = H2_X_coll_rate[ihi][ilo];
1193 
1194  /*create[ilo] +=CollRate_levn[ihi][ilo]*H2_populations[0][iVibHi][iRotHi];*/
1195  if(lgDeBug)fprintf(ioQQQ,"\t%.1e",CollRate_levn[ihi][ilo]);
1196 
1197  /* now get upward excitation rate - units s-1 */
1198  CollRate_levn[ilo][ihi] = CollRate_levn[ihi][ilo]*
1199  H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
1200  H2_stat[0][iVibHi][iRotHi] /
1201  H2_stat[0][iVib][iRot];
1202  }
1203  }
1204 
1205  if(lgDeBug)fprintf(ioQQQ,"\n");
1206 
1207  /* now do all collisions for levels within X, which are above nXLevelsMatrix,
1208  * the highest level inside the matrix */
1209  iElecHi = 0;
1210 
1211  for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
1212  {
1213  ip = H2_ipX_ener_sort[ihi];
1214  iRotHi = ipRot_H2_energy_sort[ip];
1215  iVibHi = ipVib_H2_energy_sort[ip];
1216  rateout = 0;
1217  /* first do downward deexcitation rate */
1218  /* >>chng 04 sep 14, do all levels */
1219  /* >>chng 05 may 31, use summed rate */
1220  ratein = H2_X_coll_rate[ihi][ilo];
1221  if(lgDeBug)fprintf(ioQQQ,"\t%.1e",ratein);
1222 
1223  /* now get upward excitation rate */
1224  rateout = ratein *
1225  H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
1226  H2_stat[0][iVibHi][iRotHi]/H2_stat[0][iVib][iRot];
1227 
1228  /* these are general entries and exits going into vector */
1229  create[ilo] += ratein*H2_populations[iElecHi][iVibHi][iRotHi];
1230  destroy[ilo] += rateout;
1231  }
1232  }
1233 
1234  /* H2 grain interactions
1235  * >>chng 05 apr 30,GS, Instead of hmi.H2_total, the specific populations are used because high levels have much less
1236  * populations than ground levels which consists most of the H2 population.*/
1237  if( lgDoAs )
1238  {
1239  for( ihi=2; ihi < nXLevelsMatrix; ihi++ )
1240  {
1241 
1242  ip = H2_ipX_ener_sort[ihi];
1243  iVibHi = ipVib_H2_energy_sort[ip];
1244  iRotHi = ipRot_H2_energy_sort[ip];
1245 
1246  /* collisions with grains goes to either J=1 or J=0 depending on
1247  * spin of upper level - this conserves op ratio - following
1248  * var is 1 if ortho, 0 if para, so this conserves op ratio
1249  * units are s-1 */
1250  CollRate_levn[ihi][H2_lgOrtho[0][iVibHi][iRotHi]] += hmi.rate_grain_h2_op_conserve;
1251  }
1252 
1253  /* H2 ortho - para conversion on grain surface,
1254  * rate (s-1) all v,J levels go to 0 or 1 */
1255  CollRate_levn[1][0] +=
1257  }
1258 
1259  /* now all levels in X above the matrix */
1260  for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
1261  {
1262  ip = H2_ipX_ener_sort[ihi];
1263  iVibHi = ipVib_H2_energy_sort[ip];
1264  iRotHi = ipRot_H2_energy_sort[ip];
1265 
1266  /* these collisions all go into 0 or 1 depending on whether upper level was ortho or para
1267  * units are cm-3 s-1 - rate new molecules appear in matrix */
1268  create[H2_lgOrtho[0][iVibHi][iRotHi]] += H2_populations[0][iVibHi][iRotHi]*hmi.rate_grain_h2_op_conserve;
1269  }
1270 
1271  /* debug print individual contributors to matrix elements */
1272  {
1273  enum {DEBUG_LOC=false};
1274  if( DEBUG_LOC || lgDeBug)
1275  {
1276  fprintf(ioQQQ,"DEBUG H2 matexcit");
1277  for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
1278  {
1279  fprintf(ioQQQ,"\t%li",ilo );
1280  }
1281  fprintf(ioQQQ,"\n");
1282  for(ihi=0; ihi<nXLevelsMatrix;++ihi)
1283  {
1284  fprintf(ioQQQ,"\t%.2e",excit[ihi] );
1285  }
1286  fprintf(ioQQQ,"\n");
1287  for(ihi=0; ihi<nXLevelsMatrix;++ihi)
1288  {
1289  fprintf(ioQQQ,"\t%.2e",stat_levn[ihi] );
1290  }
1291  fprintf(ioQQQ,"\n");
1292 
1293  fprintf(ioQQQ,"AulEscp[n][]\\[][n] = Aul*Pesc\n");
1294  for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
1295  {
1296  fprintf(ioQQQ,"\t%li",ilo );
1297  }
1298  fprintf(ioQQQ,"\n");
1299  for(ihi=0; ihi<nXLevelsMatrix;++ihi)
1300  {
1301  fprintf(ioQQQ,"%li", ihi);
1302  for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
1303  {
1304  fprintf(ioQQQ,"\t%.2e",AulEscp[ilo][ihi] );
1305  }
1306  fprintf(ioQQQ,"\n");
1307  }
1308 
1309  fprintf(ioQQQ,"AulPump [n][]\\[][n]\n");
1310  for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
1311  {
1312  fprintf(ioQQQ,"\t%li",ilo );
1313  }
1314  fprintf(ioQQQ,"\n");
1315  for(ihi=0; ihi<nXLevelsMatrix;++ihi)
1316  {
1317  fprintf(ioQQQ,"%li", ihi);
1318  for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
1319  {
1320  fprintf(ioQQQ,"\t%.2e",AulPump[ihi][ilo] );
1321  }
1322  fprintf(ioQQQ,"\n");
1323  }
1324 
1325  fprintf(ioQQQ,"CollRate_levn [n][]\\[][n]\n");
1326  for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
1327  {
1328  fprintf(ioQQQ,"\t%li",ilo );
1329  }
1330  fprintf(ioQQQ,"\n");
1331  for(ihi=0; ihi<nXLevelsMatrix;++ihi)
1332  {
1333  fprintf(ioQQQ,"%li", ihi);
1334  for(ilo=0; ilo<nXLevelsMatrix; ++ilo )
1335  {
1336  fprintf(ioQQQ,"\t%.2e",CollRate_levn[ihi][ilo] );
1337  }
1338  fprintf(ioQQQ,"\n");
1339  }
1340  fprintf(ioQQQ,"SOURCE");
1341  for(ihi=0; ihi<nXLevelsMatrix;++ihi)
1342  {
1343  fprintf(ioQQQ,"\t%.2e",create[ihi]);
1344  }
1345  fprintf(ioQQQ,"\nSINK");
1346  for(ihi=0; ihi<nXLevelsMatrix;++ihi)
1347  {
1348  fprintf(ioQQQ,"\t%.2e",destroy[ihi]);
1349  }
1350  fprintf(ioQQQ,"\n");
1351  }
1352  }
1353 
1354  atom_levelN(
1355  /* number of levels */
1356  nXLevelsMatrix,
1357  abundance,
1358  stat_levn,
1359  excit,
1360  'K',
1361  pops,
1362  depart,
1363  /* net transition rate, A * escape prob, s-1, indices are [upper][lower] */
1364  &AulEscp,
1365  /* col str from high to low */
1366  &col_str,
1367  &AulDest,
1368  &AulPump,
1369  &CollRate_levn,
1370  create,
1371  destroy,
1372  /* say that we have evaluated the collision rates already */
1373  true,
1374  &rot_cooling,
1375  &dCoolDT,
1376  " H2 ",
1377  /* nNegPop positive if negative pops occurred, negative if too cold */
1378  &nNegPop,
1379  &lgZeroPop,
1380  lgDeBug );/* option to print stuff - set to true for debug printout */
1381 
1382  for( i=0; i< nXLevelsMatrix; ++i )
1383  {
1384  ip = H2_ipX_ener_sort[i];
1385  iRot = ipRot_H2_energy_sort[ip];
1386  iVib = ipVib_H2_energy_sort[ip];
1387  /* >>chng 05 feb 08, do not update h2_old_populations here, since not done
1388  * like this anywhere else - only update H2_populations here, and let single loop
1389  * in main calling routine handle updating various forms of the population
1390  H2_old_populations[0][iVib][iRot] = H2_populations[0][iVib][iRot]; */
1391  H2_populations[0][iVib][iRot] = pops[i];
1392  }
1393 
1394  if( 0 && mole.nH2_TRACE >= mole.nH2_trace_full)
1395  {
1396  /*static int nn=0; ++nn; if( nn>5)cdEXIT(1);*/
1397  /* print pops that came out of matrix */
1398  fprintf(ioQQQ,"\n DEBUG H2_Level_lowJ hmi.H2_total: %.3e matrix rel pops\n",hmi.H2_total);
1399  fprintf(ioQQQ,"v\tJ\tpop\n");
1400  for( i=0; i<nXLevelsMatrix; ++i )
1401  {
1402  ip = H2_ipX_ener_sort[i];
1403  iRot = ipRot_H2_energy_sort[ip];
1404  iVib = ipVib_H2_energy_sort[ip];
1405  fprintf(ioQQQ,"%3li\t%3li\t%.3e\t%.3e\t%.3e\n",
1406  iVib , iRot , H2_populations[0][iVib][iRot]/hmi.H2_total , create[i] , destroy[i]);
1407  }
1408  }
1409 
1410  /* nNegPop positive if negative pops occurred, negative if too cold */
1411  if( nNegPop > 0 )
1412  {
1413  fprintf(ioQQQ," H2_Level_low_matrix called atom_levelN which returned negative H2_populations.\n");
1414  ConvFail( "pops" , "H2" );
1415  }
1416  return;
1417 }
1418 /* do level H2_populations for H2, called by Hydrogenic after ionization and H chemistry
1419  * has been recomputed */
1420 void H2_LevelPops( void )
1421 {
1422  static double TeUsedColl=-1.f;
1423  double H2_renorm_conserve=0.,
1424  H2_renorm_conserve_init=0. ,
1425  sumold,
1426  H2_BigH2_H2s,
1427  H2_BigH2_H2g;
1428  double old_solomon_rate=-1.;
1429  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
1430  long int i;
1431  long int n_pop_oscil = 0;
1432  int kase=0;
1433  bool lgConv_h2_soln,
1434  lgPopsConv_total,
1435  lgPopsConv_relative,
1436  lgHeatConv,
1437  lgSolomonConv,
1438  lgOrthoParaRatioConv;
1439  double quant_old=-1.,
1440  quant_new=-1.;
1441 
1442  bool lgH2_pops_oscil=false,
1443  lgH2_pops_ever_oscil=false;
1444  long int nEner,
1445  ipHi, ipLo;
1446  long int iElec , iVib , iRot,ip;
1447  double sum_pops_matrix;
1448  realnum collup;
1449  /* old and older ortho - para ratios, used to determine whether soln is converged */
1450  static double ortho_para_old=0. , ortho_para_older=0. , ortho_para_current=0.;
1451  realnum frac_new_oscil=1.f;
1452 
1453  /* keep track of changes in population */
1454  double PopChgMax_relative=0. , PopChgMaxOld_relative=0., PopChgMax_total=0., PopChgMaxOld_total=0.;
1455  long int iRotMaxChng_relative , iVibMaxChng_relative,
1456  iRotMaxChng_total , iVibMaxChng_total,
1457  nXLevelsMatrix_save;
1458  double popold_relative , popnew_relative , popold_total , popnew_total;
1459  /* reason not converged */
1460  char chReason[100];
1461 
1462  double flux_accum_photodissoc_BigH2_H2g, flux_accum_photodissoc_BigH2_H2s;
1463  long int ip_H2_level;
1464 
1465  /* these are convergence criteria - will be increased during search phase */
1466  double converge_pops_relative=0.1 ,
1467  converge_pops_total=1e-3,
1468  converge_ortho_para=1e-2;
1469 
1470  /* H2 not on, so space not allocated and return,
1471  * also return if calculation has been declared a failure */
1472  if( !h2.lgH2ON || lgAbort )
1473  return;
1474 
1475  DEBUG_ENTRY( "H2_LevelPops()" );
1476 
1478  {
1479  fprintf(ioQQQ,
1480  "\n***************H2_LevelPops call %li this iteration, zone is %.2f, H2/H:%.e Te:%e ne:%e\n",
1482  fnzone,
1484  phycon.te,
1485  dense.eden
1486  );
1487  }
1488  else if( mole.nH2_TRACE >= mole.nH2_trace_final )
1489  {
1490  static long int nzone_prt=-1;
1491  if( nzone!=nzone_prt )
1492  {
1493  nzone_prt = nzone;
1494  fprintf(ioQQQ,"DEBUG zone %li H2/H:%.3e Te:%.3e *ne:%.3e n(H2):%.3e\n",
1495  nzone,
1497  phycon.te,
1498  dense.eden,
1499  hmi.H2_total );
1500  }
1501  }
1502 
1503  /* evaluate Boltzmann factors and LTE unit population - for trivial abundances
1504  * LTE populations are used in place of full solution */
1505  mole_H2_LTE();
1506 
1507  /* zero out H2_populations and cooling, and return, if H2 fraction is small
1508  * but, if H2 has ever been done, redo irregardless of abundance -
1509  * if large H2 is ever evaluated then mole.H2_to_H_limit is ignored */
1511  || hmi.H2_total < 1e-20 )
1512  {
1513  /* will not compute the molecule */
1514  if( mole.nH2_TRACE >= mole.nH2_trace_full )
1515  fprintf(ioQQQ,
1516  " H2_LevelPops pops too small, not computing, set to LTE and return, H2/H is %.2e and mole.H2_to_H_limit is %.2e.",
1520  /* end of zero abundance branch */
1521  return;
1522  }
1523 
1524  /* check whether we need to update the H2_Boltzmann factors, LTE level H2_populations,
1525  * and partition function. LTE level pops normalized by partition function,
1526  * so sum of pops is unity */
1527 
1528  /* say that H2 has been computed, ignore previous limit to abund
1529  * in future - this is to prevent oscillations as model is engaged */
1530  hmi.lgBigH2_evaluated = true;
1531  /* end loop setting H2_Boltzmann factors, partition function, and LTE H2_populations */
1532 
1533  /* >>chng 05 jun 21,
1534  * during search phase we want to use full matrix - save number of levels so that
1535  * we can restore it */
1536  nXLevelsMatrix_save = nXLevelsMatrix;
1537  if( conv.lgSearch )
1538  {
1540  }
1541 
1542  /* 05 oct 27, had only reevaluated collision rates when 5% change in temperature
1543  * caused temp failures in large G0 sims -
1544  * do not check whether we need to update the collision rates but
1545  * reevaluate them always
1546  * >>chng 05 nov 04, above caused a 25% increase in the exec time for constant-T sims
1547  * in test suite- original code had reevaluated if > 0.05 change in T - was too much
1548  * change to 10x smaller, change > 0.005 */
1549  if( fabs(1. - TeUsedColl / phycon.te ) > 0.005 )
1550  {
1552  TeUsedColl = phycon.te;
1553  }
1554 
1555  /* set the H2_populations when this is the first call to this routine on
1556  * current iteration- will use LTE H2_populations - populations were set by
1557  * call to mole_H2_LTE before above block */
1559  {
1560  /* very first call so use LTE H2_populations */
1562  fprintf(ioQQQ,"H2 1st call - using LTE level pops\n");
1563 
1564  for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
1565  {
1566  pops_per_elec[iElec] = 0.;
1567  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
1568  {
1569  pops_per_vib[iElec][iVib] = 0.;
1570  for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
1571  {
1572  /* LTE populations are for unit H2 density, so need to multiply
1573  * by total H2 density */
1574  H2_old_populations[iElec][iVib][iRot] =
1575  (realnum)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total;
1576  H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot];
1577  }
1578  }
1579  }
1580  /* first guess at ortho and para densities */
1581  h2.ortho_density = 0.75*hmi.H2_total;
1582  h2.para_density = 0.25*hmi.H2_total;
1583  ortho_para_current = h2.ortho_density / SDIV( h2.para_density );
1584  ortho_para_older = ortho_para_current;
1585  ortho_para_old = ortho_para_current;
1586  /* this is the fraction of the H2 pops that are within the levels done with a matrix */
1587  frac_matrix = 1.;
1588  }
1589 
1590  /* find sum of all H2_populations in X */
1591  iElec = 0;
1592  pops_per_elec[0] = 0.;
1593  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
1594  {
1595  pops_per_vib[0][iVib] = 0.;
1596  for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
1597  {
1598  pops_per_elec[0] += H2_populations[iElec][iVib][iRot];
1599  pops_per_vib[0][iVib] += H2_populations[iElec][iVib][iRot];
1600  }
1601  }
1603  /* now renorm the old populations to the correct current H2 density,
1604  * At this point pops_per_elec[0] (pop of X)
1605  * is the result of the previous evaluation of the H2 population,
1606  * following is the ratio of the current chemistry solution H2 to the previous H2 */
1608 
1609  /* >>chng 05 jul 13, TE,
1610  * evaluate ratio of H2g and H2s from chemical network and big molecule model */
1611  iElec = 0;
1612  hmi.H2_chem_BigH2_H2g = 0.;
1613  hmi.H2_chem_BigH2_H2s = 0.;
1614  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
1615  {
1616  for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
1617  {
1618  if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
1619  {
1620  hmi.H2_chem_BigH2_H2s += H2_populations[iElec][iVib][iRot];
1621 
1622  }
1623  else
1624  {
1625  hmi.H2_chem_BigH2_H2g += H2_populations[iElec][iVib][iRot];
1626 
1627  }
1628  }
1629  }
1630 
1633 
1634 
1636  fprintf(ioQQQ,
1637  "H2 H2_renorm_chemistry is %.4e, hmi.H2_total is %.4e pops_per_elec[0] is %.4e\n",
1639  hmi.H2_total,
1640  pops_per_elec[0]);
1641 
1642  /* renormalize all level populations for the current chemical solution */
1643  iElec = 0;
1644  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
1645  {
1646  for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
1647  {
1648  H2_populations[iElec][iVib][iRot] *= H2_renorm_chemistry;
1649  H2_old_populations[iElec][iVib][iRot] = H2_populations[iElec][iVib][iRot];
1650  }
1651  }
1653 
1655  fprintf(ioQQQ,
1656  " H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n",
1657  pops_per_elec[0],
1658  hmi.H2_total);
1659 
1660  /* >>chng 05 feb 10, reset convergence criteria if we are in search phase */
1661  if( conv.lgSearch )
1662  {
1663  converge_pops_relative *= 2.; /*def is 0.1 */
1664  converge_pops_total *= 3.; /*def is 1e-3*/
1665  converge_ortho_para *= 3.; /*def is 1e-2*/
1666  }
1667 
1668  /* update state specific rates in H2_X_formation (cm-3 s-1) that H2 forms from grains and H- */
1669  mole_H2_form();
1670 
1671  /* evaluate total collision rates */
1673 
1674  /* this flag will say whether H2 H2_populations have converged,
1675  * by comparing old and new values */
1676  lgConv_h2_soln = false;
1677  /* this will count number of passes around following loop */
1678  loop_h2_pops = 0;
1679  {
1680  static long int nzoneEval=-1;
1681  if( nzone != nzoneEval )
1682  {
1683  nzoneEval = nzone;
1684  /* this is number of zones with H2 solution in this iteration */
1685  ++nH2_zone;
1686  }
1687  }
1688 
1689  /* begin - start level population solution
1690  * first do electronic excited states, Lyman, Werner, etc
1691  * using old solution for X
1692  * then do matrix if used, then solve for pops of rest of X
1693  * >>chng 04 apr 06, subtract number of oscillations from limit - don't waste loops
1694  * if solution is unstable */
1695  while( loop_h2_pops < LIM_H2_POP_LOOP-n_pop_oscil && !lgConv_h2_soln && !mole.lgH2_LTE )
1696  {
1697  double rate_in , rate_out;
1698  static double old_HeatH2Dexc_BigH2=0., HeatChangeOld=0. , HeatChange=0.;
1699 
1700  /* this is number of trips around loop this time */
1701  ++loop_h2_pops;
1702  /* this is number of times through this loop in entire iteration */
1703  ++nH2_pops;
1704 
1705  /* beginning solution for electronic excited states
1706  * loop over all possible pumping routes to excited electronic states
1707  * to get radiative excitation and dissociation rates */
1709 
1710  rate_out = 0.;
1711 
1712  /* these will store radiative rates between electronic excited states and X */
1713  iElec = 0;
1714  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
1715  {
1716  for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
1717  {
1718  /* radiative rates [cm-3 s-1] from electronic excited states into X vibration and rot */
1719  H2_X_rate_from_elec_excited[iVib][iRot] = 0.;
1720  /* radiative & cosmic ray rates [s-1] to electronic excited states from X */
1721  H2_X_rate_to_elec_excited[iVib][iRot] = 0.;
1722  }
1723  }
1724 
1725  iElecHi = -INT32_MAX;
1726  /* solve for population of electronic excited states by back-substitution */
1727  for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
1728  {
1729  /* for safety, these loop vars are set to insane value, will crash
1730  * if used without being set */
1731  iVib = -INT32_MAX;
1732  iRot = -INT32_MAX;
1733  iVibLo = -INT32_MAX;
1734  iRotLo = -INT32_MAX;
1735  iVibHi = -INT32_MAX;
1736  iRotHi = -INT32_MAX;
1737 
1738  /* this will be total population in each electronic state */
1739  pops_per_elec[iElecHi] = 0.;
1740 
1742  fprintf(ioQQQ," Pop(e=%li):",iElecHi);
1743 
1744  double factor = ( hmi.lgLeidenCRHack ) ? secondaries.x12tot/6.e-17 : 0.;
1745 
1746  /* electronic excited state, loop over all vibration */
1747  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
1748  {
1749  pops_per_vib[iElecHi][iVibHi] = 0.;
1750  /* ======================= EXCITED ELEC STATE POPS LOOP =====================*/
1751  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
1752  {
1753  /* Solomon process done here,
1754  * sum of all rates into and out of these upper levels
1755  * all inward rates have units cm-3 s-1 */
1756  rate_in = 0.;
1757  /* this term is spontaneous dissociation of excited electronic states into
1758  * the X continuum
1759  * all outward rates have units s-1 */
1760  rate_out = H2_dissprob[iElecHi][iVibHi][iRotHi];
1761 
1762  realnum H2gHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->g;
1763 
1764  /* now loop over all levels within X find Solomon rate */
1765  iElecLo=0;
1766 
1767  for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo )
1768  {
1769  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
1770  mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
1771  long nr = h2.nRot_hi[iElecLo][iVibLo];
1772  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1773  {
1774  /* consider all radiatively permitted transitions between X and excit states */
1775  if( *lgH2le++ )
1776  {
1777  ASSERT( H2L->ipCont > 0 );
1778 
1779  /*>>chng 07 jan 4, GS, add collisional excitation by non-thermal electrons in to
1780  * the Singlet state of H2 (B,C,B',D) from equation 5 of Liu & Dalgarno 1994 ApJ, 428, 769,
1781  * assumed incoming energy and outgoing energies are 20eV and 10eV;Eqn 5 gives downward cross
1782  * section and I multiply it by statistical weights to convert into upward cross section */
1783  double cr_cross_section = H2L->Coll.cs;
1784 # if 0
1785  /* one-time evaluation of this done in mole_h2_create.cpp */
1786  cr_cross_section =
1787  pow(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng*1e-8,3)*
1788  (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g/
1789  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g)*
1790  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*
1791  log(3.)*HPLANCK/(160.f*pow(PI,3)*0.5*1e-8*1.6e-12);
1792  ASSERT( fabs(cr_cross_section-H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cs)/
1793  SDIV(cr_cross_section) < 1e-7 );
1794 # endif
1795 
1796  /* solve electronic excited state,
1797  * rate lower level in X goes to electronic excited state, s-1
1798  * first term is direct pump, second is cosmic ray excitation */
1799  double rate_one = H2L->Emis->pump
1800  /*>>chng 05 jun 16, GS, add collisional excitation by non-thermal electrons in to
1801  * the Singlet state of H2 (B,C,B',D) from Dalgarno, Yan, & Liu 1999 ApJs;*/
1802  /*>>chng 07 jan 4, GS, add collisional excitation by non-thermal electrons in to
1803  * the Singlet state of H2 (B,C,B',D) from equation 5 of
1804  *>>refer H2 cr exc Liu & Dalgarno 1994 ApJ, 428, 769
1805  * in terms of hydrogen ionization cross-section
1806  >>refer H2 elec coll Shemansky, D.E., Hall, D.T.& Ajello, J.M. 1985ApJ, 296, 765 */
1807  +factor*cr_cross_section;
1808 
1809  /* this is a permitted electronic transition, must preserve nuclear spin */
1810  ASSERT( H2_lgOrtho[iElecHi][iVibHi][iRotHi] == H2_lgOrtho[iElecLo][iVibLo][iRotLo] );
1811 
1812  /* this is the rate [cm-3 s-1] electrons move into the upper level from X */
1813  rate_in += H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_one;
1814 
1815  /* this is total X -> excited electronic state rate, cm-3 s-1 */
1816  /*rate_tot += H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_one;*/
1817 
1818  /* rate [s-1] from levels within X to electronic excited states,
1819  * includes photoexcitation and cosmic ray excitation */
1820  H2_X_rate_to_elec_excited[iVibLo][iRotLo] += rate_one;
1821 
1822  /* excitation rate for Solomon process - this currently has units
1823  * cm-3 s-1 but will be divided by total pop of X and become s-1 */
1824  /* >>chng 04 may 25, Gargi Shaw found this bug - had been total sum */
1825  /* hmi.H2_H2g_to_H2s_rate_BigH2 += rate_one/SDIV(hmi.H2_total);*/
1826  /* this has unit s-1 and will be used for H2g->H2s */
1827  rate_one = H2L->Emis->Aul*
1828  /* escape and destruction */
1829  (H2L->Emis->Pesc +
1830  H2L->Emis->Pelec_esc +
1831  H2L->Emis->Pdest) +
1832  /* induced emission down */
1833  H2L->Emis->pump *
1834  H2L->Lo->g/H2gHi;
1835 
1836  /* this is the rate [s-1] electrons leave the excited electronic upper level
1837  * and decay into X - will be used to get pops of electronic excited states */
1838  rate_out += rate_one;
1839 
1840  ASSERT( rate_in >= 0. && rate_out >= 0. );
1841  }
1842  ++H2L;
1843  }
1844  }
1845 
1846  /* update population [cm-3] of the electronic excited state this only includes
1847  * radiative processes between X and excited electronic states, and cosmic rays -
1848  * thermal collisions are neglected
1849  * X is done below and includes all processes */
1850  H2_rad_rate_out[iElecHi][iVibHi][iRotHi] = rate_out;
1851  double H2popHi = rate_in / SDIV( rate_out );
1852  H2_populations[iElecHi][iVibHi][iRotHi] = H2popHi;
1853  if( H2_old_populations[iElecHi][iVibHi][iRotHi]==0. )
1854  H2_old_populations[iElecHi][iVibHi][iRotHi] = H2popHi;
1855 
1856  /* for this population, find rate electronic states decay into X */
1857  /* now loop over all levels within X find Solomon rate */
1858  iElecLo=0;
1859  for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo )
1860  {
1861  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
1862  mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
1863 
1864  long nr = h2.nRot_hi[iElecLo][iVibLo];
1865  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1866  {
1867  if( *lgH2le++ )
1868  {
1869  ASSERT( H2L->ipCont > 0 );
1870 
1871  double rate_one =
1872  /* >>chng 05 may 31, from old to new as in matrix */
1873  /*H2_old_populations[iElecHi][iVibHi][iRotHi]* */
1874  H2popHi*
1875  (H2L->Emis->Aul*
1876  /* escape and destruction */
1877  (H2L->Emis->Pesc + H2L->Emis->Pelec_esc + H2L->Emis->Pdest) +
1878  /* induced emission down */
1879  H2L->Emis->pump * H2L->Lo->g/H2gHi);
1880 
1881  /* radiative rates [cm-3 s-1] from electronic excited states to X */
1882  H2_X_rate_from_elec_excited[iVibLo][iRotLo] += rate_one;
1883  }
1884  ++H2L;
1885  }
1886  }
1887 
1888  ASSERT( H2popHi >= 0. && H2popHi <= hmi.H2_total );
1889 
1890  /* this is total pop in this vibration state */
1891  pops_per_vib[iElecHi][iVibHi] += H2popHi;
1892 
1893  /* ======================= POPS EXCITED ELEC CONVERGE LOOP =====================*/
1894  }/* end excit electronic state rot pops loop */
1895 
1897  fprintf(ioQQQ,"\t%.2e",pops_per_vib[iElecHi][iVibHi]/hmi.H2_total);
1898 
1899  /* total pop in each electronic state */
1900  pops_per_elec[iElecHi] += pops_per_vib[iElecHi][iVibHi];
1901  }/* end excit electronic state all vibration loop */
1902 
1903  /* end excited electronic pops loop */
1905  fprintf(ioQQQ,"\n");
1906  } /* end loop over all electronic excited states */
1907  /*fprintf(ioQQQ,"DEBUG\t%.3e\t%.3e\n",
1908  H2_X_rate_from_elec_excited[0][0],
1909  pops_per_elec[0] );*/
1910  /* above set pops of excited electronic levels and found rates between them and X -
1911  * now solve highly excited levels within the X state by back-substitution */
1912  /* these will do convergence check */
1913  PopChgMaxOld_relative = PopChgMax_relative;
1914  PopChgMaxOld_total = PopChgMax_total;
1915  PopChgMax_relative = 0.;
1916  PopChgMax_total = 0.;
1917  iElec = 0;
1918  iElecHi = 0;
1919  iRotMaxChng_relative =-1;
1920  iVibMaxChng_relative = -1;
1921  iRotMaxChng_total =-1;
1922  iVibMaxChng_total = -1;
1923  popold_relative = 0.;
1924  popnew_relative = 0.;
1925  popold_total = 0.;
1926  popnew_total = 0.;
1927 
1928  /* now evaluate total rates for all levels within X */
1929  for( nEner=0; nEner<nLevels_per_elec[0]; ++nEner )
1930  {
1931  /* array of energy sorted indices within X */
1932  ip = H2_ipX_ener_sort[nEner];
1933  iVib = ipVib_H2_energy_sort[ip];
1934  iRot = ipRot_H2_energy_sort[ip];
1935 
1936  realnum H2stat = H2_stat[0][iVib][iRot];
1937  double H2boltz = H2_Boltzmann[0][iVib][iRot];
1938 
1939  /* these will be total rates into and out of the level */
1940  double col_rate_in = 0.;
1941  double col_rate_out = 0.;
1942  for( ipLo=0; ipLo<nEner; ++ipLo )
1943  {
1944  ip = H2_ipX_ener_sort[ipLo];
1945  iVibLo = ipVib_H2_energy_sort[ip];
1946  iRotLo = ipRot_H2_energy_sort[ip];
1947 
1948  /* this is rate from this level down to lower level, units s-1 */
1949  col_rate_out += H2_X_coll_rate[nEner][ipLo];
1950  /*if( nEner==288 ) fprintf(ioQQQ,"DEBUG %3li %3li %.3e\n", nEner,ipLo,H2_X_coll_rate[nEner][ipLo]);*/
1951 
1952  /* inverse, rate up, cm-3 s-1 */
1953  collup = (realnum)(H2_old_populations[0][iVibLo][iRotLo] * H2_X_coll_rate[nEner][ipLo] *
1954  H2stat / H2_stat[0][iVibLo][iRotLo] *
1955  H2boltz / SDIV( H2_Boltzmann[0][iVibLo][iRotLo] ) );
1956 
1957  col_rate_in += collup;
1958  }
1959 
1960  for( ipHi=nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi )
1961  {
1962  double colldn;
1963  ip = H2_ipX_ener_sort[ipHi];
1964  iVibHi = ipVib_H2_energy_sort[ip];
1965  iRotHi = ipRot_H2_energy_sort[ip];
1966 
1967  /* this is rate from this level up to higher level, units s-1 */
1968  col_rate_out += H2_X_coll_rate[ipHi][nEner] *
1969  H2_stat[0][iVibHi][iRotHi] / H2stat *
1970  (realnum)(H2_Boltzmann[0][iVibHi][iRotHi] / SDIV( H2boltz ) );
1971 
1972  /* rate down from higher level, cm-3 s-1 */
1973  colldn = H2_old_populations[0][iVibHi][iRotHi] * H2_X_coll_rate[ipHi][nEner];
1974  /*if( nEner==288 ) fprintf(ioQQQ,"DEBUG %3li %3li %.3e\n", nEner,ipHi,H2_X_coll_rate[ipHi][nEner] *
1975  H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVib][iRot] *
1976  (realnum)(H2_Boltzmann[0][iVibHi][iRotHi] /
1977  SDIV( H2_Boltzmann[0][iVib][iRot] ) ));*/
1978 
1979  col_rate_in += colldn;
1980  }
1981 
1982  H2_col_rate_in[iVib][iRot] = col_rate_in;
1983  H2_col_rate_out[iVib][iRot] = col_rate_out;
1984  }
1985 
1986  /* =======================INSIDE X POPULATIONS CONVERGE LOOP =====================*/
1987  /* begin solving for X by back-substitution
1988  * this is the main loop that determines H2_populations within X
1989  * units of all rates in are cm-3 s-1, all rates out are s-1
1990  * nLevels_per_elec is number of levels within electronic 0 - so nEner is one
1991  * beyond end of array here - but will be decremented at start of loop
1992  * this starts at the highest energy wihtin X and moves down to lower energies */
1993  nEner = nLevels_per_elec[0];
1994  while( (--nEner) >= nXLevelsMatrix )
1995  {
1996 
1997  /* array of energy sorted indices within X - we are moving down
1998  * starting from highest level within X */
1999  ip = H2_ipX_ener_sort[nEner];
2000  iVib = ipVib_H2_energy_sort[ip];
2001  iRot = ipRot_H2_energy_sort[ip];
2002 
2003  if( nEner+1 < nLevels_per_elec[0] )
2005 
2006  /* >>chng 05 apr 30,GS, Instead of hmi.H2_total, the specific populations are used because high levels have much less
2007  * populations than ground levels which consists most of the H2 population.
2008  * only do this if working level is not v=0, J=0, 1 */
2009  if( nEner >1 )
2010  {
2011  H2_col_rate_out[iVib][iRot] +=
2012  /* H2 grain interactions
2013  * rate (s-1) all v,J levels go to 0 or 1 preserving spin */
2015 
2016  /* this goes into v=0, and J=0 or 1 depending on whether initial
2017  * state is ortho or para */
2018  H2_col_rate_in[0][H2_lgOrtho[0][iVib][iRot]] +=
2019  /* H2 grain interactions
2020  * rate (cm-3 s-1) all v,J levels go to 0 or 1 preserving spin,
2021  * in above lgOrtho says whether should go to 0 or 1 */
2023  }
2024  else if( nEner == 1 )
2025  {
2026  /* this is special J=1 to J=0 collision, which is only fast at
2027  * very low grain temperatures */
2028  H2_col_rate_out[0][1] +=
2029  /* H2 grain interactions
2030  * H2 ortho - para conversion on grain surface,
2031  * rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
2033 
2034  H2_col_rate_in[0][0] +=
2035  /* H2 grain interactions
2036  * H2 ortho - para conversion on grain surface,
2037  * rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
2039  }
2040 
2041  /* will become rate (cm-3 s-1) other levels have radiative transitions to here */
2042  H2_rad_rate_in[iVib][iRot] = 0.;
2043  H2_rad_rate_out[0][iVib][iRot] = 0.;
2044 
2045  /* the next two account for the Solomon process,
2046  * the first is the sum of decays from electronic excited into X
2047  * second is X going into all excited electronic states
2048  * units cm-3 s-1 */
2049  H2_rad_rate_in[iVib][iRot] += H2_X_rate_from_elec_excited[iVib][iRot];
2050 
2051  /* radiative & cosmic ray rates [s-1] to electronic excited states from X only vibration and rot */
2052  H2_rad_rate_out[0][iVib][iRot] += H2_X_rate_to_elec_excited[iVib][iRot];
2053 
2054  /* now sum over states within X which are higher than current state */
2055  iElecHi = 0;
2056  for( ipHi = nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi )
2057  {
2058  ip = H2_ipX_ener_sort[ipHi];
2059  iVibHi = ipVib_H2_energy_sort[ip];
2060  iRotHi = ipRot_H2_energy_sort[ip];
2061  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2062  /* the rate we enter this state from more highly excited states within X
2063  * by radiative decays, which have delta J = 0 or 2 */
2064  /* note test on vibration is needed - iVibHi<iVib, energy order ok and space not allocated */
2065  /* >>chng 05 feb 07, tried to use use lgH2_line_exists but cant */
2066  if( ( abs(iRotHi-iRot) ==2 || iRotHi==iRot ) && (iVib <= iVibHi) &&
2067  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].ipCont > 0 )
2068  {
2069  double rateone;
2070  rateone =
2071  H2_old_populations[iElecHi][iVibHi][iRotHi]*
2072  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Aul*
2073  (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pesc +
2074  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pelec_esc +
2075  H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pdest);
2076  ASSERT( rateone >=0 );
2077 
2078  /* units cm-3 s-1 */
2079  H2_rad_rate_in[iVib][iRot] += rateone;
2080  }
2081  }
2082 
2083  /* =======================INSIDE POPULATIONS X CONVERGE LOOP =====================*/
2084  /* we now have total rate this state is populated from above, now get rate
2085  * this state interacts with levels that are below */
2086  iElecLo = 0;
2087  for( ipLo = 0; ipLo<nEner; ++ipLo )
2088  {
2089  ip = H2_ipX_ener_sort[ipLo];
2090  iVibLo = ipVib_H2_energy_sort[ip];
2091  iRotLo = ipRot_H2_energy_sort[ip];
2092  /* radiative interactions between this level and lower levels */
2093  /* the test on vibration is needed - the energies are ok but the space does not exist */
2094  /* >>chng 05 feb 07, can't use lgH2_line_exists */
2095  if( ( abs(iRotLo-iRot) == 2 || iRotLo == iRot ) && (iVibLo <= iVib) &&
2096  H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].ipCont > 0 )
2097  {
2098  H2_rad_rate_out[0][iVib][iRot] +=
2099  H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Aul*
2100  (H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pesc +
2101  H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pelec_esc +
2102  H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pdest);
2103  }
2104  }
2105  /* =======================INSIDE X POPULATIONS CONVERGE LOOP =====================*/
2106 
2107  /* we now have the total rates into and out of this level, get its population
2108  * units cm-3 */
2109  H2_populations[iElec][iVib][iRot] =
2110  (H2_col_rate_in[iVib][iRot]+ H2_rad_rate_in[iVib][iRot]+H2_X_source[nEner]) /
2111  SDIV(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[0][iVib][iRot]+H2_X_sink[nEner]);
2112 
2113  ASSERT( H2_populations[iElec][iVib][iRot] >= 0. );
2114  }
2115  /* >>chng 05 may 10, move to following back substitution part within X */
2116  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2117  /* now do lowest levels H2_populations with matrix,
2118  * these should be collisionally dominated */
2119  if( nXLevelsMatrix )
2120  {
2122  /* the total abundance - frac_matrix is fraction of pop that was in these
2123  * levels the last time this was done */
2124  hmi.H2_total * (realnum)frac_matrix );
2125  }
2126  iElecHi = 0;
2128  {
2129  fprintf(ioQQQ," Rel pop(e=%li)" ,iElecHi);
2130  }
2131 
2132  /* find ortho and para densites, sum of pops in each vibration */
2133  /* this will become total pop is X, which will be renormed to equal hmi.H2_total */
2134  pops_per_elec[0] = 0.;
2135  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
2136  {
2137  double sumv;
2138  sumv = 0.;
2139  pops_per_vib[0][iVibHi] = 0.;
2140 
2141  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
2142  {
2143  pops_per_elec[0] += H2_populations[iElecHi][iVibHi][iRotHi];
2144  sumv += H2_populations[iElecHi][iVibHi][iRotHi];
2145  pops_per_vib[0][iVibHi] += H2_populations[iElecHi][iVibHi][iRotHi];
2146  }
2147  /* print sum of H2_populations in each vibration if trace on */
2149  fprintf(ioQQQ,"\t%.2e",sumv/hmi.H2_total);
2150  }
2152  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2154  {
2155  fprintf(ioQQQ,"\n");
2156  /* print the ground vibration state */
2157  fprintf(ioQQQ," Rel pop(0,J)");
2158  iElecHi = 0;
2159  iVibHi = 0;
2160  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
2161  {
2162  fprintf(ioQQQ,"\t%.2e",H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total);
2163  }
2164  fprintf(ioQQQ,"\n");
2165  }
2166 
2167  /* now find population in states done with matrix - this is only used to pass
2168  * to matrix solver */
2169  iElec = 0;
2170  sum_pops_matrix = 0.;
2171  ip =0;
2172  for( i=0; i<nXLevelsMatrix; ++i )
2173  {
2174  ip = H2_ipX_ener_sort[i];
2175  iVib = ipVib_H2_energy_sort[ip];
2176  iRot = ipRot_H2_energy_sort[ip];
2177  sum_pops_matrix += H2_populations[iElec][iVib][iRot];
2178  }
2179  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2180  /* this is self consistent since pops_per_elec[0] came from current soln,
2181  * as did the matrix. pops will be renormalized by results from the chemistry
2182  * a few lines down */
2183  frac_matrix = sum_pops_matrix / SDIV(pops_per_elec[0]);
2184 
2185  /* assuming that all H2 population is in X, this is the
2186  * ratio of H2 that came out of the chemistry network to what we just obtained -
2187  * we need to multiply the pops by renorm to agree with the chemistry,
2188  * this routine does not alter hmi.H2_total, but does change pops_per_elec */
2189  H2_renorm_conserve = hmi.H2_total/ SDIV(pops_per_elec[0]);
2190  /*pops_per_elec[0] = hmi.H2_total;*/
2191 
2192  /* renormalize H2_populations - the H2_populations were updated by renorm when routine entered,
2193  * before pops determined - but population determinations above do not have a sum rule on total
2194  * population - this renorm is to preserve total population */
2195  H2_sum_excit_elec_den = 0.;
2196  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
2197  {
2198  pops_per_elec[iElecHi] *= H2_renorm_conserve;
2199  if( iElecHi > 0 )
2201 
2202  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
2203  {
2204  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
2205  {
2206  H2_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve;
2207  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2208  }
2209  }
2210  }
2211 
2212  /* this loop first checks for largest changes in populations, to determine whether
2213  * we have converged, then updates the population array with a new value,
2214  * which may be a mean of old and new
2215  * update populations check convergence converged */
2216  sumold = 0.;
2217  for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
2218  {
2219  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
2220  {
2221  for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
2222  {
2223  double rel_change;
2224  /* keep track of largest relative change in H2_populations to
2225  * determines convergence */
2226  if( fabs(H2_populations[iElec][iVib][iRot] -
2227  H2_old_populations[iElec][iVib][iRot])/
2228  /* on first call some very high J states can have zero pop ,
2229  * hence the SDIV, will retain sign for checks on oscilations,
2230  * hence the fabs */
2231  SDIV(H2_populations[iElec][iVib][iRot]) > fabs(PopChgMax_relative) &&
2232  /* >>chng 03 jul 19, this had simply been H2_populations > SMALLFLOAT,
2233  * change to relative pops > 1e-15, spent too much time converging
2234  * levels at pops = 1e-37 */
2235  /* >>chng 03 dec 27, from rel pop 1e-15 to 1e-6 since converging heating will
2236  * be main convergence criteria check convergence */
2237  /*H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(hmi.H2_total)>1e-15 )*/
2238  H2_populations[iElec][iVib][iRot]/SDIV(hmi.H2_total)>1e-6 )
2239  {
2240  PopChgMax_relative =
2241  (H2_populations[iElec][iVib][iRot] -
2242  H2_old_populations[iElec][iVib][iRot])/
2243  SDIV(H2_populations[iElec][iVib][iRot]);
2244  iRotMaxChng_relative = iRot;
2245  iVibMaxChng_relative = iVib;
2246  popold_relative = H2_old_populations[iElec][iVib][iRot];
2247  popnew_relative = H2_populations[iElec][iVib][iRot];
2248  }
2249  /* >>chng 05 feb 08, add largest rel change in total, this will be converged
2250  * down to higher accuracy than above
2251  * keep track of largest change in H2_populations relative to total H2 to
2252  * determine convergence check convergence */
2253  rel_change = (H2_populations[iElec][iVib][iRot] -
2254  H2_old_populations[iElec][iVib][iRot])/SDIV(hmi.H2_total);
2255  /* retain sign for checks on oscillations hence the fabs */
2256  if( fabs(rel_change) > fabs(PopChgMax_total) )
2257  {
2258  PopChgMax_total = rel_change;
2259  iRotMaxChng_total = iRot;
2260  iVibMaxChng_total = iVib;
2261  popold_total = H2_old_populations[iElec][iVib][iRot];
2262  popnew_total = H2_populations[iElec][iVib][iRot];
2263  }
2264 
2265  kase = -1;
2266  /* update populations - we used the old populations to update the
2267  * current new populations - will do another iteration if they changed
2268  * by much. here old populations are updated for next sweep through molecule */
2269  /* pop oscillations have occurred - use small changes */
2270  /* >>chng 04 may 10, turn this back on - now with min on how small frac new
2271  * can become */
2272  rel_change = fabs( H2_old_populations[iElec][iVib][iRot] -
2273  H2_populations[iElec][iVib][iRot] )/
2274  SDIV( H2_populations[iElec][iVib][iRot] );
2275 
2276  /* this branch very large changes, use mean of logs but onlly if both are positive*/
2277  if( rel_change > 3. &&
2278  H2_old_populations[iElec][iVib][iRot]*H2_populations[iElec][iVib][iRot]>0 )
2279  {
2280  /* large changes or oscillations - take average in the log */
2281  H2_old_populations[iElec][iVib][iRot] = pow( 10. ,
2282  log10(H2_old_populations[iElec][iVib][iRot])/2. +
2283  log10(H2_populations[iElec][iVib][iRot])/2. );
2284  kase = 2;
2285  }
2286 
2287  /* modest change, use means of old and new */
2288  else if( rel_change> 0.1 )
2289  {
2290  realnum frac_old=0.25f;
2291  /* large changes or oscillations - take average */
2292  H2_old_populations[iElec][iVib][iRot] =
2293  frac_old*H2_old_populations[iElec][iVib][iRot] +
2294  (1.f-frac_old)*H2_populations[iElec][iVib][iRot];
2295  kase = 3;
2296  }
2297  else
2298  {
2299  /* small changes, use new value */
2300  H2_old_populations[iElec][iVib][iRot] =
2301  H2_populations[iElec][iVib][iRot];
2302  kase = 4;
2303  }
2304  sumold += H2_old_populations[iElec][iVib][iRot];
2305  }
2306  }
2307  }
2308  /* will renormalize so that total population is correct */
2309  H2_renorm_conserve_init = hmi.H2_total/sumold;
2310 
2311  /* renormalize H2_populations - the H2_populations were updated by renorm when routine entered,
2312  * before pops determined - but population determinations above do not have a sum rule on total
2313  * population - this renorm is to preserve total population */
2314  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
2315  {
2316  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
2317  {
2318  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
2319  {
2320  H2_old_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve_init;
2321  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2322  }
2323  }
2324  }
2325  /* get current ortho-para ratio, will be used as test on convergence */
2326  iElecHi = 0;
2327  h2.ortho_density = 0.;
2328  h2.para_density = 0.;
2329  H2_den_s = 0.;
2330  H2_den_g = 0.;
2331 
2332  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
2333  {
2334  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
2335  {
2336  /* find current population in H2s and H2g */
2337  if( energy_wn[0][iVibHi][iRotHi]> ENERGY_H2_STAR )
2338  {
2339  H2_den_s += H2_populations[iElecHi][iVibHi][iRotHi];
2340  }
2341  else
2342  {
2343  H2_den_g += H2_populations[iElecHi][iVibHi][iRotHi];
2344  }
2345  if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] )
2346  {
2347  h2.ortho_density += H2_populations[iElecHi][iVibHi][iRotHi];
2348  }
2349  else
2350  {
2351  h2.para_density += H2_populations[iElecHi][iVibHi][iRotHi];
2352  }
2353  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2354  }
2355  }
2356  /* these will be used to determine whether solution has converged */
2357  ortho_para_older = ortho_para_old;
2358  ortho_para_old = ortho_para_current;
2359  ortho_para_current = h2.ortho_density / SDIV( h2.para_density );
2360 
2361  /* this will be evaluated in call to routine that follows - will check
2362  * whether this has converged */
2363  old_solomon_rate = hmi.H2_Solomon_dissoc_rate_BigH2_H2g;
2364 
2365  /* >>chng 05 jul 24, break code out into separate routine for clarify
2366  * located in mole_h2_etc.c - true says to only do Solomon rate */
2367  H2_Solomon_rate( );
2368 
2369  /* are changes too large? must decide whether population shave converged,
2370  * will check whether H2_populations themselves have changed by much,
2371  * but also change in heating by collisional deexcitation is stable */
2372  HeatChangeOld = HeatChange;
2373  HeatChange = old_HeatH2Dexc_BigH2 - hmi.HeatH2Dexc_BigH2;
2374  {
2375  static long int loop_h2_oscil=-1;
2376  /* check whether pops are oscillating, as evidenced by change in
2377  * heating changing sign */
2378  if( loop_h2_pops>2 && (
2379  (HeatChangeOld*HeatChange<0. ) ||
2380  (PopChgMax_relative*PopChgMaxOld_relative<0. ) ) )
2381  {
2382  lgH2_pops_oscil = true;
2383  if( loop_h2_pops > 6 )
2384  {
2385  loop_h2_oscil = loop_h2_pops;
2386  lgH2_pops_ever_oscil = true;
2387  /* make this smaller in attempt to damp out oscillations,
2388  * but don't let get too small*/
2389  frac_new_oscil *= 0.8f;
2390  frac_new_oscil = MAX2( frac_new_oscil , 0.1f);
2391  ++n_pop_oscil;
2392  }
2393  }
2394  else
2395  {
2396  lgH2_pops_oscil = false;
2397  /* turn off flag if no oscillations for a while */
2398  if( loop_h2_pops - loop_h2_oscil > 4 )
2399  {
2400  frac_new_oscil = 1.f;
2401  lgH2_pops_ever_oscil = false;
2402  }
2403  }
2404  }
2405 
2406  /* reevaluate heating - cooling if H2 molecule is significant source or either,
2407  * since must have stable heating cooling rate */
2408  old_HeatH2Dexc_BigH2 = hmi.HeatH2Dexc_BigH2;
2410  hmi.HeatH2Dexc_BigH2==0. )
2411  H2_Cooling("H2lup");
2412 
2413  /* begin check on whether solution is converged */
2414  lgConv_h2_soln = true;
2415  lgPopsConv_total = true;
2416  lgPopsConv_relative = true;
2417  lgHeatConv = true;
2418  lgSolomonConv = true;
2419  lgOrthoParaRatioConv = true;
2420 
2421  /* these are all the convergence tests
2422  * check convergence converged */
2423  if( fabs(PopChgMax_relative)>converge_pops_relative )
2424  {
2425  /*lgPopsConv = (fabs(PopChgMax_relative)<=0.1);*/
2426  lgConv_h2_soln = false;
2427  lgPopsConv_relative = false;
2428  /* >>chng 04 sep 08, set quant_new to new chng max gs */
2429  /*quant_old = PopChgMax_relative;*/
2430  quant_old = PopChgMaxOld_relative;
2431  /*quant_new = 0.;*/
2432  quant_new = PopChgMax_relative;
2433 
2434  strcpy( chReason , "rel pops changed" );
2435  }
2436 
2437  /* check largest change in a level population relative to total h2
2438  * population convergence converged check */
2439  else if( fabs(PopChgMax_total)>converge_pops_total)
2440  {
2441  lgConv_h2_soln = false;
2442  lgPopsConv_total = false;
2443  /* >>chng 04 sep 08, set quant_new to new chng max gs */
2444  /*quant_old = PopChgMax_relative;*/
2445  quant_old = PopChgMaxOld_total;
2446  /*quant_new = 0.;*/
2447  quant_new = PopChgMax_total;
2448 
2449  strcpy( chReason , "tot pops changed" );
2450  }
2451 
2452  /* >>chng 04 apr 30, look at change in ortho-para ratio, also that is not
2453  * oscillating */
2454  /* >>chng 04 dec 15, only look at change, and don't make allowed change so tiny -
2455  * these were attempts at fixing problems that were due to shielding not thin*/
2456  else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> converge_ortho_para )
2457  /* else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> 1e-3
2458  && (ortho_para_current-ortho_para_old)*(ortho_para_old-ortho_para_older)>0. )*/
2459  {
2460  lgConv_h2_soln = false;
2461  lgOrthoParaRatioConv = false;
2462  quant_old = ortho_para_old;
2463  quant_new = ortho_para_current;
2464  strcpy( chReason , "ortho/para ratio changed" );
2465  }
2466  /* >>chng 04 dec 16, reduce error allowed fm /5 to /2, to be similar to
2467  * logic in conv_base */
2468  else if( !thermal.lgTemperatureConstant &&
2469  fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/MAX2(thermal.ctot,thermal.htot) >
2471  /* >>chng 04 may 09, do not check on error in heating if constant temperature */
2472  /*&& !(thermal.lgTemperatureConstant || phycon.te <= phycon.TEMP_LIMIT_LOW )*/ )
2473  {
2474  /* default on HeatCoolRelErrorAllowed is 0.02 */
2475  /*lgHeatConv = (fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot <=
2476  * conv.HeatCoolRelErrorAllowed/5.);*/
2477  lgConv_h2_soln = false;
2478  lgHeatConv = false;
2479  quant_old = old_HeatH2Dexc_BigH2/MAX2(thermal.ctot,thermal.htot);
2481  strcpy( chReason , "heating changed" );
2482  /*fprintf(ioQQQ,"DEBUG old new trip \t%.4e \t %.4e\n",
2483  old_HeatH2Dexc_BigH2,
2484  hmi.HeatH2Dexc_BigH2);*/
2485  }
2486 
2487  /* check on Solomon rate,
2488  * >>chng 04 aug 28, do not do this check if induced processes are disabled,
2489  * since Solomon process is then irrelevant */
2490  /* >>chng 04 sep 21, GS*/
2491  else if( rfield.lgInducProcess &&
2492  /* this is check that H2 abundance has not been set - if it has been
2493  * then we don't care what the Solomon rate is doing */
2494  hmi.H2_frac_abund_set==0 &&
2495  /*>>chng 05 feb 10, rather than checking change in Solomon relative to Solomon,
2496  * check it relative to total h2 destruction rate */
2497  fabs( hmi.H2_Solomon_dissoc_rate_BigH2_H2g - old_solomon_rate)/SDIV(hmi.H2_rate_destroy) >
2499  {
2500  lgConv_h2_soln = false;
2501  lgSolomonConv = false;
2502  quant_old = old_solomon_rate;
2504  strcpy( chReason , "Solomon rate changed" );
2505  }
2506 
2507  /* did we pass all the convergence test */
2508  if( !lgConv_h2_soln )
2509  {
2510  /* this branch H2 H2_populations within X are not converged,
2511  * print diagnostic */
2512 
2514  {
2515  /*fprintf(ioQQQ,"temppp\tnew\t%.4e\tnew\t%.4e\t%.4e\n",
2516  hmi.HeatH2Dexc_BigH2,
2517  old_HeatH2Dexc_BigH2,
2518  fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot );*/
2519  fprintf(ioQQQ," loop %3li no conv oscl?%c why:%s ",
2520  loop_h2_pops,
2521  TorF(lgH2_pops_ever_oscil),
2522  chReason );
2523  if( !lgPopsConv_relative )
2524  fprintf(ioQQQ," PopChgMax_relative:%.4e v:%li J:%li old:%.4e new:%.4e",
2525  PopChgMax_relative,
2526  iVibMaxChng_relative,
2527  iRotMaxChng_relative ,
2528  popold_relative ,
2529  popnew_relative );
2530  else if( !lgPopsConv_total )
2531  fprintf(ioQQQ," PopChgMax_total:%.4e v:%li J:%li old:%.4e new:%.4e",
2532  PopChgMax_total,
2533  iVibMaxChng_total,
2534  iRotMaxChng_total ,
2535  popold_total ,
2536  popnew_total );
2537  else if( !lgHeatConv )
2538  fprintf(ioQQQ," heat:%.4e old:%.4e new:%.4e",
2539  (hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/MAX2(thermal.ctot,thermal.htot),
2540  quant_old ,
2541  quant_new);
2542  /* Solomon rate changed */
2543  else if( !lgSolomonConv )
2544  fprintf(ioQQQ," d(sol rate)/tot dest\t%2e",(old_solomon_rate - hmi.H2_Solomon_dissoc_rate_BigH2_H2g)/SDIV(hmi.H2_rate_destroy));
2545  else if( !lgOrthoParaRatioConv )
2546  fprintf(ioQQQ," current, old, older ratios are %.4e %.4e %.4e",
2547  ortho_para_current , ortho_para_old, ortho_para_older );
2548  else
2549  TotalInsanity();
2550  fprintf(ioQQQ,"\n");
2551  }
2552  }
2553  /* end convergence criteria */
2554 
2555  /*fprintf(ioQQQ,"DEBUG h2 heat\t%3li\t%.2f\t%.4e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2556  loop_h2_pops,
2557  fnzone,
2558  phycon.te,
2559  dense.eden,
2560  hmi.HeatH2Dexc_BigH2,
2561  hmi.HeatH2Dexc_BigH2/thermal.ctot ,
2562  hmi.H2_total,
2563  H2_renorm_chemistry ,
2564  H2_renorm_conserve,
2565  hmi.H2_H2g_to_H2s_rate_BigH2);*/
2566  if( trace.nTrConvg >= 5 )
2567  {
2568  fprintf( ioQQQ,
2569  " H2 5lev %li Conv?%c",
2570  loop_h2_pops ,
2571  TorF(lgConv_h2_soln) );
2572 
2573  if( fabs(PopChgMax_relative)>0.1 )
2574  fprintf(ioQQQ," pops, rel chng %.3e",PopChgMax_relative);
2575  else
2576  fprintf(ioQQQ," rel heat %.3e rel chng %.3e H2 heat/cool %.2e",
2578  fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot ,
2580 
2581  fprintf( ioQQQ,
2582  " Oscil?%c Ever Oscil?%c",
2583  TorF(lgH2_pops_oscil) ,
2584  TorF(lgH2_pops_ever_oscil) );
2585  if( lgH2_pops_ever_oscil )
2586  fprintf(ioQQQ," frac_new_oscil %.4f",frac_new_oscil);
2587  fprintf(ioQQQ,"\n");
2588  }
2589 
2590  if( mole.nH2_TRACE >= mole.nH2_trace_full )
2591  {
2592  fprintf(ioQQQ,
2593  "H2 loop\t%li\tkase pop chng\t%i\tchem renorm fac\t%.4e\tortho/para ratio:\t%.3e\tfrac of pop in matrix: %.3f\n",
2594  loop_h2_pops,
2595  kase,
2598  frac_matrix);
2599 
2600  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
2601  if( iVibMaxChng_relative>=0 && iRotMaxChng_relative>=0 )
2602  fprintf(ioQQQ,
2603  "end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n",
2604  loop_h2_pops,
2605  PopChgMax_relative ,
2606  H2_old_populations[0][iVibMaxChng_relative][iRotMaxChng_relative],
2607  H2_populations[0][iVibMaxChng_relative][iRotMaxChng_relative],
2608  iVibMaxChng_relative , iRotMaxChng_relative
2609  );
2610  }
2611  }
2612  /* =======================END POPULATIONS CONVERGE LOOP =====================*/
2613 
2614  /* evaluate H2 rates over H2g and H2s for use in chemistry network */
2615  H2_gs_rates();
2616 
2617  /* >>chng 05 feb 08, do not print if we are in search phase */
2618  if( !lgConv_h2_soln && !conv.lgSearch )
2619  {
2620  conv.lgConvPops = false;
2621  strcpy( conv.chConvIoniz, "H2 pop cnv" );
2622  fprintf(ioQQQ,
2623  " H2_LevelPops: H2_populations not converged in %li tries; due to %s, old, new are %.4e %.4e, iteration %li zone %.2f.\n",
2624  loop_h2_pops,
2625  chReason,
2626  quant_old,
2627  quant_new ,
2628  iteration ,
2629  fnzone );
2630  ConvFail("pops","H2");
2631  }
2632 
2633  /* loop over all possible lines and set H2_populations,
2634  * and quantities that depend on them */
2635  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
2636  {
2637  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
2638  {
2639  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
2640  {
2641  long int lim_elec_lo = 0;
2642  double H2popHi = H2_populations[iElecHi][iVibHi][iRotHi];
2643  // this will update Lo->Pop as well as Lo and Hi point to the same set of data structures
2644  // the loops below are OK since Lo->Pop is only accessed after Hi->Pop has been set.
2645  H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop = H2popHi;
2646  ASSERT( H2popHi >= 0. );
2647  realnum H2gHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->g;
2648  /* now the lower levels */
2649  /* NB - X is the only lower level considered here, since we are only
2650  * concerned with excited electronic levels as a photodissociation process
2651  * code exists to relax this assumption - simply change following to iElecHi */
2652  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
2653  {
2654  /* want to include all vibration states in lower level if different electronic level,
2655  * but only lower vibration levels if same electronic level */
2656  long int nv = h2.nVib_hi[iElecLo];
2657  if( iElecLo==iElecHi )
2658  nv = iVibHi;
2659  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
2660  {
2661  long nr = h2.nRot_hi[iElecLo][iVibLo];
2662  if( iElecLo==iElecHi && iVibHi==iVibLo )
2663  nr = iRotHi-1;
2664 
2665  mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
2666  mt6i H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]);
2667  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
2668  {
2669  if( *lgH2le++ )
2670  {
2671  /* following two heat exchange excitation, deexcitation */
2672  H2L->Coll.cool = 0.;
2673  H2L->Coll.heat = 0.;
2674 
2675  H2L->Emis->PopOpc =
2676  H2L->Lo->Pop - H2popHi * H2L->Lo->g / H2gHi;
2677 
2678  /* number of photons in the line */
2679  H2L->Emis->phots = H2L->Emis->Aul *
2680  (H2L->Emis->Pesc + H2L->Emis->Pelec_esc) *
2681  H2popHi;
2682 
2683  /* intensity of line */
2684  H2L->Emis->xIntensity =
2685  H2L->Emis->phots *
2686  H2L->EnergyErg;
2687 
2688  if( iElecHi==0 )
2689  {
2691  /* the ground electronic state, most excitations are not direct pumping
2692  * (rather indirect, which does not count for ColOvTot) */
2693  H2L->Emis->ColOvTot = 1.;
2694  }
2695  else
2696  {
2697  /* these are excited electronic states, mostly pumped, except for supras */
2699  H2L->Emis->ColOvTot = 0.;
2700  }
2701  }
2702  ++H2L;
2703  }
2704  }
2705  }
2706  }
2707  }
2708  }
2709 
2710  /* add up H2 + hnu => 2H, continuum photodissociation,
2711  * this is not the Solomon process, true continuum */
2712  /* >>chng 05 jun 16, GS, add dissociation to triplet states*/
2717  hmi.H2_BigH2_H2g_av = 0.;
2718  hmi.H2_BigH2_H2s_av = 0.;
2719  /* >>chng 05 jul 20, GS, add dissociation by H2 g and H2s*/
2722 
2723  iElec = 0;
2724  H2_BigH2_H2s = 0.;
2725  H2_BigH2_H2g = 0.;
2726  hmi.H2g_BigH2 =0;
2727  hmi.H2s_BigH2 = 0;
2728  hmi.H2_total_BigH2 =0;
2729  hmi.H2g_LTE_bigH2 =0.;
2730  hmi.H2s_LTE_bigH2 = 0.;
2731  /* >>chng 05 oct 20, no need to reset this var here */
2732  /*exp_disoc = sexp(H2_DissocEnergies[0]/phycon.te_wn);*/
2733 
2734  /* >>chng 05 sep 12, TE, define a cutoff wavelength of 800 Angstrom
2735  * this is chosen as the cross sections given by
2736  *>>refer H2 photo cs Allison, A.C. & Dalgarno, A. 1969, Atomic Data, 1, 91
2737  * show a sharp decline in the cross section*/
2738  {
2739  static long ip_cut_off = -1;
2740  if( ip_cut_off < 0 )
2741  {
2742  /* one-time initialization of this pointer */
2743  ip_cut_off = ipoint( 1.14 );
2744  }
2745 
2746  /* >>chng 05 sep 12, TE, assume all H2s is at 2.5 eV
2747  * the dissociation threshold is at 1.07896 Rydberg*/
2748  flux_accum_photodissoc_BigH2_H2s = 0;
2749  ip_H2_level = ipoint( 1.07896 - 2.5 / EVRYD);
2750  for( i= ip_H2_level; i < ip_cut_off; ++i )
2751  {
2752  flux_accum_photodissoc_BigH2_H2s += ( rfield.flux[i-1] + rfield.ConInterOut[i-1]+
2753  rfield.outlin[i-1]+ rfield.outlin_noplot[i-1] );
2754  }
2755 
2756  /* sum over all levels to obtain s and g populations and dissociation rates */
2757  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
2758  {
2759  for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
2760  {
2761  /* >>chng 05 mar 22, TE, moved H2_photodissoc_BigH2_H2s in this statement and divide by the
2762  * density of H2s not total H2, we consider direct photodissociation only for H2s */
2763  /* >>chng 05 mar 22, TE, this should be for H2* rather than total */
2764  /* this is the total rate of direct photo-dissociation of excited electronic states into
2765  * the X continuum - this is continuum photodissociation, not the Solomon process */
2766  /* >>chng 03 sep 03, make sum of pops of excited states */
2767  if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
2768  {
2769  double arg_ratio;
2771  H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2s;
2772 
2773  /* cosmic ray & secondary electron excitation to triplets
2774  * physics described where similar process is done for
2775  * big molecule
2776  *>>chng 07 apr 08, from 3 to 10 to better capture results
2777  * of Dalgarno et al 99 */
2779  H2_populations[iElec][iVib][iRot] * 10.f*secondaries.x12tot;
2780 
2781  /* sum of pops in levels in H2* for use in chemistry network */
2782  H2_BigH2_H2s += H2_populations[iElec][iVib][iRot];
2783 
2784  /* >>chng 05 jun 28, TE, determine average energy level in H2s */
2785  hmi.H2_BigH2_H2s_av += (H2_populations[iElec][iVib][iRot] * energy_wn[0][iVib][iRot]);
2786 
2787  /* >>chng 05 july 20, GS, collisional dissociation by H2s, unit s-1*/
2788  hmi.Average_collH2s_dissoc += H2_populations[iElec][iVib][iRot] * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
2789 
2790  /* >>chng 05 oct 17, GS, LTE populations of H2s*/
2791  arg_ratio = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]);
2792  if( arg_ratio > 0. )
2793  {
2794  /* >>chng 05 oct 21, GS, only add ratio if Boltzmann factor > 0 */
2795  hmi.H2s_LTE_bigH2 += H2_populations[0][iVib][iRot]*SAHA/SDIV(phycon.te32*arg_ratio)*
2796  (H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5;
2797  }
2798  }
2799  else
2800  {
2801  double arg_ratio;
2802  /* >>chng 05 sep 12, TE, for H2g do the sum explicitly for every level*/
2803  flux_accum_photodissoc_BigH2_H2g = 0;
2804  /* this is the dissociation energy needed for the level*/
2805  ip_H2_level = ipoint( 1.07896 - energy_wn[0][iVib][iRot] * WAVNRYD);
2806 
2807  for( i= ip_H2_level; i < ip_cut_off; ++i )
2808  {
2809  flux_accum_photodissoc_BigH2_H2g += ( rfield.flux[i-1] + rfield.ConInterOut[i-1]+
2810  rfield.outlin[i-1]+ rfield.outlin_noplot[i-1] );
2811  }
2812 
2814  H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2g;
2815 
2816 
2817  /* cosmic ray & secondary electron excitation to triplets
2818  * physics described where similar process is done for
2819  * big molecule
2820  *>>chng 07 apr 08, from 3 to 10 to better capture results
2821  * of Dalgarno et al 99 */
2823  H2_populations[iElec][iVib][iRot] * 10.f*secondaries.x12tot;
2824 
2825  /* sum of pops in levels in H2g for use in chemistry network */
2826  H2_BigH2_H2g += H2_populations[iElec][iVib][iRot];
2827 
2828  /* >>chng 05 jun 28, TE, determine average energy level in H2g */
2829  hmi.H2_BigH2_H2g_av += (H2_populations[iElec][iVib][iRot] * energy_wn[0][iVib][iRot]);
2830 
2831  /* >>chng 05 jul 20, GS, collisional dissociation by H2s, unit s-1*/
2832  hmi.Average_collH2g_dissoc += H2_populations[iElec][iVib][iRot] * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
2833 
2834  /* >>chng 05 oct 17, GS, LTE populations of H2g*/
2835  arg_ratio = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]);
2836  if( arg_ratio > 0. )
2837  {
2838  hmi.H2g_LTE_bigH2 += H2_populations[0][iVib][iRot]*(SAHA/SDIV(phycon.te32*arg_ratio)*
2839  (H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5);
2840  }
2841  }
2842  }
2843  }
2844  }
2845  hmi.H2g_BigH2 = (realnum)H2_BigH2_H2g;
2846  hmi.H2s_BigH2 = (realnum)H2_BigH2_H2s;
2848 
2849  /* average energy in H2s */
2850  hmi.H2_BigH2_H2s_av = hmi.H2_BigH2_H2s_av / H2_BigH2_H2s;
2851  /* average energy in H2g */
2852  hmi.H2_BigH2_H2g_av = hmi.H2_BigH2_H2g_av / H2_BigH2_H2g;
2853 
2854  /* above sum was rate per unit vol since mult by H2 density, now div by H2* density to get rate s-1 */
2855  /* 0.25e-18 is wild guess of typical photodissociation cross section, from
2856  * >>refer H2 dissoc Allison, A.C. & Dalgarno, A. 1969, Atomic Data, 1, 91
2857  * this is based on an average of the highest v values they gave. unfortunately, we want
2858  * the highest J values -
2859  * final units are s-1*/
2864  hmi.Average_collH2g_dissoc = hmi.Average_collH2g_dissoc /SDIV(H2_BigH2_H2g);/* unit cm3s-1*/
2865  hmi.Average_collH2s_dissoc = hmi.Average_collH2s_dissoc /SDIV(H2_BigH2_H2s);/* unit cm3s-1*/
2866  hmi.H2s_LTE_bigH2 = hmi.H2s_LTE_bigH2/SDIV(H2_BigH2_H2s);
2867  hmi.H2g_LTE_bigH2 = hmi.H2g_LTE_bigH2/SDIV(H2_BigH2_H2g);
2868 
2869  /* >>chng 05 jul 09, GS*/
2870  /* average Einstein value for H2* to H2g, GS*/
2871  double sumpop1 = 0.;
2872  double sumpopA1 = 0.;
2873  double sumpopcollH2O_deexcit = 0.;
2874  double sumpopcollH2p_deexcit = 0.;
2875  double sumpopcollH_deexcit = 0.;
2876  double popH2s = 0.;
2877  double sumpopcollH2O_excit = 0.;
2878  double sumpopcollH2p_excit = 0.;
2879  double sumpopcollH_excit = 0.;
2880  double popH2g = 0.;
2881 
2882 
2883  iElecLo = 0;
2884  for( iVibHi=0; iVibHi<=h2.nVib_hi[0]; ++iVibHi )
2885  {
2886  long nr1 = h2.nRot_hi[0][iVibHi];
2887  for( iRotHi=h2.Jlowest[0]; iRotHi<=nr1; ++iRotHi )
2888  {
2889  double ewnHi = energy_wn[0][iVibHi][iRotHi];
2890  for( iVibLo=0; iVibLo<=h2.nVib_hi[0]; ++iVibLo )
2891  {
2892  long nr2 = h2.nRot_hi[0][iVibLo];
2893  md3ci ewnLo = energy_wn.ptr(0,iVibLo,h2.Jlowest[0]);
2894  for( iRotLo=h2.Jlowest[0]; iRotLo<=nr2; ++iRotLo )
2895  {
2896  double ewnLo2 = energy_wn[0][iVibLo][iRotLo];
2897  /*if( ewnHi > ENERGY_H2_STAR && *ewnLo++ < ENERGY_H2_STAR )*/
2898  if( ewnHi > ENERGY_H2_STAR && ewnLo2 < ENERGY_H2_STAR )
2899  {
2900  /* >>chng 05 jul 10, GS*/
2901  /* average collisional rate for H2* to H2g, GS*/
2902  if( H2_lgOrtho[0][iVibHi][iRotHi] == H2_lgOrtho[0][iVibLo][iRotLo] )
2903  {
2904  /* sums of populations */
2905  popH2s += H2_populations[0][iVibHi][iRotHi];
2906  popH2g += H2_populations[0][iVibLo][iRotLo];
2907 
2908  /* sums of deexcitation rates - H2* to H2g */
2909  sumpopcollH_deexcit += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0];
2910  sumpopcollH2O_deexcit += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][2];
2911  sumpopcollH2p_deexcit += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][3];
2912 
2913  /* sums of excitation rates - H2g to H2* */
2914  sumpopcollH_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] *
2915  H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
2916  sumpopcollH2O_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][2]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] *
2917  H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
2918  sumpopcollH2p_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][3]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] *
2919  H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
2920 
2921 
2922 
2923  /* if( (abs((iRotHi-iRotLo) ))==2 || (iRotHi==iRotLo ) ) */
2924  if( lgH2_line_exists[0][iVibHi][iRotHi][0][iVibLo][iRotLo] )
2925  {
2926  sumpop1 += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Hi->Pop;
2927  sumpopA1 += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Hi->Pop*
2928  H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Aul;
2929  }
2930  }
2931  }
2932  }
2933  }
2934  }
2935  }
2936  hmi.Average_A = sumpopA1/SDIV(sumpop1);
2937 
2938  /* collisional excitation and deexcitation of H2g and H2s */
2939  hmi.Average_collH2_deexcit = (sumpopcollH2O_deexcit+sumpopcollH2p_deexcit)/SDIV(popH2s);
2940  hmi.Average_collH2_excit = (sumpopcollH2O_excit+sumpopcollH2p_excit)/SDIV(popH2g);
2941  hmi.Average_collH_excit = sumpopcollH_excit/SDIV(popH2g);
2942  hmi.Average_collH_deexcit = sumpopcollH_deexcit/SDIV(popH2s);
2943  /* populations are badly off during search phase */
2944  if( conv.lgSearch )
2945  {
2946  hmi.Average_collH_excit /=10.;
2948  }
2949 
2950  /*fprintf(ioQQQ,
2951  "DEBUG Average_collH_excit sumpop = %.2e %.2e %.2e %.2e %.2e %.2e \n",
2952  popH2g,popH2s,sumpopcollH_deexcit ,sumpopcollH_excit ,
2953  sumpopcollH_deexcit/SDIV(popH2s) ,sumpopcollH_excit/SDIV(popH2g));*/
2954  /*fprintf(ioQQQ,"sumpop = %le sumpopA = %le Av= %le\n",
2955  sumpop1,sumpopA1 , hmi.Average_A );*/
2956 
2958  {
2959  fprintf(ioQQQ," H2_LevelPops exit2 Sol dissoc %.2e (TH85 %.2e)",
2963 
2964  /* Solomon process rate from X into the X continuum with units s-1
2965  * rates are total rate, and rates from H2g and H2s */
2966  fprintf(ioQQQ," H2g Sol %.2e H2s Sol %.2e",
2969 
2970  /* photoexcitation from H2g to H2s */
2971  fprintf(ioQQQ," H2g->H2s %.2e (TH85 %.2e)",
2974 
2975  /* add up H2s + hnu => 2H, continuum photodissociation,
2976  * this is not the Solomon process, true continuum, units s-1 */
2977  fprintf(ioQQQ," H2 con diss %.2e (TH85 %.2e)\n",
2980  }
2981  else if( mole.nH2_TRACE )
2982  {
2983  fprintf(ioQQQ," H2_LevelPops exit1 %8.2f loops:%3li H2/H:%.3e Sol dis old %.3e new %.3e",
2984  fnzone ,
2985  loop_h2_pops ,
2987  old_solomon_rate,
2989  fprintf(ioQQQ,"\n");
2990  }
2991 
2992  /* >>chng 03 sep 01, add this population - before had just used H2star from chem network */
2993  /* if big H2 molecule is turned on and used for this zone, use its
2994  * value of H2* (pops of all states with v > 0 ) rather than simple network */
2995 
2996  /* update number of times we have been called */
2998 
2999  /* this will say how many times the large H2 molecule has been called in this zone -
3000  * if not called (due to low H2 abundance) then not need to update its line arrays */
3002 
3003  /* >>chng 05 jun 21,
3004  * during search phase we want to use full matrix - save number of levels so that
3005  * we can restore it */
3006  nXLevelsMatrix = nXLevelsMatrix_save;
3007 
3008  /* >>chng 05 jan 19, check how many levels should be in the matrix if first call on
3009  * new zone, and we have a solution */
3010  /* end loop setting very first LTE H2_populations */
3012  {
3013  /* this is fraction of populations to include in matrix */
3014 # define FRAC 0.99999
3015  /* this loop is over increasing energy */
3016  double sum_pop = 0.;
3017  nEner = 0;
3018  iElec = 0;
3019 # define PRT false
3020  if( PRT ) fprintf(ioQQQ,"DEBUG pops ");
3021  while( nEner < nLevels_per_elec[0] && sum_pop/hmi.H2_total < FRAC )
3022  {
3023 
3024  /* array of energy sorted indices within X */
3025  ip = H2_ipX_ener_sort[nEner];
3026  iVib = ipVib_H2_energy_sort[ip];
3027  iRot = ipRot_H2_energy_sort[ip];
3028  sum_pop += H2_old_populations[iElec][iVib][iRot];
3029  if( PRT ) fprintf(ioQQQ,"\t%.3e ", H2_old_populations[iElec][iVib][iRot]);
3030  ++nEner;
3031  }
3032  if( PRT ) fprintf(ioQQQ,"\n ");
3034  /*fprintf(ioQQQ,"DEBUG zone %.2f old nmatrix %li proposed nmatrix %li sum_pop %.4e H2_total %.4e\n",
3035  fnzone , nXLevelsMatrix ,nEner , sum_pop,hmi.H2_total);
3036  nXLevelsMatrix = nEner;*/
3037 # undef FRAC
3038 # undef PRT
3039  }
3040  return;
3041 }
3042 /*lint -e802 possible bad pointer */
3043 
3044 /*H2_cooling evaluate cooling and heating due to H2 molecule, called by
3045  * H2_LevelPops in convergence loop when h2 heating is important, also
3046  * called by CoolEvaluate to get final heating - argument is name of
3047  * routine that called it */
3048 #if defined(__ICC) && defined(__i386)
3049 #pragma optimization_level 1
3050 #endif
3052  /* string saying who called this routine,
3053  * "H2lup" call within H2 level populations solver
3054  * "CoolEvaluate" call from main cooling routine */
3055  const char *chRoutine)
3056 {
3057  long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
3058  double heatone ,
3059  rate_dn_heat,
3060  rate_up_cool;
3061  long int nColl,
3062  ipHi, ipLo;
3063  double Big1_heat , Big1_cool,
3064  H2_X_col_cool , H2_X_col_heat;
3065  long int ipVib_big_heat_hi,ipVib_big_heat_lo ,ipRot_big_heat_hi ,
3066  ipRot_big_heat_lo ,ipVib_big_cool_hi,ipVib_big_cool_lo ,
3067  ipRot_big_cool_hi , ipRot_big_cool_lo;
3068 /*# define DEBUG_DIS_DEAT*/
3069 # ifdef DEBUG_DIS_DEAT
3070  double heatbig;
3071  long int iElecBig , iVibBig , iRotBig;
3072 # endif
3073 
3074  /* option to keep track of strongest single heating agent due to collisions
3075  * within X */
3076  enum {DEBUG_H2_COLL_X_HEAT=false };
3077 
3078  DEBUG_ENTRY( "H2_Cooling()" );
3079 
3080  /* possible debug counters, counter itself, counter to turn on
3081  * debug output, and counter for stopping */
3082  static long int nCount=-1;
3083  long int nCountDebug = 930,
3084  nCountStop = 940;
3085  ++nCount;
3086 
3087  /* nCallH2_this_iteration is not incremented until after the level
3088  * populations have converged the first time. so for the first n calls
3089  * this will return zero, a good idea since populations will be wildly
3090  * incorrect during search for first valid pops */
3091  if( !h2.lgH2ON || !nCallH2_this_iteration )
3092  {
3093  hmi.HeatH2Dexc_BigH2 = 0.;
3094  hmi.HeatH2Dish_BigH2 = 0.;
3096  return;
3097  }
3098 
3099  hmi.HeatH2Dish_BigH2 = 0.;
3100  heatone = 0.;
3101 # ifdef DEBUG_DIS_DEAT
3102  heatbig = 0.;
3103  iElecBig = -1;
3104  iVibBig = -1;
3105  iRotBig = -1;
3106 # endif
3107  /* heating due to dissociation of electronic excited states */
3108  for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
3109  {
3110  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
3111  {
3112  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
3113  {
3114  /* population, cm-3, in excited state */
3115  heatone = H2_populations[iElecHi][iVibHi][iRotHi] *
3116  H2_dissprob[iElecHi][iVibHi][iRotHi] *
3117  H2_disske[iElecHi][iVibHi][iRotHi];
3118  hmi.HeatH2Dish_BigH2 += heatone;
3119 # ifdef DEBUG_DIS_DEAT
3120  if( heatone > heatbig )
3121  {
3122  heatbig = heatone;
3123  iElecBig = iElecHi;
3124  iVibBig = iVibHi;
3125  iRotBig = iRotHi;
3126  }
3127 # endif
3128  }
3129  }
3130  }
3131 # ifdef DEBUG_DIS_DEAT
3132  fprintf(ioQQQ,"DEBUG H2 dis heat\t%.2f\t%.3f\t%li\t\t%li\t%li\n",
3133  fnzone ,
3134  heatbig / SDIV( hmi.HeatH2Dish_BigH2 ) ,
3135  iElecBig ,
3136  iVibBig ,
3137  iRotBig );
3138 # endif
3139  /* dissociation heating HeatH2Dish_BigH2 was in eV -
3140  * convert to ergs */
3142 
3143  /*fprintf(ioQQQ,"DEBUG H2 heat/dissoc %.3e\n",
3144  hmi.HeatH2Dish_BigH2/ SDIV(hmi.H2_rate_destroy*hmi.H2_total) );*/
3145  /* now work on collisional heating due to bound-bound
3146  * collisional transitions within X */
3147  hmi.HeatH2Dexc_BigH2 = 0.;
3148  H2_X_col_cool = 0.;
3149  H2_X_col_heat = 0.;
3150  /* these are the colliders that will be considered as depopulating agents */
3151  /* the colliders are H, He, H2 ortho, H2 para, H+ */
3152  /* atomic hydrogen */
3153  long int nBug1 = 200 , nBug2 = 201;
3154 # if 0
3155  /* fudge(-1) returns number of fudge factors entered on fudge command */
3156  if( fudge(-1) > 0 )
3157  {
3158  nBug1 = (long)fudge(0);
3159  nBug2 = (long)fudge(1);
3160  }
3161 # endif
3162  if( DEBUG_H2_COLL_X_HEAT && (nCount == nBug1 || nCount==nBug2) )
3163  {
3164  FILE *ioBAD=NULL;
3165  if( nCount==nBug1 )
3166  {
3167  ioBAD = fopen("firstpop.txt" , "w" );
3168  }
3169  else if( nCount==nBug2 )
3170  {
3171  ioBAD = fopen("secondpop.txt" , "w" );
3172  }
3173  for( ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
3174  {
3175  long int ip = H2_ipX_ener_sort[ipHi];
3176  iVibHi = ipVib_H2_energy_sort[ip];
3177  iRotHi = ipRot_H2_energy_sort[ip];
3178  fprintf(ioBAD , "%li\t%li\t%.2e\n",
3179  iVibHi , iRotHi ,
3180  H2_populations[0][iVibHi][iRotHi] );
3181  }
3182  fclose(ioBAD);
3183  }
3184 
3185  /* now make sum of all collisions within X itself */
3186  iElecHi = 0;
3187  iElecLo = 0;
3188  Big1_heat = 0.;
3189  Big1_cool = 0.;
3190  ipVib_big_heat_hi = -1;
3191  ipVib_big_heat_lo = -1;
3192  ipRot_big_heat_hi = -1;
3193  ipRot_big_heat_lo = -1;
3194  ipVib_big_cool_hi = -1;
3195  ipVib_big_cool_lo = -1;
3196  ipRot_big_cool_hi = -1;
3197  ipRot_big_cool_lo = -1;
3198  /* this will be derivative */
3200  for( ipHi=1; ipHi<nLevels_per_elec[iElecHi]; ++ipHi )
3201  {
3202  long int ip = H2_ipX_ener_sort[ipHi];
3203  iVibHi = ipVib_H2_energy_sort[ip];
3204  iRotHi = ipRot_H2_energy_sort[ip];
3205  if( iVibHi > VIB_COLLID )
3206  continue;
3207 
3208  realnum H2statHi = H2_stat[iElecHi][iVibHi][iRotHi];
3209  double H2boltzHi = H2_Boltzmann[iElecHi][iVibHi][iRotHi];
3210  double H2popHi = H2_populations[iElecHi][iVibHi][iRotHi];
3211  double ewnHi = energy_wn[iElecHi][iVibHi][iRotHi];
3212 
3213  for( ipLo=0; ipLo<ipHi; ++ipLo )
3214  {
3215  double coolone , oneline;
3216  ip = H2_ipX_ener_sort[ipLo];
3217  iVibLo = ipVib_H2_energy_sort[ip];
3218  iRotLo = ipRot_H2_energy_sort[ip];
3219  if( iVibLo > VIB_COLLID)
3220  continue;
3221 
3222  rate_dn_heat = 0.;
3223 
3224  /* this sum is total downward heating summed over all colliders */
3225  mr5ci H2cr = H2_CollRate.begin(iVibHi,iRotHi,iVibLo,iRotLo);
3226  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
3227  /* downward collision rate */
3228  rate_dn_heat += H2cr[nColl]*collider_density[nColl];
3229 
3230  /* now get upward collisional cooling by detailed balance */
3231  rate_up_cool = rate_dn_heat * H2_populations[iElecLo][iVibLo][iRotLo] *
3232  /* rest converts into upward collision rate */
3233  H2statHi / H2_stat[iElecLo][iVibLo][iRotLo] *
3234  H2boltzHi / SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] );
3235 
3236  rate_dn_heat *= H2popHi;
3237 
3238  /* net heating due to collisions within X -
3239  * positive if heating, negative is cooling
3240  * this will usually be heating if X is photo pumped
3241  * in printout and in punch heating this is called "H2cX" */
3242  double conversion = (ewnHi - energy_wn[iElecLo][iVibLo][iRotLo]) * ERG1CM;
3243  heatone = rate_dn_heat * conversion;
3244  coolone = rate_up_cool * conversion;
3245  /* this is net heating, negative if cooling */
3246  oneline = heatone - coolone;
3247  hmi.HeatH2Dexc_BigH2 += oneline;
3248 
3249  /* keep track of heating and cooling separately */
3250  H2_X_col_cool += coolone;
3251  H2_X_col_heat += heatone;
3252 
3253  if( 0 && DEBUG_H2_COLL_X_HEAT && (nCount == 692 || nCount==693) )
3254  {
3255  static FILE *ioBAD=NULL;
3256  if(ipHi == 1 && ipLo == 0)
3257  {
3258  if( nCount==692 )
3259  {
3260  ioBAD = fopen("firstheat.txt" , "w" );
3261  }
3262  else if( nCount==693 )
3263  {
3264  ioBAD = fopen("secondheat.txt" , "w" );
3265  }
3266  fprintf(ioBAD,"DEBUG start \n");
3267  }
3268  fprintf(ioBAD,"DEBUG BAD DAY %li %li %li %li %.3e %.3e\n",
3269  iVibHi , iRotHi, iVibLo , iRotLo , heatone , coolone );
3270  if( ipHi==nLevels_per_elec[iElecHi]-1 &&
3271  ipLo==ipHi-1 )
3272  fclose(ioBAD);
3273  }
3274 
3275  /* derivative wrt temperature - assume exp wrt ground -
3276  * this needs to be divided by square of temperature in wn -
3277  * done at end of loop */
3278  hmi.deriv_HeatH2Dexc_BigH2 += (realnum)(oneline * ewnHi);
3279 
3280  /* oneline is net heating, positive for heating, negative
3281  * for cooling */
3282  if( DEBUG_H2_COLL_X_HEAT )
3283  {
3284  if( oneline < Big1_cool )
3285  {
3286  Big1_cool = oneline;
3287  ipVib_big_cool_hi = iVibHi;
3288  ipVib_big_cool_lo = iVibLo;
3289  ipRot_big_cool_hi = iRotHi;
3290  ipRot_big_cool_lo = iRotLo;
3291  }
3292  else if( oneline > Big1_heat )
3293  {
3294  Big1_heat = oneline;
3295  ipVib_big_heat_hi = iVibHi;
3296  ipVib_big_heat_lo = iVibLo;
3297  ipRot_big_heat_hi = iRotHi;
3298  ipRot_big_heat_lo = iRotLo;
3299  }
3300  }
3301 
3302  /* this would be a major logical error */
3303  ASSERT(
3304  (rate_up_cool==0 && rate_dn_heat==0) ||
3305  (energy_wn[iElecHi][iVibHi][iRotHi] > energy_wn[iElecLo][iVibLo][iRotLo]) );
3306  }/* end loop over lower levels, all collisions within X */
3307  }/* end loop over upper levels, all collisions within X */
3308 
3309  /* this debug statement will identify the single strongest heating or
3310  * cooling agent within X and give the lowest 7 level populations */
3311  if( DEBUG_H2_COLL_X_HEAT )
3312  {
3313  /* nCount counts number of calls through here,
3314  * nCountDebug is count to turn on debug prints,
3315  * nCountStop is count to abort code */
3316  if(nCount>nCountDebug )
3317  {
3318  fprintf(ioQQQ,
3319  "DEBUG H2_Cooling A %li %15s, Te %.3e net Heat(Xcol) %.2e "
3320  "heat %.2e cool %.2e H+/0 "
3321  "%.2e n(H2)%.3e Sol rat %.3e grn J1->0%.2e frac heat 1 line "
3322  "%.2e Hi(v,j)%li %li Lo(v,J)%li %li frac cool 1 line %.2e "
3323  "Hi(v,j)%li %li Lo(v,J)%li %li POP(J=1,13)",
3324  nCount , chRoutine,
3325  phycon.te,
3327  H2_X_col_cool ,
3328  H2_X_col_heat ,
3330  hmi.H2_total,
3333  Big1_heat/hmi.HeatH2Dexc_BigH2 ,
3334  ipVib_big_heat_hi , ipRot_big_heat_hi , ipVib_big_heat_lo , ipRot_big_heat_lo ,
3335  Big1_cool/hmi.HeatH2Dexc_BigH2 ,
3336  ipVib_big_cool_hi , ipRot_big_cool_hi , ipVib_big_cool_lo , ipRot_big_cool_lo );
3337 
3338  for( iRotLo=0; iRotLo<14; ++iRotLo )
3339  {
3340  fprintf(ioQQQ,"\t%.2e" ,
3341  H2_populations[0][0][iRotLo]/hmi.H2_total );
3342  }
3343  fprintf(ioQQQ,"\t%li\n",nCount);
3344  /* now give collision rates for strongest heat/cool level */
3345  fprintf(ioQQQ,"DEBUG H2_Cooling B heat Coll Rate (lrg col) dn,up" );
3346  double HeatNet = 0.;
3347  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
3348  {
3349  fprintf(ioQQQ,"\t%.2e" ,
3350  H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3351  collider_density[nColl] );
3352  fprintf(ioQQQ,"\t%.2e" ,
3353  /* downward collision rate */
3354  H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3355  collider_density[nColl]*
3356  /* rest converts into upward collision rate */
3357  H2_stat[0][ipVib_big_heat_hi][ipRot_big_heat_hi] / H2_stat[0][ipVib_big_heat_lo][ipRot_big_heat_lo] *
3358  H2_Boltzmann[0][ipVib_big_heat_hi][ipRot_big_heat_hi] /
3359  SDIV( H2_Boltzmann[0][ipVib_big_heat_lo][ipRot_big_heat_lo] ) );
3360  HeatNet +=
3361  H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3362  collider_density[nColl]*H2_populations[iElecLo][ipRot_big_heat_hi][ipVib_big_heat_lo]
3363  -
3364  H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3365  collider_density[nColl]*
3366  /* rest converts into upward collision rate */
3367  H2_stat[0][ipVib_big_heat_hi][ipRot_big_heat_hi] / H2_stat[0][ipVib_big_heat_lo][ipRot_big_heat_lo] *
3368  H2_Boltzmann[0][ipVib_big_heat_hi][ipRot_big_heat_hi] /
3369  SDIV( H2_Boltzmann[0][ipVib_big_heat_lo][ipRot_big_heat_lo] ) *
3370  H2_populations[iElecLo][ipVib_big_heat_lo][ipRot_big_heat_lo] ;
3371  }
3372  /* HeatNet includes the level populations, rates do not */
3373  fprintf( ioQQQ , " HeatNet %.2e",HeatNet);
3374  fprintf(ioQQQ,"\n");
3375  fprintf(ioQQQ,"DEBUG H2_Cooling C cool Coll Rate (lrg col) dn,up" );
3376  double CoolNet = 0.;
3377  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
3378  {
3379  fprintf(ioQQQ,"\t%.2e" ,
3380  H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3381  collider_density[nColl] );
3382  fprintf(ioQQQ,"\t%.2e" ,
3383  /* downward collision rate */
3384  H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3385  collider_density[nColl]*
3386  /* rest converts into upward collision rate */
3387  H2_stat[0][ipVib_big_cool_hi][ipRot_big_cool_hi] / H2_stat[0][ipVib_big_cool_lo][ipRot_big_cool_lo] *
3388  H2_Boltzmann[0][ipVib_big_cool_hi][ipRot_big_cool_hi] /
3389  SDIV( H2_Boltzmann[0][ipVib_big_cool_lo][ipRot_big_cool_lo] ) );
3390  CoolNet +=
3391  H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3392  collider_density[nColl]*H2_populations[iElecLo][ipRot_big_cool_hi][ipVib_big_cool_lo]
3393  -
3394  H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3395  collider_density[nColl]*
3396  /* rest converts into upward collision rate */
3397  H2_stat[0][ipVib_big_cool_hi][ipRot_big_cool_hi] / H2_stat[0][ipVib_big_cool_lo][ipRot_big_cool_lo] *
3398  H2_Boltzmann[0][ipVib_big_cool_hi][ipRot_big_cool_hi] /
3399  SDIV( H2_Boltzmann[0][ipVib_big_cool_lo][ipRot_big_cool_lo] ) *
3400  H2_populations[iElecLo][ipVib_big_cool_lo][ipRot_big_cool_lo] ;
3401  }
3402  /* CoolNet includes the level populations, rates do not */
3403  fprintf( ioQQQ , " CoolNet %.2e",CoolNet);
3404  fprintf(ioQQQ,"\n");
3405  }
3406  }
3407 
3408  /* this is inside h2 cooling, and is called extra times when H2 heating is important */
3409  if( PRT_POPS )
3410  fprintf(ioQQQ,
3411  " DEBUG H2 heat fnzone\t%.2f\trenrom\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n",
3412  fnzone ,
3414  phycon.te ,
3417 
3418  /* this is derivative of collisional heating wrt temperature - needs
3419  * to be divided by square of temperature in wn */
3421 
3422  {
3423  enum {DEBUG_LOC=false };
3424  if( DEBUG_H2_COLL_X_HEAT && DEBUG_LOC &&
3425  (fabs(hmi.HeatH2Dexc_BigH2) > SMALLFLOAT) )
3426  {
3427  int iVib = 0;
3428 
3429  /*fprintf(ioQQQ," H2_cooling pops\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
3430  H2_populations[0][iVib][0]/hmi.H2_total,
3431  H2_populations[0][iVib][1]/hmi.H2_total,
3432  H2_populations[0][iVib][2]/hmi.H2_total,
3433  H2_populations[0][iVib][3]/hmi.H2_total,
3434  H2_populations[0][iVib][4]/hmi.H2_total,
3435  H2_populations[0][iVib][5]/hmi.H2_total);*/
3436 
3437  iElecHi = iElecLo = 0;
3438  iVibHi = iVibLo = 0;
3439  iRotHi = 7;
3440  iRotLo = 5;
3441  rate_dn_heat = rate_up_cool = 0.;
3442  /* this sum is total downward heating */
3443  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
3444  {
3445  rate_dn_heat +=
3446  H2_populations[iElecHi][iVibHi][iRotHi] *
3447  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl]*
3448  collider_density[nColl];
3449 
3450  /* now get upward collisional cooling by detailed balance */
3451  rate_up_cool +=
3452  H2_populations[iElecLo][iVibLo][iRotLo] *
3453  /* downward collision rate */
3454  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl]*
3455  collider_density[nColl]*
3456  /* rest converts into upward collision rate */
3457  H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElecLo][iVibLo][iRotLo] *
3458  H2_Boltzmann[iElecHi][iVibHi][iRotHi] /
3459  SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] );
3460  }
3461 
3462  fprintf(ioQQQ,"DEBUG H2_cooling D pop %li ov %li\t%.3e\tdn up 31\t%.3e\t%.3e\n",
3463  iRotHi , iRotLo ,
3464  H2_populations[0][iVib][iRotHi]/H2_populations[0][iVib][iRotLo],
3465  rate_dn_heat,
3466  rate_up_cool);
3467  if( nCount>= nCountStop )
3468  cdEXIT(EXIT_FAILURE);
3469  }
3470  }
3471  {
3472  enum {DEBUG_LOC=false};
3473  if( DEBUG_LOC )
3474  {
3475  static long nzdone=-1 , nzincre;
3476  if( nzone!=nzdone )
3477  {
3478  nzdone = nzone;
3479  nzincre = -1;
3480  }
3481  ++nzincre;
3482  fprintf(ioQQQ," H2 nz\t%.2f\tnzinc\t%li\tTe\t%.4e\tH2\t%.3e\tcXH\t%.2e\tdcXH/dt%.2e\tDish\t%.2e \n",
3483  fnzone,
3484  nzincre,
3485  phycon.te,
3486  hmi.H2_total ,
3490 
3491  }
3492  }
3493 
3494 # if 0
3495  /* this can be noisy due to finite accuracy of solution, so take average with
3496  * previous value */
3497  /*>>chng 04 mar 01, do not take average */
3498  if( 1 || nzone <1 || old_HeatH2Dexc==0. || nCallH2_this_iteration <2)
3499  {
3500  old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;
3501  }
3502  else
3503  {
3504  hmi.HeatH2Dexc_BigH2 = (hmi.HeatH2Dexc_BigH2+old_HeatH2Dexc)/2.f;
3505  old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2;
3506  }
3507 # endif
3508  {
3509  enum {DEBUG_LOC=false};
3510  if( DEBUG_LOC /*&& DEBUG_H2_COLL_X_HEAT*/ )
3511  {
3512  fprintf(ioQQQ,"DEBUG H2_cooling E %15s %c vib deex %li Te %.3e net heat %.3e cool %.3e heat %.3e\n",
3513  chRoutine ,
3514  TorF(conv.lgSearch),
3515  nCount,
3517  H2_X_col_cool ,
3518  H2_X_col_heat /*,
3519  H2_populations[0][0][7]/SDIV(H2_populations[0][0][5]) ,
3520  H2_populations[0][0][13]/SDIV(H2_populations[0][0][11])*/ );
3521  if( 0 && nCount > nCountStop )
3522  {
3523  cdEXIT( EXIT_FAILURE );
3524  }
3525  }
3526  }
3527  if( mole.nH2_TRACE >= mole.nH2_trace_full )
3528  fprintf(ioQQQ,
3529  " H2_Cooling Ctot\t%.4e\t HeatH2Dish_BigH2 \t%.4e\t HeatH2Dexc_BigH2 \t%.4e\n" ,
3530  thermal.ctot ,
3533 
3534  /* when we are very far from solution, during search phase, collisions within
3535  * X can be overwhelmingly large heating and cooling terms, which nearly
3536  * cancel out. Some dense cosmic ray heated clouds could not find correct
3537  * initial solution due to noise introduced by large net heating which was
3538  * the very noisy tiny difference between very large heating and cooling
3539  * terms. Do not include collisions with x as heat/cool during the
3540  * initial search phase */
3541  if( conv.lgSearch )
3542  hmi.HeatH2Dexc_BigH2 = 0.;
3543  return;
3544 }
3545 
3546 
3547 /*cdH2_colden return column density in H2, negative -1 if cannot find state,
3548  * header is cdDrive */
3549 double cdH2_colden( long iVib , long iRot )
3550 {
3551 
3552  /*if iVib is negative, return
3553  * total column density - iRot=0
3554  * ortho column density - iRot 1
3555  * para column density - iRot 2
3556  * else return column density in iVib, iRot */
3557  if( iVib < 0 )
3558  {
3559  if( iRot==0 )
3560  {
3561  /* return total H2 column density */
3562  return( h2.ortho_colden + h2.para_colden );
3563  }
3564  else if( iRot==1 )
3565  {
3566  /* return ortho H2 column density */
3567  return h2.ortho_colden;
3568  }
3569  else if( iRot==2 )
3570  {
3571  /* return para H2 column density */
3572  return h2.para_colden;
3573  }
3574  else
3575  {
3576  fprintf(ioQQQ," iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n");
3577  return -1.;
3578  }
3579  }
3580  else if( h2.lgH2ON )
3581  {
3582  /* this branch want state specific column density, which can only result from
3583  * evaluation of big molecule */
3584  int iElec = 0;
3585  if( iRot <0 || iVib >h2.nVib_hi[iElec] || iRot > h2.nRot_hi[iElec][iVib])
3586  {
3587  fprintf(ioQQQ," iVib and iRot must lie within X, returning -2.\n");
3588  fprintf(ioQQQ," iVib must be <= %li and iRot must be <= %li.\n",
3589  h2.nVib_hi[iElec],h2.nRot_hi[iElec][iVib]);
3590  return -2.;
3591  }
3592  else
3593  {
3594  return H2_X_colden[iVib][iRot];
3595  }
3596  }
3597  /* error condition - no valid parameter */
3598  else
3599  return -1;
3600 }
3601 
3602 /*H2_Colden maintain H2 column densities within X */
3603 void H2_Colden( const char *chLabel )
3604 {
3605  long int iVib , iRot;
3606 
3607  /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */
3608  if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ )
3609  return;
3610 
3611  DEBUG_ENTRY( "H2_Colden()" );
3612 
3613  if( strcmp(chLabel,"ZERO") == 0 )
3614  {
3615  /* zero out formation rates and column densites */
3616  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
3617  {
3618  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
3619  {
3620  /* space for the rotation quantum number */
3621  H2_X_colden[iVib][iRot] = 0.;
3622  H2_X_colden_LTE[iVib][iRot] = 0.;
3623  }
3624  }
3625  }
3626 
3627  else if( strcmp(chLabel,"ADD ") == 0 )
3628  {
3629  /* add together column densities */
3630  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
3631  {
3632  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
3633  {
3634  /* state specific H2 column density */
3635  H2_X_colden[iVib][iRot] += (realnum)(H2_populations[0][iVib][iRot]*radius.drad_x_fillfac);
3636  /* LTE state specific H2 column density - H2_populations_LTE is normed to unity
3637  * so must be multiplied by total H2 density */
3638  H2_X_colden_LTE[iVib][iRot] += (realnum)(H2_populations_LTE[0][iVib][iRot]*
3640  }
3641  }
3642  }
3643 
3644  /* we will not print column densities so skip that - if not print then we have a problem */
3645  else if( strcmp(chLabel,"PRIN") != 0 )
3646  {
3647  fprintf( ioQQQ, " H2_Colden does not understand the label %s\n",
3648  chLabel );
3649  cdEXIT(EXIT_FAILURE);
3650  }
3651 
3652  return;
3653 }
3654 
3655 /*H2_DR choose next zone thickness based on H2 big molecule */
3656 double H2_DR(void)
3657 {
3658  return BIGFLOAT;
3659 }
3660 
3661 /*H2_RT_OTS - add H2 ots fields */
3662 void H2_RT_OTS( void )
3663 {
3664 
3665  long int iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo;
3666 
3667  /* do not compute if H2 not turned on, or not computed for these conditions */
3668  if( !h2.lgH2ON || !h2.nCallH2_this_zone )
3669  return;
3670 
3671  DEBUG_ENTRY( "H2_RT_OTS()" );
3672 
3673  /* loop over all possible lines and set H2_populations, and quantities that depend on escape prob, dest, etc */
3674  long int lim_elec_hi = 0;
3675  for( iElecHi=0; iElecHi<=lim_elec_hi; ++iElecHi )
3676  {
3677  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
3678  {
3679  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
3680  {
3681  double H2popHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop;
3682  long int lim_elec_lo = 0;
3683  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
3684  {
3685  /* want to include all vibration states in lower level if different electronic level,
3686  * but only lower vibration levels if same electronic level */
3687  long int nv = h2.nVib_hi[iElecLo];
3688  if( iElecLo==iElecHi )
3689  nv = iVibHi;
3690  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
3691  {
3692  long nr = h2.nRot_hi[iElecLo][iVibLo];
3693  if( iElecLo==iElecHi && iVibHi==iVibLo )
3694  nr = iRotHi-1;
3695 
3696  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
3697  {
3698  /* >>chng 05 feb 07, use lgH2_line_exists */
3699  if( iElecHi==0 && lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
3700  {
3701  /* ots destruction rate */
3702  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots =
3703  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul * H2popHi *
3704  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Pdest;
3705 
3706  /* dump the ots rate into the stack - but only for ground electronic state*/
3708  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots,
3709  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont );
3710  }
3711  }
3712  }
3713  }
3714  }
3715  }
3716  }
3717 
3718  return;
3719 }
3720 /*lint +e802 possible bad pointer */
3721 

Generated for cloudy by doxygen 1.8.1.1