cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /* mole_H2_LTE sets Boltzmann factors and LTE unit population of large H2 molecular */ 00004 /* H2_init - called to initialize things from cdInit */ 00005 /*H2_init_coreload one time initialization */ 00006 /* H2_Zero zero out vars in the large H2 molecule, called from zero 00007 * before any commands are parsed */ 00008 /* H2_zero_pops_too_low - zero out some H2 variables if we decide not to compute 00009 * the full sim, called by H2_LevelPops*/ 00010 /* H2_Solomon_rate find rates between H2s and H2g and other levels, 00011 * for eventual use in the chemistry */ 00012 /* H2_gs_rates evaluate rates between ground and star states of H2 for use in chemistry */ 00013 /* H2_He_coll_init initialize H2 - He collision data set 00014 * H2_He_coll interpolate on h2 - He collision data set to return rate at temp*/ 00015 #include "cddefines.h" 00016 #include "phycon.h" 00017 #include "mole.h" 00018 #include "hmi.h" 00019 #include "taulines.h" 00020 #include "h2.h" 00021 #include "h2_priv.h" 00022 00023 /*H2_Solomon_rate find rates between H2s and H2g and other levels, 00024 * for eventual use in the chemistry */ 00025 void H2_Solomon_rate( void ) 00026 { 00027 00028 long int iElecLo , iElecHi , iVibLo , iVibHi , iRotLo , iRotHi; 00029 00030 DEBUG_ENTRY( "H2_Solomon_rate()" ); 00031 00032 /* iElecLo will always be X in this routine */ 00033 iElecLo = 0; 00034 00035 /* find rate (s-1) h2 dissociation into X continuum by Solomon process and 00036 * assign to the TH85 g and s states 00037 * these will go back into the chemistry network */ 00038 00039 /* rates [s-1] for dissociation from s or g, into electronic excited states 00040 * followed by dissociation */ 00041 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.; 00042 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.; 00043 00044 /* these are used in a print statement - are they needed? */ 00045 hmi.H2_Solomon_elec_decay_H2g = 0.; 00046 hmi.H2_Solomon_elec_decay_H2s = 0.; 00047 00048 /* at this point we have already evaluated the sum of the radiative rates out 00049 * of the electronic excited states - this is H2_rad_rate_out[electronic][vib][rot] 00050 * and this includes both decays into the continuum and bound states of X */ 00051 00052 /* sum over all electronic states, finding dissociation rate */ 00053 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00054 { 00055 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00056 { 00057 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00058 { 00059 double factor = (double)H2_dissprob[iElecHi][iVibHi][iRotHi]/ 00060 H2_rad_rate_out[iElecHi][iVibHi][iRotHi]; 00061 double H2popHi = H2_populations[iElecHi][iVibHi][iRotHi]; 00062 00063 /* loop over all levels within X to find 00064 * decay rates from H2g and H2s to continuum 00065 * distinction between H2g and H2s is determined 00066 * by ENERGY_H2_STAR */ 00067 iElecLo = 0; 00068 00069 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo ) 00070 { 00071 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00072 00073 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00074 mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00075 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00076 { 00077 if( *lgH2le++ ) 00078 { 00079 /* this is the rate [cm-3 s-1] that mole goes from 00080 * lower level into electronic excited states then 00081 * into continuum */ 00082 double rate_up_cont = 00083 H2_populations[iElecLo][iVibLo][iRotLo]* 00084 H2L->Emis->pump*factor; 00085 00086 /* rate electronic state decays into H2g */ 00087 double elec_decay = 00088 H2popHi* 00089 H2L->Emis->Aul*(H2L->Emis->Pesc+H2L->Emis->Pdest+H2L->Emis->Pelec_esc); 00090 00091 if( energy_wn[0][iVibLo][iRotLo] > ENERGY_H2_STAR ) 00092 { 00093 /* this is H2g up to excited then to continuum - 00094 * cm-3 s-1 at this point */ 00095 hmi.H2_Solomon_dissoc_rate_BigH2_H2s += rate_up_cont; 00096 /* rate electronic state decays into H2g */ 00097 hmi.H2_Solomon_elec_decay_H2s += elec_decay; 00098 } 00099 else 00100 { 00101 /* this is H2g up to excited then to continuum - 00102 * cm-3 s-1 at this point */ 00103 hmi.H2_Solomon_dissoc_rate_BigH2_H2g += rate_up_cont; 00104 /* rate electronic state decays into H2g */ 00105 hmi.H2_Solomon_elec_decay_H2g += elec_decay; 00106 } 00107 } 00108 ++H2L; 00109 } 00110 } 00111 }/* end iRotHi */ 00112 }/* end iVibHi */ 00113 }/* end iElecHi */ 00114 /* at this point units of hmi.H2_Solomon_elec_decay_H2g, H2s are cm-3 s-1 00115 * since H2_populations are included - 00116 * div by pops to get actual dissocation rate, s-1 */ 00117 if( hmi.H2_total > SMALLFLOAT ) 00118 { 00119 hmi.H2_Solomon_elec_decay_H2g /= SDIV( H2_sum_excit_elec_den ); 00120 hmi.H2_Solomon_elec_decay_H2s /= SDIV( H2_sum_excit_elec_den ); 00121 00122 /* will be used for H2s-> H + H */ 00123 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_H2s / SDIV(H2_den_s); 00124 00125 /* will be used for H2g-> H + H */ 00126 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_H2g / SDIV(H2_den_g); 00127 00128 } 00129 else 00130 { 00131 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0; 00132 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0; 00133 00134 } 00135 /*fprintf(ioQQQ,"DEBUG H2 new %.2e %.2e %.2e %.2e \n", 00136 hmi.H2_Solomon_elec_decay_H2g , 00137 hmi.H2_Solomon_elec_decay_H2s , 00138 hmi.H2_Solomon_dissoc_rate_BigH2_H2s, 00139 hmi.H2_Solomon_dissoc_rate_BigH2_H2g );*/ 00140 00141 /* rate g goes to s */ 00142 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.; 00143 return; 00144 } 00145 00146 /*H2_gs_rates evaluate rates between ground and star states of H2 for use in chemistry */ 00147 void H2_gs_rates( void ) 00148 { 00149 long int ipLoX , ip , iRotLoX , iVibLoX , 00150 iElecHi , iVibHi , ipOther , iRotOther, 00151 iRotHi , iVibOther; 00152 00153 DEBUG_ENTRY( "H2_gs_rates()" ); 00154 00155 /* rate g goes to s */ 00156 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.; 00157 00158 /* loop over all levels in H2g */ 00159 for( ipLoX=0; ipLoX < nEner_H2_ground; ++ipLoX ) 00160 { 00161 ip = H2_ipX_ener_sort[ipLoX]; 00162 iRotLoX = ipRot_H2_energy_sort[ip]; 00163 iVibLoX = ipVib_H2_energy_sort[ip]; 00164 /* now find all pumps up to electronic excited states */ 00165 /* sum over all electronic states, finding dissociation rate */ 00166 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00167 { 00168 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00169 { 00170 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00171 { 00172 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX] ) 00173 { 00174 /* this is the rate {cm-3 s-1] that mole goes from 00175 * lower level into electronic excited states then 00176 * into continuum */ 00177 double rate_up_cont = 00178 H2_populations[0][iVibLoX][iRotLoX]* 00179 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX].Emis->pump; 00180 00181 double decay_star = H2_rad_rate_out[iElecHi][iVibHi][iRotHi] - H2_dissprob[iElecHi][iVibHi][iRotHi]; 00182 /* loop over all other levels in H2g, subtracting 00183 * their rate - remainder is rate into star, this is 00184 * usually only a few levels */ 00185 for( ipOther=0; ipOther < nEner_H2_ground; ++ipOther ) 00186 { 00187 ip = H2_ipX_ener_sort[ipOther]; 00188 iRotOther = ipRot_H2_energy_sort[ip]; 00189 iVibOther = ipVib_H2_energy_sort[ip]; 00190 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther] ) 00191 { 00192 decay_star -= 00193 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Aul*( 00194 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Pesc+ 00195 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Pdest+ 00196 H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Emis->Pelec_esc); 00197 } 00198 } 00199 /* MAX because may underflow to negative numbers is rates very large 00200 * this is fraction that returns to H2s */ 00201 decay_star = MAX2(0., decay_star)/SDIV(H2_rad_rate_out[iElecHi][iVibHi][iRotHi]); 00202 hmi.H2_H2g_to_H2s_rate_BigH2 += rate_up_cont*decay_star; 00203 00204 }/* end if line exists */ 00205 }/* end loop rot electronic excited */ 00206 }/* end loop vib electronic excited */ 00207 }/* end loop electronic electronic excited */ 00208 } 00209 00210 /* at this point units are cm-3 s-1 - convert to rate s-1 */ 00211 hmi.H2_H2g_to_H2s_rate_BigH2 /= SDIV( H2_den_g ); 00212 return; 00213 } 00214 00215 /* H2_zero_pops_too_low - zero out some H2 variables if we decide not to compute 00216 * the full sim, called by H2_LevelPops*/ 00217 void H2_zero_pops_too_low( void ) 00218 { 00219 00220 long int iElec, iElecHi, iElecLo, iVib , iVibLo , iVibHi , 00221 iRot , iRotHi , iRotLo; 00222 00223 DEBUG_ENTRY( "H2_zero_pops_too_low()" ); 00224 00225 /* >>chng 05 jan 26, add this block to set populations to LTE value */ 00226 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 00227 { 00228 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 00229 { 00230 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 00231 { 00232 /* LTE populations are for unit H2 density, so need to multiply 00233 * by total H2 density */ 00234 H2_old_populations[iElec][iVib][iRot] = 00235 (realnum)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total; 00236 H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot]; 00237 } 00238 } 00239 } 00240 /* zero everything out - loop over all possible lines */ 00241 /* >>chng 05 jan 26, set to LTE values, since we still need to accumulate Lyman line 00242 * optical depths to have correct self-shielding when large h2 does come on */ 00243 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00244 { 00245 pops_per_elec[iElecHi] = 0.; 00246 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00247 { 00248 pops_per_vib[iElecHi][iVibHi] = 0.; 00249 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00250 { 00251 long int lim_elec_lo = 0; 00252 // this will update Lo->Pop as well as Lo and Hi point to the same set of data structures 00253 // the loops below are OK since Lo->Pop is only accessed after Hi->Pop has been set. 00254 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop = H2_populations[iElecHi][iVibHi][iRotHi]; 00255 /* now the lower levels 00256 * NB - iElecLo the lower electronic level is only X 00257 * we don't consider excited electronic to excited electronic trans here, since we are only 00258 * concerned with excited electronic levels as a photodissociation process 00259 * code exists to relax this assumption - simply change following to iElecHi */ 00260 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00261 { 00262 /* want to include all vib states in lower level if different electronic level, 00263 * but only lower vib levels if same electronic level */ 00264 long int nv = h2.nVib_hi[iElecLo]; 00265 if( iElecLo==iElecHi ) 00266 nv = iVibHi; 00267 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00268 { 00269 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00270 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00271 nr = iRotHi-1; 00272 00273 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00274 { 00275 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00276 { 00277 /* population of lower level with correction for stim emission */ 00278 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->PopOpc = 00279 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop - 00280 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->Pop* 00281 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g / 00282 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g; 00283 00284 /* following two heat exchange excitation, deexcitation */ 00285 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cool = 0.; 00286 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.heat = 0.; 00287 00288 /* intensity of line */ 00289 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->xIntensity = 0.; 00290 00291 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->phots = 0.; 00292 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots = 0.; 00293 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ColOvTot = 0.; 00294 } 00295 } 00296 } 00297 } 00298 } 00299 } 00300 } 00301 hmi.H2_photodissoc_BigH2_H2s = 0.; 00302 hmi.H2_photodissoc_BigH2_H2g = 0.; 00303 hmi.HeatH2Dish_BigH2 = 0.; 00304 hmi.HeatH2Dexc_BigH2 = 0.; 00305 hmi.deriv_HeatH2Dexc_BigH2 = 0.; 00306 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.; 00307 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.; 00308 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.; 00309 return; 00310 } 00311 00312 /*mole_H2_LTE sets Boltzmann factors and LTE unit population of large H2 molecular */ 00313 void mole_H2_LTE( void ) 00314 { 00315 /* used to recall the temperature used for last set of Boltzmann factors */ 00316 static double TeUsedBoltz = -1.; 00317 double part_fun; 00318 long int iElec , iVib , iRot; 00319 00320 DEBUG_ENTRY( "mole_H2_LTE()" ); 00321 00322 /* do we need to update the Boltzmann factors and unit LTE populations? */ 00323 if( ! fp_equal( phycon.te, TeUsedBoltz ) ) 00324 { 00325 part_fun = 0.; 00326 TeUsedBoltz = phycon.te; 00327 /* loop over all levels setting H2_Boltzmann and deriving partition function */ 00328 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 00329 { 00330 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 00331 { 00332 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 00333 { 00334 H2_Boltzmann[iElec][iVib][iRot] = 00335 /* energy is relative to lowest level in the molecule, v=0, J=0, 00336 * so Boltzmann factor is relative to this level */ 00337 sexp( energy_wn[iElec][iVib][iRot] / phycon.te_wn ); 00338 /* sum the partition function - Boltzmann factor times statistical weight */ 00339 part_fun += H2_Boltzmann[iElec][iVib][iRot] * H2_stat[iElec][iVib][iRot]; 00340 ASSERT( part_fun > 0 ); 00341 } 00342 } 00343 } 00344 /* have partition function, set H2_populations_LTE (populations for unit H2 density) */ 00345 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 00346 { 00347 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 00348 { 00349 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 00350 { 00351 /* these are the H2 LTE populations for a unit H2 density - 00352 * these populations will sum up to unity */ 00353 H2_populations_LTE[iElec][iVib][iRot] = 00354 H2_Boltzmann[iElec][iVib][iRot] * 00355 H2_stat[iElec][iVib][iRot] / part_fun; 00356 /*if( iElec==0 && iVib < 2) 00357 fprintf(ioQQQ,"DEBUG LTE pop\t%i\t%i\t%e\n", 00358 iVib,iRot,H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total);*/ 00359 } 00360 } 00361 } 00362 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 00363 fprintf(ioQQQ, 00364 "mole_H2_LTE set H2_Boltzmann factors, T=%.2f, partition function is %.2f\n", 00365 phycon.te, 00366 part_fun); 00367 } 00368 00369 return; 00370 } 00371 00372 /*H2_init_coreload one time initialization */ 00373 void H2_init_coreload( void ) 00374 { 00375 /* the order of the electronic states is 00376 * X, B, C+, C-, B', D+, and D- */ 00377 /* this will be the number of vibration levels within each electronic */ 00378 /* number of vib states within electronic states from 00379 * >>refer H2 energies Abgrall, */ 00380 long int nVib_hi_init[N_H2_ELEC] = {14 , 37 , 13 , 13, 9, 2 , 2}; 00381 00382 /* this gives the first rotational state for each electronic state - J=0 does 00383 * not exist when Lambda = 1 */ 00384 long int Jlowest_init[N_H2_ELEC] = {0 , 0 , 1 , 1 , 0 , 1 , 1 }; 00385 00386 /* number of rotation levels within each electronic - vib */ 00387 /*lint -e785 too few init for aggregate */ 00388 long int nRot_hi_init[N_H2_ELEC][50]= 00389 /* ground, X */ 00390 { {31, 30, 28, 27, 25, 00391 23, 22, 20, 18, 16, 00392 14, 12, 10, 7, 3 } , 00393 /* B */ 00394 {25,25,25,25,25,25,25,25, 25,25, 00395 25,25,25,25,25,25,25,25, 25,25, 00396 25,25,25,25,25,25,25,25, 23,21, 00397 19,17,15,15,11,9,7, 7}, 00398 /* C plus */ 00399 { 25, 25, 25, 25, 24, 23, 21, 19, 17, 14, 12, 10, 6, 2 }, 00400 /* C minus (the same) */ 00401 { 25, 25, 25, 25, 24, 23, 21, 19, 17, 15, 13, 10, 7, 2 }, 00402 /* B primed */ 00403 {19,17, 14, 12, 9, 8, 7, 7, 4, 1 }, 00404 /* D plus */ 00405 {13, 10, 5}, 00406 /* D minus */ 00407 {25 , 25 ,25 } } 00408 ; 00409 /*lint +e785 too few init for aggregate */ 00410 /* one time init */ 00411 00412 long int iElec; 00413 00414 DEBUG_ENTRY( "H2_init_coreload()" ); 00415 00416 /* the order of the electronic states is 00417 * X, B, C+, C-, B', D+, and D- */ 00418 /* this will be the number of vibration levels within each electronic */ 00419 /* number of vib states within electronic states from 00420 * >>refer H2 energies Abgrall, */ 00421 for( iElec=0; iElec<N_H2_ELEC; ++iElec ) 00422 { 00423 int iVib; 00424 h2.nVib_hi[iElec] = nVib_hi_init[iElec]; 00425 h2.Jlowest[iElec] = Jlowest_init[iElec]; 00426 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 00427 { 00428 h2.nRot_hi[iElec][iVib] = nRot_hi_init[iElec][iVib]; 00429 /*fprintf(ioQQQ,"DEBUG h2 set %li %li %li \n", 00430 iElec , 00431 iVib , 00432 h2.nRot_hi[iElec][iVib] );*/ 00433 } 00434 } 00435 strcpy( chH2ColliderLabels[0] , "H0" ); 00436 strcpy( chH2ColliderLabels[1] , "He" ); 00437 strcpy( chH2ColliderLabels[2] , "H2 o" ); 00438 strcpy( chH2ColliderLabels[3] , "H2 p" ); 00439 strcpy( chH2ColliderLabels[4] , "H+" ); 00440 00441 return; 00442 } 00443 00444 /*H2_init - called to initialize things from cdInit */ 00445 void H2_Init(void) 00446 { 00447 00448 DEBUG_ENTRY( "H2_Init()" ); 00449 00450 /* the number of electronic quantum states to include. 00451 * To do both Lyman and Werner bands want nelec = 3, 00452 * default is to do all bands included */ 00453 mole.n_h2_elec_states = N_H2_ELEC; 00454 h2.nCallH2_this_zone = 0; 00455 return; 00456 } 00457 00458 00459 /*H2_Reset called by IterRestart to reset variables that are needed after an iteration */ 00460 void H2_Reset( void ) 00461 { 00462 00463 DEBUG_ENTRY( "H2_Reset()" ); 00464 00465 if(mole.nH2_TRACE) 00466 fprintf(ioQQQ, 00467 "\n***************H2_Reset called, resetting nCallH2_this_iteration, zone %.2f iteration %li\n", 00468 fnzone, 00469 iteration ); 00470 00471 /* number of times large molecules evaluated in this iteration, 00472 * is false if never evaluated, on next evaluation will start with LTE populations */ 00473 nCallH2_this_iteration = 0; 00474 00475 /* these remember the largest and smallest factors needed to 00476 * renormalize the H2 chemistry */ 00477 h2.renorm_max = 1.; 00478 h2.renorm_min = 1.; 00479 00480 /* counters used by H2_itrzn to find number of calls of h2 per zone */ 00481 nH2_pops = 0; 00482 nH2_zone = 0; 00483 /* this is used to establish zone number for evaluation of number of levels in matrix */ 00484 nzone_nlevel_set = 0; 00485 00486 nzoneAsEval = -1; 00487 iterationAsEval = -1; 00488 00489 /* zero out array used to save emission line intensities */ 00490 H2_SaveLine.zero(); 00491 00492 return; 00493 } 00494 00495 /*H2_Zero zero out vars in the large H2 molecule, called from zero 00496 * before any commands are parsed 00497 * NB - this routine is called before space allocated - must not zero out 00498 * any allocated arrays */ 00499 void H2_Zero( void ) 00500 { 00501 00502 DEBUG_ENTRY( "H2_Zero()" ); 00503 00504 /* this is the smallest ratio of H2/H where we will bother with the large H2 molecule 00505 * this value was chosen so that large mole used at very start of TH85 standard PDR, 00506 * NB - this appears in headinfo and must be updated there if changed here */ 00507 /* >>chng 03 jun 02, from 1e-6 to 1e-8 - in orion veil large H2 turned on half way 00508 * across, and Solomon process was very fast since all lines optically thin. correct 00509 * result has some shielding, so reset to lower value so that H2 comes on sooner. */ 00510 mole.H2_to_H_limit = 1e-8; 00511 00512 h2.lgH2ON = false; 00513 /* flag to force using LTE level populations */ 00514 mole.lgH2_LTE = false; 00515 00516 /* counters used by H2_itrzn to find number of calls of h2 per zone */ 00517 nH2_pops = 0; 00518 nH2_zone = 0; 00519 /* this is used to establish zone number for evaluation of number of levels in matrix */ 00520 nzone_nlevel_set = 0; 00521 00522 /* these remember the largest and smallest factors needed to 00523 * renormalize the H2 chemistry */ 00524 h2.renorm_max = 1.; 00525 h2.renorm_min = 1.; 00526 00527 nCallH2_this_iteration = 0; 00528 h2.ortho_density = 0.; 00529 h2.para_density = 0.; 00530 00531 hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.; 00532 hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.; 00533 00534 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.; 00535 hmi.H2_photodissoc_BigH2_H2s = 0.; 00536 hmi.H2_photodissoc_BigH2_H2g = 0.; 00537 00538 hmi.HeatH2Dexc_BigH2 = 0.; 00539 00540 /* say that H2 has never been computed */ 00541 hmi.lgBigH2_evaluated = false; 00542 00543 hmi.lgH2_Thermal_BigH2 = true; 00544 hmi.lgH2_Chemistry_BigH2 = true; 00545 00546 if( !lgH2_READ_DATA ) 00547 { 00548 /* the number of electronic levels in the H2 molecule, 00549 * to just do the Lyman and Werner bands set to 3 - 00550 * reset with atom h2 levels command, 00551 * default is all levels with data */ 00552 mole.n_h2_elec_states = N_H2_ELEC; 00553 } 00554 00555 /* the number of levels used in the matrix solution 00556 * of the level H2_populations - set with atom h2 matrix nlevel, 00557 * >>chng 04 oct 05, make default 30 levels 00558 * >>chng 04 dec 23, make default 70 levels */ 00559 nXLevelsMatrix = 70; 00560 00561 /* this is used to establish zone number for evaluation of number of levels in matrix */ 00562 nzone_nlevel_set = -1; 00563 00564 nzoneAsEval = -1; 00565 iterationAsEval = -1; 00566 return; 00567 }