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 /*H2_ContPoint set the ipCont struc element for the H2 molecule, called by ContCreatePointers */ 00004 /*H2_Accel radiative acceleration due to H2 */ 00005 /*H2_RadPress rad pressure due to h2 lines called in PresTotCurrent */ 00006 /*H2_InterEnergy internal energy of H2 called in PresTotCurrent */ 00007 /*H2_RT_diffuse do emission from H2 - called from RT_diffuse */ 00008 /*H2_itrzn - average number of H2 pop evaluations per zone */ 00009 /*H2_RTMake do RT for H2 - called from RT_line_all */ 00010 /*H2_RT_tau_inc increment optical depth for the H2 molecule, called from RT_tau_inc */ 00011 /*H2_LineZero initialize optical depths in H2, called from RT_tau_init */ 00012 /*H2_RT_tau_reset the large H2 molecule, called from RT_tau_reset */ 00013 /*H2_Colden maintain H2 column densities within X */ 00014 /*H2_LevelPops do level H2_populations for H2, called by Hydrogenic */ 00015 /*H2_Level_low_matrix evaluate CO rotation cooling */ 00016 /*H2_cooling evaluate cooling and heating due to H2 molecule */ 00017 /*H2_X_coll_rate_evaluate find collisional rates within X */ 00018 /*cdH2_colden return column density in H2, negative -1 if cannot find state, 00019 * header is cddrive */ 00020 /*H2_DR choose next zone thickness based on H2 big molecule */ 00021 /* turn this flag on to do minimal debug print of pops */ 00022 #define PRT_POPS false 00023 /* this is limit to number of loops over H2 pops before giving up */ 00024 #define LIM_H2_POP_LOOP 100 00025 /* this is a typical dissociation cross section (cm2) for H2 + Hnu -> 2H + ke */ 00026 /* >>chng 05 may 11, had been 2.5e-19 */ 00027 #define H2_DISS_ALLISON_DALGARNO 6e-19f 00028 #include "cddefines.h" 00029 #include "cddrive.h" 00030 #include "physconst.h" 00031 #include "taulines.h" 00032 #include "atoms.h" 00033 #include "conv.h" 00034 #include "secondaries.h" 00035 #include "pressure.h" 00036 #include "trace.h" 00037 #include "hmi.h" 00038 #include "hextra.h" 00039 #include "rt.h" 00040 #include "radius.h" 00041 #include "ipoint.h" 00042 #include "phycon.h" 00043 #include "thermal.h" 00044 #include "dense.h" 00045 #include "rfield.h" 00046 #include "lines_service.h" 00047 #include "mole.h" 00048 #include "h2.h" 00049 #include "h2_priv.h" 00050 00051 /* this counts how many times we go through the H2 level H2_populations loop */ 00052 static long int loop_h2_pops; 00053 00054 realnum H2_te_hminus[nTE_HMINUS] = {10.,30.,100.,300.,1000.,3000.,10000.}; 00055 00056 /* this will contain a vector for collisions within the X ground electronic state, 00057 * CollRateFit[vib_up][rot_up][vib_lo][rot_lo][coll_type][3] */ 00058 static realnum collider_density[N_X_COLLIDER]; 00059 static realnum collider_density_total_not_H2; 00060 00061 /* this integer is added to rotation quantum number J for the test of whether 00062 * a particular J state is ortho or para - the state is ortho if J+below is odd, 00063 * and para if J+below is even */ 00064 int H2_nRot_add_ortho_para[N_H2_ELEC] = {0 , 1 , 1 , 0, 1, 1 , 0}; 00065 00066 /* dissociation energies (cm-1) for each electronic state, from 00067 * >>refer H2 energies Sharp, T. E., 1971, Atomic Data, 2, 119 00068 * energy for H_2 + hnu -> 2H, 00069 * note that energies for excited states are in the Lyman continuum 00070 * (1 Ryd = 109670 cm-1) */ 00071 /* >>chng 02 oct 08, improved energies */ 00072 double H2_DissocEnergies[N_H2_ELEC] = 00073 { 36118.11, 118375.6, 118375.6, 118375.6, 118375.6,133608.6,133608.6 }; 00074 /* original values 00075 { 36113., 118372., 118372., 118372., 118372.,0.,0. };*/ 00076 00077 double exp_disoc; 00078 00079 /*H2_X_coll_rate_evaluate find collisional rates within X - 00080 * this is one time upon entry into H2_LevelPops */ 00081 STATIC void H2_X_coll_rate_evaluate( void ) 00082 { 00083 long int nColl, 00084 ipHi , 00085 ipLo, 00086 ip, 00087 iVibHi, 00088 iRotHi, 00089 iVibLo, 00090 iRotLo; 00091 realnum colldown; 00092 double exph2_big, 00093 exph2_big_0, 00094 rel_pop_LTE_H2_big; 00095 00096 DEBUG_ENTRY( "H2_X_coll_rate_evaluate()" ); 00097 00098 /* set collider density 00099 * the colliders are: 00100 * [0] = H 00101 * [1], [5] = He (old and new cs data) 00102 * [2] = H2 ortho 00103 * [3] = H2 para 00104 * [4] = H+ + H3+ */ 00105 /* atomic hydrogen */ 00106 collider_density[0] = dense.xIonDense[ipHYDROGEN][0]; 00107 /* all ortho h2 */ 00108 /* He - H2 */ 00109 collider_density[1] = dense.xIonDense[ipHELIUM][0]; 00110 /* H2 - H2(ortho) */ 00111 collider_density[2] = (realnum)h2.ortho_density; 00112 /* all para H2 */ 00113 collider_density[3] = (realnum)h2.para_density; 00114 /* protons - ionized hydrogen */ 00115 collider_density[4] = dense.xIonDense[ipHYDROGEN][1]; 00116 /* H3+ - assume that H3+ has same rates as proton */ 00117 collider_density[4] += hmi.Hmolec[ipMH3p]; 00118 00119 ASSERT( fp_equal(hmi.H2_total ,collider_density[2]+collider_density[3]) ); 00120 00121 /* this is total density of all colliders, is only used for collisional dissociation 00122 * rates for H2 are not included here, will be added separately*/ 00123 collider_density_total_not_H2 = collider_density[0] + 00124 collider_density[1] + collider_density[4] + 00125 (realnum)dense.eden; 00126 00127 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 00128 { 00129 fprintf(ioQQQ," Collider densities are:"); 00130 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 00131 { 00132 fprintf(ioQQQ,"\t%.3e", collider_density[nColl]); 00133 } 00134 fprintf(ioQQQ,"\n"); 00135 } 00136 00137 for( ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi ) 00138 { 00139 H2_X_source[ipHi] = 0.; 00140 H2_X_sink[ipHi] = 0.; 00141 } 00142 H2_X_coll_rate.zero(); 00143 00144 /*>>chng 05 sep 18, GS, determine LTE populations for H2+ e = H- + H*/ 00145 exp_disoc = sexp(H2_DissocEnergies[0]/phycon.te_wn); 00146 exph2_big_0 = exp_disoc/SDIV(H2_Boltzmann[0][0][0]); 00147 /*>>chng 05 oct 17, GS, (2*m_e/m_p)^3/2 = 3.634e-5 */ 00148 rel_pop_LTE_H2_big =SAHA/SDIV((phycon.te32*exph2_big_0))*(H2_stat[0][0][0]/(2.*2.))*3.634e-5; 00149 /* now find rates for all collisions within X */ 00150 /* this is special J=1 to J=0 collision, which is only fast at 00151 * very low grain temperatures 00152 H2_X_coll_rate[1][0] = 00153 (realnum)(hmi.rate_grain_h2_J1_to_J0);*/ 00154 00155 /* count formation from grains and H- as a collisional formation process */ 00156 /* cm-3 s-1, evaluated in mole_H2_form */ 00157 H2_X_source[0] += H2_X_formation[0][0]; 00158 /*>>chng 05 sep 18, GS, H2+ e = H- + H, unit s-1 00159 * H2_X_Hmin_back[iVib][iRot] is resolved formation rate, units cm3 s-1 */ 00160 H2_X_sink[0] += (realnum)(H2_X_Hmin_back[0][0]*hmi.rel_pop_LTE_Hmin/ 00161 SDIV(rel_pop_LTE_H2_big)*dense.eden); 00162 /* this represents collisional dissociation into continuum of X, 00163 * rates are just guesses */ 00164 H2_X_sink[0] += collider_density_total_not_H2 * 00165 H2_coll_dissoc_rate_coef[0][0] * mole.lgColl_deexec_Calc; 00166 /*>>chng 05 jul 20, GS, collisional dissociation with H2g and H2s are added here*/ 00167 H2_X_sink[0] += hmi.H2_total* 00168 H2_coll_dissoc_rate_coef_H2[0][0] * mole.lgColl_deexec_Calc; 00169 00170 /* this is cosmic ray or secondary photodissociation into mainly H2+ */ 00171 H2_X_sink[0] += secondaries.csupra[ipHYDROGEN][0]*2.02f * mole.lgColl_deexec_Calc; 00172 00173 /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H- */ 00174 H2_X_sink[0] += (realnum)(3.9e-21 * hextra.cryden_ov_background * co.lgUMISTrates); 00175 00176 /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H + e */ 00177 H2_X_sink[0] += (realnum)(2.2e-19 * hextra.cryden_ov_background * co.lgUMISTrates); 00178 00179 /*>>chng 05 sep 26, GS, H2 + CR -> H+ + H + e */ 00180 H2_X_sink[0] += (realnum)(secondaries.csupra[ipHYDROGEN][0]*0.0478); 00181 00182 /* collisional dissociation by non-thermal electrons 00183 * and cosmic rays to the triplet state of H2 (a,b,c ) from 00184 *>>>refer Dalgarno, Yan, & Liu 1999 ApJs 00185 * rates depend on energy of secondary electrons, this is an estimate 00186 * from fig 4b around incident electron energy 15 to 20 eV. 00187 * The cross section of state b is divided by the cross-section of Lya 00188 * (the ratio is ~ 10)and then multiplied by the cosmic ray excitation 00189 * rate of Lya (secondaries.x12tot) the triplet state of H2 (b ) 00190 *>>chng 07 apr 08, GS from 3 to 10 after study of Dalgarno et al */ 00191 H2_X_sink[0] += 10.f*secondaries.x12tot; 00192 00193 /*>>chng 05 jul 07, GS, ionization by hard photons are added to be consistent with chemistry-network */ 00194 H2_X_sink[0] += (realnum)hmi.H2_photoionize_rate; 00195 00196 /* rate (s-1) out of this level, this is dissociation into the continuum of singlet B and C states*/ 00197 H2_X_sink[0] += rfield.flux_accum[H2_ipPhoto[0][0]-1]*H2_DISS_ALLISON_DALGARNO; 00198 00199 /*>>chng 05 sept 28, GS, H2 + H+ = H2+ + H; note that H2 is destroyed from v=0, J=0, 00200 * forward and backward reactions are not in detailed balance, astro-ph/0404288, final unit s-1 00201 * ihmi.rh2h2p is a rate coefficient, cm3 s-1, and needs a density to become 00202 * the sink inverse lifetime. The process is H2(v=0, J=0) + H+ -> H0 + H2+ */ 00203 H2_X_sink[0] += (realnum)(hmi.rh2h2p*dense.xIonDense[ipHYDROGEN][1]); 00204 00205 /* now find total coll rates - this loop is over the energy-sorted levels */ 00206 ipHi = nLevels_per_elec[0]; 00207 while( (--ipHi) > 0 ) 00208 { 00209 /* array of energy sorted indices within X */ 00210 ip = H2_ipX_ener_sort[ipHi]; 00211 iVibHi = ipVib_H2_energy_sort[ip]; 00212 iRotHi = ipRot_H2_energy_sort[ip]; 00213 00214 /*>>chng 05 sep 18, GS, determine LTE populations for H2+ e = H- + H, unit s-1*/ 00215 exph2_big = exp_disoc/SDIV(H2_Boltzmann[0][iVibHi][iRotHi]); 00216 /* >>chng 05 oct 13, replace 1 with stat wght of level */ 00217 /*>>chng 05 oct 17, GS, (2*m_e/m_p)^3/2 = 3.634e-5 */ 00218 rel_pop_LTE_H2_big = SAHA/SDIV((phycon.te32*exph2_big))*(H2_stat[0][iVibHi][iRotHi]/(2.*2.))*3.634e-5; 00219 /* formation */ 00220 /* count formation from grains and H- as a collisional formation process */ 00221 /* cm-3 s-1, evaluated in mole_H2_form */ 00222 H2_X_source[ipHi] += H2_X_formation[iVibHi][iRotHi]; 00223 /*>>chng 05 sep 18, GS, H2+ e = H- + H*, H2_X_Hmin_back has units cm3 s-1 so final units s-1 */ 00224 H2_X_sink[ipHi] += (realnum)(H2_X_Hmin_back[iVibHi][iRotHi]*hmi.rel_pop_LTE_Hmin/ 00225 SDIV(rel_pop_LTE_H2_big)*dense.eden); 00226 /* this represents collisional dissociation into continuum of X, 00227 * rates are just guesses */ 00228 H2_X_sink[ipHi] += collider_density_total_not_H2 * 00229 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] * mole.lgColl_deexec_Calc; 00230 00231 /*>>chng 05 jul 20, GS, collisional dissociation with H2g and H2s are added here*/ 00232 H2_X_sink[ipHi] += hmi.H2_total* 00233 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] * mole.lgColl_deexec_Calc; 00234 00235 /* this is cosmic ray or secondary photodissociation into mainly H2+ */ 00236 H2_X_sink[ipHi] += secondaries.csupra[ipHYDROGEN][0]*2.02f * mole.lgColl_deexec_Calc; 00237 00238 /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H- */ 00239 H2_X_sink[ipHi] += (realnum)(3.9e-21 * hextra.cryden_ov_background * co.lgUMISTrates); 00240 00241 /*>>chng 05 sep 26, GS, H2 + CRP -> H+ + H + e */ 00242 H2_X_sink[ipHi] += (realnum)(2.2e-19 * hextra.cryden_ov_background * co.lgUMISTrates); 00243 00244 /*>>chng 05 sep 26, GS, H2 + CR -> H+ + H + e */ 00245 H2_X_sink[ipHi] += (realnum)(secondaries.csupra[ipHYDROGEN][0]*0.0478); 00246 00247 /* collisional dissociation by non-thermal electrons 00248 * and cosmic rays to the triplet state of H2 (a,b,c ) from 00249 *>>>refer Dalgarno, Yan, & Liu 1999 ApJs 00250 * rates depend on energy of secondary electrons, this is an estimate 00251 *from fig 4b around incident electron energy 15 to 20 eV. 00252 * The cross section of state b is divided by the cross-section of Lya 00253 * (the ratio is ~ 10)and then multiplied by the cosmic ray excitation 00254 * rate of Lya (secondaries.x12tot) the triplet state of H2 (b ) 00255 *>>chng 07 apr 08, GS from 3 to 10 after study of Dalgarno et al */ 00256 H2_X_sink[ipHi] += 10.f*secondaries.x12tot; 00257 00258 /*>>chng 05 jul 07, GS, ionization by hard photons are added to be consistent with chemistry-network */ 00259 H2_X_sink[ipHi] += (realnum)hmi.H2_photoionize_rate; 00260 00261 /* rate (s-1) out of this level */ 00262 H2_X_sink[ipHi] += rfield.flux_accum[H2_ipPhoto[iVibHi][iRotHi]-1]*H2_DISS_ALLISON_DALGARNO; 00263 00264 if( mole.lgColl_deexec_Calc ) 00265 { 00266 /* excitation within X due to thermal particles */ 00267 for( ipLo=0; ipLo<ipHi; ++ipLo ) 00268 { 00269 /* array of energy sorted indices within X */ 00270 ip = H2_ipX_ener_sort[ipLo]; 00271 iVibLo = ipVib_H2_energy_sort[ip]; 00272 iRotLo = ipRot_H2_energy_sort[ip]; 00273 00274 /* collisional interactions with upper levels within X */ 00275 colldown = 0.; 00276 mr5ci H2CollRate = H2_CollRate.begin(iVibHi,iRotHi,iVibLo,iRotLo); 00277 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 00278 { 00279 /* downward collision rate, units s-1 */ 00280 colldown += H2CollRate[nColl]*collider_density[nColl]; 00281 ASSERT( H2CollRate[nColl]*collider_density[nColl] >= 0. ); 00282 } 00283 /* rate in from upper level, units cm-3 s-1 */ 00284 H2_X_coll_rate[ipHi][ipLo] += colldown; 00285 }/* end loop over ipLo */ 00286 } 00287 }/* end loop over ipHi */ 00288 return; 00289 } 00290 00291 /*H2_itrzn - average number of H2 pop evaluations per zone */ 00292 double H2_itrzn( void ) 00293 { 00294 if( h2.lgH2ON && nH2_zone>0 ) 00295 { 00296 return( (double)nH2_pops / (double)nH2_zone ); 00297 } 00298 else 00299 { 00300 return 0.; 00301 } 00302 } 00303 00304 /* set the ipCont struc element for the H2 molecule, called by ContCreatePointers */ 00305 void H2_ContPoint( void ) 00306 { 00307 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00308 00309 if( !h2.lgH2ON ) 00310 return; 00311 00312 DEBUG_ENTRY( "H2_ContPoint()" ); 00313 00314 /* set array index for line energy within continuum array */ 00315 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00316 { 00317 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00318 { 00319 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00320 { 00321 /* now the lower levels */ 00322 /* NB - X is the only lower level considered here, since we are only 00323 * concerned with excited electronic levels as a photodissociation process 00324 * code exists to relax this assumption - simply change following to iElecHi */ 00325 long int lim_elec_lo = 0; 00326 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00327 { 00328 /* want to include all vibration states in lower level if different electronic level, 00329 * but only lower vibration levels if same electronic level */ 00330 long int nv = h2.nVib_hi[iElecLo]; 00331 if( iElecLo==iElecHi ) 00332 nv = iVibHi; 00333 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00334 { 00335 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00336 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00337 nr = iRotHi-1; 00338 00339 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00340 { 00341 /* >>chng 05 feb 07, use lgH2_line_exists */ 00342 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00343 { 00344 ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul > 0. ); 00345 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont = 00346 ipLineEnergy( 00347 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD , 00348 "H2 " , 0 ); 00349 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ipFine = 00350 ipFineCont( 00351 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD ); 00352 } 00353 } 00354 } 00355 } 00356 } 00357 } 00358 } 00359 return; 00360 } 00361 00362 /* ===================================================================== */ 00363 /* radiative acceleration due to H2 called in rt_line_driving */ 00364 double H2_Accel(void) 00365 { 00366 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00367 double h2_drive; 00368 00369 /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */ 00370 if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ ) 00371 return 0.; 00372 00373 DEBUG_ENTRY( "H2_Accel()" ); 00374 00375 /* this routine computes the line driven radiative acceleration 00376 * due to H2 molecule*/ 00377 00378 h2_drive = 0.; 00379 /* loop over all possible lines */ 00380 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00381 { 00382 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00383 { 00384 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00385 { 00386 /* now the lower levels */ 00387 /* NB - X is the only lower level considered here, since we are only 00388 * concerned with excited electronic levels as a photodissociation process 00389 * code exists to relax this assumption - simply change following to iElecHi */ 00390 long int lim_elec_lo = 0; 00391 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00392 { 00393 /* want to include all vibration states in lower level if different electronic level, 00394 * but only lower vibration levels if same electronic level */ 00395 long int nv = h2.nVib_hi[iElecLo]; 00396 if( iElecLo==iElecHi ) 00397 nv = iVibHi; 00398 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00399 { 00400 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00401 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00402 nr = iRotHi-1; 00403 00404 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00405 mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00406 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00407 { 00408 /* >>chng 03 feb 14, from !=0 to >0 */ 00409 /* >>chng 05 feb 07, use lgH2_line_exists */ 00410 if( *lgH2le++ ) 00411 { 00412 ASSERT( H2L->ipCont > 0 ); 00413 h2_drive += H2L->Emis->pump*H2L->EnergyErg*H2L->Emis->PopOpc; 00414 } 00415 ++H2L; 00416 } 00417 } 00418 } 00419 } 00420 } 00421 } 00422 return h2_drive; 00423 } 00424 00425 /* ===================================================================== */ 00426 /* rad pressure due to H2 lines called in PresTotCurrent */ 00427 double H2_RadPress(void) 00428 { 00429 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00430 double press; 00431 00432 /* will be used to check on size of opacity, was capped at this value */ 00433 realnum smallfloat=SMALLFLOAT*10.f; 00434 00435 /* radiation pressure sum is expensive - do not evaluate if we did not 00436 * bother evaluating large molecule */ 00437 if( !h2.lgH2ON || !h2.nCallH2_this_zone ) 00438 return 0.; 00439 00440 DEBUG_ENTRY( "H2_RadPress()" ); 00441 00442 press = 0.; 00443 /* loop over all possible lines */ 00444 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00445 { 00446 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00447 { 00448 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00449 { 00450 /* now the lower levels - X is the only lower level considered 00451 * here, we only do excited electronic levels as a 00452 * photodissociation process - code exists to relax this 00453 * assumption - simply change following to iElecHi */ 00454 long int lim_elec_lo = 0; 00455 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00456 { 00457 /* want to include all vibration states in lower level if different electronic level, 00458 * but only lower vibration levels if same electronic level */ 00459 long int nv = h2.nVib_hi[iElecLo]; 00460 if( iElecLo==iElecHi ) 00461 nv = iVibHi; 00462 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00463 { 00464 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00465 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00466 nr = iRotHi-1; 00467 00468 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00469 mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00470 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00471 { 00472 if( *lgH2le++ ) 00473 { 00474 ASSERT( H2L->ipCont > 0 ); 00475 00476 if( H2L->Hi->Pop > smallfloat && H2L->Emis->PopOpc > smallfloat ) 00477 { 00478 double RadPres1 = PressureRadiationLine( &(*H2L), 1.0 ); 00479 press += RadPres1; 00480 } 00481 } 00482 ++H2L; 00483 } 00484 } 00485 } 00486 } 00487 } 00488 } 00489 00490 if(mole.nH2_TRACE >= mole.nH2_trace_full) 00491 fprintf(ioQQQ, 00492 " H2_RadPress returns, radiation pressure is %.2e\n", 00493 press ); 00494 return press; 00495 } 00496 00497 #if 0 00498 /* ===================================================================== */ 00499 /* internal energy of H2 called in PresTotCurrent */ 00500 double H2_InterEnergy(void) 00501 { 00502 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00503 double energy; 00504 00505 /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */ 00506 if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ ) 00507 return 0.; 00508 00509 DEBUG_ENTRY( "H2_InterEnergy()" ); 00510 00511 broken(); /* the code below contains serious bugs! It is supposed to implement a loop 00512 * over all states in order to evaluate the potential energy stored in 00513 * excited states of H2. The code below however implements loops over all 00514 * combinations of upper and lower levels! */ 00515 00516 energy = 0.; 00517 /* loop over all possible lines */ 00518 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00519 { 00520 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00521 { 00522 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00523 { 00524 double H2popHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop; 00525 00526 /* now the lower levels 00527 * NB - X is the only lower level considered here, since we are only 00528 * concerned with excited electronic levels as a photodissociation process 00529 * code exists to relax this assumption - simply change following to iElecHi */ 00530 long int lim_elec_lo = 0; 00531 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00532 { 00533 /* want to include all vibration states in lower level if different electronic level, 00534 * but only lower vibration levels if same electronic level */ 00535 long int nv = h2.nVib_hi[iElecLo]; 00536 if( iElecLo==iElecHi ) 00537 nv = iVibHi; 00538 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00539 { 00540 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00541 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00542 nr = iRotHi-1; 00543 00544 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00545 mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00546 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00547 { 00548 /* >>chng 03 feb 14, from !=0 to >0 */ 00549 /* >>chng 05 feb 07, use lgH2_line_exists */ 00550 if( *lgH2le++ ) 00551 { 00552 energy += H2popHi*H2L->EnergyErg; 00553 } 00554 ++H2L; 00555 } 00556 } 00557 } 00558 } 00559 } 00560 } 00561 return energy; 00562 } 00563 #endif 00564 00565 /*H2_RT_diffuse do emission from H2 - called from RT_diffuse */ 00566 void H2_RT_diffuse(void) 00567 { 00568 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00569 00570 if( !h2.lgH2ON || !h2.nCallH2_this_zone ) 00571 return; 00572 00573 DEBUG_ENTRY( "H2_RT_diffuse()" ); 00574 00575 /* loop over all possible lines */ 00576 /* NB - this loop does not include the electronic lines */ 00577 for( iElecHi=0; iElecHi<1; ++iElecHi ) 00578 { 00579 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00580 { 00581 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00582 { 00583 /* now the lower levels */ 00584 /* NB - X is the only lower level considered here, since we are only 00585 * concerned with excited electronic levels as a photodissociation process 00586 * code exists to relax this assumption - simply change following to iElecHi */ 00587 long int lim_elec_lo = 0; 00588 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00589 { 00590 /* want to include all vibration states in lower level if different electronic level, 00591 * but only lower vibration levels if same electronic level */ 00592 long int nv = h2.nVib_hi[iElecLo]; 00593 if( iElecLo==iElecHi ) 00594 nv = iVibHi; 00595 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00596 { 00597 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00598 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00599 nr = iRotHi-1; 00600 00601 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00602 { 00603 /* >>chng 03 feb 14, from !=0 to > 0 */ 00604 /* >>chng 05 feb 07, use lgH2_line_exists */ 00605 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00606 { 00607 outline( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ); 00608 } 00609 } 00610 } 00611 } 00612 } 00613 } 00614 } 00615 return; 00616 } 00617 00618 /* do RT for H2 - called from RT_line_all */ 00619 void H2_RTMake( bool lgDoEsc , bool lgUpdateFineOpac ) 00620 { 00621 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00622 00623 if( !h2.lgH2ON ) 00624 return; 00625 00626 DEBUG_ENTRY( "H2_RTMake()" ); 00627 00628 /* this routine drives calls to make RT relations for H2 molecule */ 00629 /* loop over all possible lines */ 00630 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00631 { 00632 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00633 { 00634 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00635 { 00636 /* now the lower levels */ 00637 /* NB - X is the only lower level considered here, since we are only 00638 * concerned with excited electronic levels as a photodissociation process 00639 * code exists to relax this assumption - simply change following to iElecHi */ 00640 long int lim_elec_lo = 0; 00641 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00642 { 00643 /* want to include all vibration states in lower level if different electronic level, 00644 * but only lower vibration levels if same electronic level */ 00645 long int nv = h2.nVib_hi[iElecLo]; 00646 if( iElecLo==iElecHi ) 00647 nv = iVibHi; 00648 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00649 { 00650 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00651 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00652 nr = iRotHi-1; 00653 00654 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00655 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00656 { 00657 /* >>chng 03 feb 14, change test from !=0 to >0 */ 00658 /* >>chng 05 feb 07, use lgH2_line_exists */ 00659 if( *lgH2le++ ) 00660 { 00661 /* >>chng 03 jun 18, added 4th parameter in call to this routine - says to not 00662 * include self-shielding of line across this zone. This introduces a dr dependent 00663 * variation in the line pumping rate, which made H2 abundance fluctuate due to 00664 * Solomon process having slight dr-caused mole. */ 00665 RT_line_one( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo], lgDoEsc, lgUpdateFineOpac, false ); 00666 } 00667 } 00668 } 00669 } 00670 } 00671 } 00672 } 00673 00674 /* debug print population weighted transition probability */ 00675 { 00676 enum {DEBUG_LOC=false}; 00677 if( DEBUG_LOC ) 00678 { 00679 double sumpop = 0.; 00680 double sumpopA = 0.; 00681 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00682 { 00683 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00684 { 00685 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00686 { 00687 /* now the lower levels */ 00688 /* NB - X is the only lower level considered here, since we are only 00689 * concerned with excited electronic levels as a photodissociation process 00690 * code exists to relax this assumption - simply change following to iElecHi */ 00691 iElecLo = 0; 00692 /* want to include all vibration states in lower level if different electronic level, 00693 * but only lower vibration levels if same electronic level */ 00694 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo ) 00695 { 00696 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00697 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00698 { 00699 /* >>chng 03 feb 14, change test from !=0 to >0 */ 00700 /* >>chng 05 feb 07, use lgH2_line_exists */ 00701 /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )*/ 00702 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00703 { 00704 ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul > 0. ); 00705 00706 sumpop += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop; 00707 sumpopA += H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop* 00708 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul; 00709 } 00710 } 00711 } 00712 } 00713 } 00714 } 00715 fprintf(ioQQQ,"DEBUG sumpop = %.3e sumpopA= %.3e A=%.3e\n", sumpop, sumpopA, 00716 sumpopA/SDIV(sumpop) ); 00717 } 00718 } 00719 return; 00720 } 00721 00722 /* increment optical depth for the H2 molecule, called from RT_tau_inc which is called by cloudy, 00723 * one time per zone */ 00724 void H2_RT_tau_inc(void) 00725 { 00726 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00727 00728 /* >>chng 05 jan 26, now use LTE populations for small H2 abundance case, since electronic 00729 * lines become self-shielding surprisingly quickly */ 00730 if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ ) 00731 return; 00732 00733 DEBUG_ENTRY( "H2_RT_tau_inc()" ); 00734 00735 /* remember largest and smallest chemistry renorm factor - 00736 * if both networks are parallel will be unity, 00737 * but only do this after we have stable solution */ 00738 if( nzone > 0 && nCallH2_this_iteration>2 ) 00739 { 00740 h2.renorm_max = MAX2( H2_renorm_chemistry , h2.renorm_max ); 00741 h2.renorm_min = MIN2( H2_renorm_chemistry , h2.renorm_min ); 00742 } 00743 00744 /* loop over all possible lines */ 00745 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00746 { 00747 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00748 { 00749 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00750 { 00751 /* now the lower levels */ 00752 /* NB - X is the only lower level considered here, since we are only 00753 * concerned with excited electronic levels as a photodissociation process 00754 * code exists to relax this assumption - simply change following to iElecHi */ 00755 long int lim_elec_lo = 0; 00756 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00757 { 00758 /* want to include all vibration states in lower level if different electronic level, 00759 * but only lower vibration levels if same electronic level */ 00760 long int nv = h2.nVib_hi[iElecLo]; 00761 if( iElecLo==iElecHi ) 00762 nv = iVibHi; 00763 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00764 { 00765 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00766 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00767 nr = iRotHi-1; 00768 00769 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 00770 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00771 { 00772 /* >>chng 03 feb 14, from !=0 to >0 */ 00773 /* >>chng 05 feb 07, use lgH2_line_exists */ 00774 if( *lgH2le++ ) 00775 { 00776 ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 ); 00777 RT_line_one_tauinc( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , 00778 -9 , iRotHi , iVibLo , iRotLo); 00779 } 00780 }/* iRotLo loop */ 00781 }/* iVibLo loop */ 00782 }/* iElecLo loop */ 00783 }/* iRotHi loop */ 00784 }/* iVibHi loop */ 00785 }/* iElecHi loop */ 00786 /*fprintf(ioQQQ,"\t%.3e\n",H2Lines[1][0][0][0][0][1].TauCon);*/ 00787 return; 00788 } 00789 00790 00791 /* initialize optical depths in H2, called from RT_tau_init */ 00792 void H2_LineZero( void ) 00793 { 00794 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00795 00796 if( !h2.lgH2ON ) 00797 return; 00798 00799 DEBUG_ENTRY( "H2_LineZero()" ); 00800 00801 /* loop over all possible lines */ 00802 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00803 { 00804 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00805 { 00806 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00807 { 00808 /* now the lower levels */ 00809 /* NB - X is the only lower level considered here, since we are only 00810 * concerned with excited electronic levels as a photodissociation process 00811 * code exists to relax this assumption - simply change following to iElecHi */ 00812 long int lim_elec_lo = 0; 00813 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00814 { 00815 /* want to include all vibration states in lower level if different electronic level, 00816 * but only lower vibration levels if same electronic level */ 00817 long int nv = h2.nVib_hi[iElecLo]; 00818 if( iElecLo==iElecHi ) 00819 nv = iVibHi; 00820 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00821 { 00822 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00823 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00824 nr = iRotHi-1; 00825 00826 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00827 { 00828 /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul != 0. )*/ 00829 /* >>chng 03 feb 14, from !=0 to > 0 */ 00830 /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )*/ 00831 /* >>chng 05 feb 07, use lgH2_line_exists */ 00832 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00833 { 00834 /* >>chng 03 feb 14, use TransitionZero rather than explicit sets */ 00835 TransitionZero( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ); 00836 } 00837 } 00838 } 00839 } 00840 } 00841 } 00842 } 00843 return; 00844 } 00845 00846 /* the large H2 molecule, called from RT_tau_reset */ 00847 void H2_RT_tau_reset( void ) 00848 { 00849 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 00850 00851 if( !h2.lgH2ON ) 00852 return; 00853 00854 DEBUG_ENTRY( "H2_RT_tau_reset()" ); 00855 00856 /* loop over all possible lines */ 00857 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00858 { 00859 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00860 { 00861 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00862 { 00863 /* now the lower levels */ 00864 /* NB - X is the only lower level considered here, since we are only 00865 * concerned with excited electronic levels as a photodissociation process 00866 * code exists to relax this assumption - simply change following to iElecHi */ 00867 long int lim_elec_lo = 0; 00868 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00869 { 00870 /* want to include all vibration states in lower level if different electronic level, 00871 * but only lower vibration levels if same electronic level */ 00872 long int nv = h2.nVib_hi[iElecLo]; 00873 if( iElecLo==iElecHi ) 00874 nv = iVibHi; 00875 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00876 { 00877 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00878 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00879 nr = iRotHi-1; 00880 00881 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00882 { 00883 /* >>chng 03 feb 14, change test from !=0 to >0 */ 00884 /*if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul > 0. )*/ 00885 /* >>chng 05 feb 07, use lgH2_line_exists */ 00886 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00887 { 00888 /* inward optical depth */ 00889 RT_line_one_tau_reset( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] , 0.5); 00890 } 00891 } 00892 } 00893 } 00894 } 00895 } 00896 } 00897 return; 00898 } 00899 00900 /* this is fraction of population that is within levels done with matrix */ 00901 static double frac_matrix; 00902 00903 /*H2_Level_low_matrix evaluate lower populations within X */ 00904 STATIC void H2_Level_low_matrix( 00905 /* total abundance within matrix */ 00906 realnum abundance ) 00907 { 00908 00909 /* will need to MALLOC space for these but only on first call */ 00910 static double **AulEscp , 00911 **col_str , 00912 **AulDest, 00913 /* AulPump[low][high] is rate (s^-1) from lower to upper level */ 00914 **AulPump, 00915 **CollRate_levn, 00916 *pops, 00917 *create, 00918 *destroy, 00919 *depart, 00920 /* statistical weight */ 00921 *stat_levn , 00922 /* excitation energies in kelvin */ 00923 *excit; 00924 bool lgDoAs; 00925 static long int levelAsEval=-1; 00926 long int ip; 00927 static bool lgFirst=true; 00928 long int i, 00929 j, 00930 ilo , 00931 ihi, 00932 iElec, 00933 iElecHi, 00934 iVib, 00935 iRot, 00936 iVibHi, 00937 iRotHi; 00938 int nNegPop; 00939 bool lgDeBug, 00940 lgZeroPop; 00941 double rot_cooling , dCoolDT; 00942 static long int ndimMalloced = 0; 00943 double rateout , ratein; 00944 00945 DEBUG_ENTRY( "H2_Level_low_matrix()" ); 00946 00947 /* option to not use the matrix */ 00948 if( nXLevelsMatrix <= 1 ) 00949 { 00950 return; 00951 } 00952 00953 if( lgFirst ) 00954 { 00955 /* check that not more levels than there are in X */ 00956 if( nXLevelsMatrix > nLevels_per_elec[0] ) 00957 { 00958 /* number is greater than number of levels within X */ 00959 fprintf( ioQQQ, 00960 " The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n", 00961 nLevels_per_elec[0]); 00962 cdEXIT(EXIT_FAILURE); 00963 } 00964 /* will never do this again */ 00965 lgFirst = false; 00966 /* remember how much space we malloced in case ever called with more needed */ 00967 /* >>chng 05 jan 19, allocate max number of levels 00968 ndimMalloced = nXLevelsMatrix;*/ 00969 ndimMalloced = nLevels_per_elec[0]; 00970 /* allocate the 1D arrays*/ 00971 excit = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) ); 00972 stat_levn = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) ); 00973 pops = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) ); 00974 create = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) ); 00975 destroy = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) ); 00976 depart = (double *)MALLOC( sizeof(double)*(size_t)(ndimMalloced) ); 00977 /* create space for the 2D arrays */ 00978 AulPump = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *))); 00979 CollRate_levn = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *))); 00980 AulDest = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *))); 00981 AulEscp = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *))); 00982 col_str = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *))); 00983 for( i=0; i<(ndimMalloced); ++i ) 00984 { 00985 AulPump[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double ))); 00986 CollRate_levn[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double ))); 00987 AulDest[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double ))); 00988 AulEscp[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double ))); 00989 col_str[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double ))); 00990 } 00991 00992 for( j=0; j < ndimMalloced; j++ ) 00993 { 00994 stat_levn[j]=0; 00995 excit[j] =0; 00996 } 00997 /* the statistical weights of the levels 00998 * and excitation potentials of each level relative to ground */ 00999 for( j=0; j < ndimMalloced; j++ ) 01000 { 01001 /* obtain the proper indices for the upper level */ 01002 ip = H2_ipX_ener_sort[j]; 01003 iVib = ipVib_H2_energy_sort[ip]; 01004 iRot = ipRot_H2_energy_sort[ip]; 01005 01006 /* statistical weights for each level */ 01007 stat_levn[j] = H2_stat[0][iVib][iRot]; 01008 /* excitation energy of each level relative to ground, in K */ 01009 excit[j] = energy_wn[0][iVib][iRot]*T1CM; 01010 } 01011 01012 for( j=0; j < ndimMalloced-1; j++ ) 01013 { 01014 /* make sure that the energies are ok */ 01015 ASSERT( excit[j+1] > excit[j] ); 01016 } 01017 } 01018 /* end malloc space and creating constant terms */ 01019 01020 /* this is test for call with too many rotation levels to handle - 01021 * logic needs for largest model atom to be called first */ 01022 if( nXLevelsMatrix > ndimMalloced ) 01023 { 01024 fprintf(ioQQQ," H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n"); 01025 cdEXIT(EXIT_FAILURE); 01026 } 01027 01028 /* all elements are used, and must be set to zero */ 01029 for( i=0; i < nXLevelsMatrix; i++ ) 01030 { 01031 pops[i] = 0.; 01032 depart[i] = 0; 01033 for( j=0; j < nXLevelsMatrix; j++ ) 01034 { 01035 col_str[j][i] = 0.; 01036 } 01037 } 01038 01039 /* do we need to reevaluate radiative quantities? only do this one time per zone */ 01040 if( nzone!=nzoneAsEval || iteration!=iterationAsEval || nXLevelsMatrix!=levelAsEval) 01041 { 01042 lgDoAs = true; 01043 nzoneAsEval = nzone; 01044 iterationAsEval = iteration; 01045 levelAsEval = nXLevelsMatrix; 01046 ASSERT( levelAsEval <= ndimMalloced ); 01047 } 01048 else 01049 { 01050 lgDoAs = false; 01051 } 01052 01053 /* all elements are used, and must be set to zero */ 01054 if( lgDoAs ) 01055 { 01056 for( i=0; i < nXLevelsMatrix; i++ ) 01057 { 01058 pops[i] = 0.; 01059 depart[i] = 0; 01060 for( j=0; j < nXLevelsMatrix; j++ ) 01061 { 01062 AulEscp[j][i] = 0.; 01063 AulDest[j][i] = 0.; 01064 AulPump[j][i] = 0.; 01065 CollRate_levn[j][i] = 0.; 01066 } 01067 } 01068 } 01069 01070 /* find all radiative interactions within matrix, and between 01071 * matrix and upper X and excited electronic states */ 01072 iElec = 0; 01073 for( ilo=0; ilo < nXLevelsMatrix; ilo++ ) 01074 { 01075 ip = H2_ipX_ener_sort[ilo]; 01076 iRot = ipRot_H2_energy_sort[ip]; 01077 iVib = ipVib_H2_energy_sort[ip]; 01078 01079 /* H2_X_sink[ilo] includes all processes that destroy H2 in one step, 01080 * these include cosmic ray ionization and dissociation, photodissociation, 01081 * BUT NOT THE SOLOMON process, which, directly, only goes to excited 01082 * electronic states */ 01083 destroy[ilo] = H2_X_sink[ilo]; 01084 01085 /* rates H2 is created from grains and H- units cm-3 s-1, evaluated in mole_H2_form */ 01086 create[ilo] = H2_X_source[ilo]; 01087 01088 /* this loop does radiative decays from upper states inside matrix, 01089 * and upward pumps within matrix region into this lower level */ 01090 if( lgDoAs ) 01091 { 01092 for( ihi=ilo+1; ihi<nXLevelsMatrix; ++ihi ) 01093 { 01094 ip = H2_ipX_ener_sort[ihi]; 01095 iRotHi = ipRot_H2_energy_sort[ip]; 01096 iVibHi = ipVib_H2_energy_sort[ip]; 01097 ASSERT( H2_energies[ip] <= H2_energies[H2_ipX_ener_sort[nXLevelsMatrix-1]] ); 01098 /* general case - but line may not actually exist */ 01099 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) ) 01100 { 01101 /* >>chng 05 feb 07, use lgH2_line_exists */ 01102 if( lgH2_line_exists[0][iVibHi][iRotHi][0][iVib][iRot] ) 01103 { 01104 ASSERT( H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].ipCont > 0 ); 01105 01106 /* NB - the destruction probability is included in 01107 * the total and the destruction is set to zero 01108 * since we want to only count one ots rate, in 01109 * main calling routine, and do not want matrix 01110 * solver below to include it */ 01111 AulEscp[ihi][ilo] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Aul*( 01112 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pesc + 01113 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pdest + 01114 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pelec_esc); 01115 AulDest[ilo][ihi] = 0.; 01116 AulPump[ilo][ihi] = H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->pump; 01117 } 01118 } 01119 } 01120 } 01121 01122 iElecHi = 0; 01123 iElec = 0; 01124 rateout = 0.; 01125 ratein = 0.; 01126 /* now do all levels within X, which are above nXLevelsMatrix, 01127 * the highest level inside the matrix */ 01128 for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi ) 01129 { 01130 ip = H2_ipX_ener_sort[ihi]; 01131 iRotHi = ipRot_H2_energy_sort[ip]; 01132 iVibHi = ipVib_H2_energy_sort[ip]; 01133 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) ) 01134 { 01135 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot] ) 01136 { 01137 ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].ipCont > 0 ); 01138 01139 /* these will enter as net creation terms in creation vector, with 01140 * units cm-3 s-1 01141 * radiative transitions from above the matrix within X */ 01142 ratein += 01143 H2_populations[iElecHi][iVibHi][iRotHi] * 01144 (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Aul* 01145 (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pesc + 01146 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pelec_esc + 01147 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pdest)+H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->pump * 01148 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Lo->g/ 01149 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Hi->g); 01150 /* rate out has units s-1 - destroys current lower level */ 01151 rateout += 01152 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->pump; 01153 } 01154 } 01155 } 01156 01157 /* all states above the matrix but within X */ 01158 create[ilo] += ratein; 01159 01160 /* rates out of matrix into levels in X but above matrix */ 01161 destroy[ilo] += rateout; 01162 01163 /* Solomon process, this sum dos all pump and decays from all electronic excited states */ 01164 /* radiative rates [cm-3 s-1] from electronic excited states into X only vibration and rot */ 01165 create[ilo] += H2_X_rate_from_elec_excited[iVib][iRot]; 01166 01167 /* radiative & cosmic ray rates [s-1] to electronic excited states from X only vibration and rot */ 01168 destroy[ilo] += H2_X_rate_to_elec_excited[iVib][iRot]; 01169 } 01170 01171 /* this flag set with atom H2 trace matrix */ 01172 if( mole.nH2_TRACE >= mole.nH2_trace_matrix ) 01173 lgDeBug = true; 01174 else 01175 lgDeBug = false; 01176 01177 /* now evaluate the rates for all transitions within matrix */ 01178 for( ilo=0; ilo < nXLevelsMatrix; ilo++ ) 01179 { 01180 ip = H2_ipX_ener_sort[ilo]; 01181 iRot = ipRot_H2_energy_sort[ip]; 01182 iVib = ipVib_H2_energy_sort[ip]; 01183 if( lgDoAs ) 01184 { 01185 if(lgDeBug)fprintf(ioQQQ,"DEBUG H2_Level_low_matrix, ilo=%li",ilo); 01186 for( ihi=ilo+1; ihi < nXLevelsMatrix; ihi++ ) 01187 { 01188 ip = H2_ipX_ener_sort[ihi]; 01189 iRotHi = ipRot_H2_energy_sort[ip]; 01190 iVibHi = ipVib_H2_energy_sort[ip]; 01191 /* >>chng 05 may 31, replace with simple expresion */ 01192 CollRate_levn[ihi][ilo] = H2_X_coll_rate[ihi][ilo]; 01193 01194 /*create[ilo] +=CollRate_levn[ihi][ilo]*H2_populations[0][iVibHi][iRotHi];*/ 01195 if(lgDeBug)fprintf(ioQQQ,"\t%.1e",CollRate_levn[ihi][ilo]); 01196 01197 /* now get upward excitation rate - units s-1 */ 01198 CollRate_levn[ilo][ihi] = CollRate_levn[ihi][ilo]* 01199 H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])* 01200 H2_stat[0][iVibHi][iRotHi] / 01201 H2_stat[0][iVib][iRot]; 01202 } 01203 } 01204 01205 if(lgDeBug)fprintf(ioQQQ,"\n"); 01206 01207 /* now do all collisions for levels within X, which are above nXLevelsMatrix, 01208 * the highest level inside the matrix */ 01209 iElecHi = 0; 01210 01211 for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi ) 01212 { 01213 ip = H2_ipX_ener_sort[ihi]; 01214 iRotHi = ipRot_H2_energy_sort[ip]; 01215 iVibHi = ipVib_H2_energy_sort[ip]; 01216 rateout = 0; 01217 /* first do downward deexcitation rate */ 01218 /* >>chng 04 sep 14, do all levels */ 01219 /* >>chng 05 may 31, use summed rate */ 01220 ratein = H2_X_coll_rate[ihi][ilo]; 01221 if(lgDeBug)fprintf(ioQQQ,"\t%.1e",ratein); 01222 01223 /* now get upward excitation rate */ 01224 rateout = ratein * 01225 H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])* 01226 H2_stat[0][iVibHi][iRotHi]/H2_stat[0][iVib][iRot]; 01227 01228 /* these are general entries and exits going into vector */ 01229 create[ilo] += ratein*H2_populations[iElecHi][iVibHi][iRotHi]; 01230 destroy[ilo] += rateout; 01231 } 01232 } 01233 01234 /* H2 grain interactions 01235 * >>chng 05 apr 30,GS, Instead of hmi.H2_total, the specific populations are used because high levels have much less 01236 * populations than ground levels which consists most of the H2 population.*/ 01237 if( lgDoAs ) 01238 { 01239 for( ihi=2; ihi < nXLevelsMatrix; ihi++ ) 01240 { 01241 01242 ip = H2_ipX_ener_sort[ihi]; 01243 iVibHi = ipVib_H2_energy_sort[ip]; 01244 iRotHi = ipRot_H2_energy_sort[ip]; 01245 01246 /* collisions with grains goes to either J=1 or J=0 depending on 01247 * spin of upper level - this conserves op ratio - following 01248 * var is 1 if ortho, 0 if para, so this conserves op ratio 01249 * units are s-1 */ 01250 CollRate_levn[ihi][H2_lgOrtho[0][iVibHi][iRotHi]] += hmi.rate_grain_h2_op_conserve; 01251 } 01252 01253 /* H2 ortho - para conversion on grain surface, 01254 * rate (s-1) all v,J levels go to 0 or 1 */ 01255 CollRate_levn[1][0] += 01256 (realnum)(hmi.rate_grain_h2_J1_to_J0); 01257 } 01258 01259 /* now all levels in X above the matrix */ 01260 for( ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi ) 01261 { 01262 ip = H2_ipX_ener_sort[ihi]; 01263 iVibHi = ipVib_H2_energy_sort[ip]; 01264 iRotHi = ipRot_H2_energy_sort[ip]; 01265 01266 /* these collisions all go into 0 or 1 depending on whether upper level was ortho or para 01267 * units are cm-3 s-1 - rate new molecules appear in matrix */ 01268 create[H2_lgOrtho[0][iVibHi][iRotHi]] += H2_populations[0][iVibHi][iRotHi]*hmi.rate_grain_h2_op_conserve; 01269 } 01270 01271 /* debug print individual contributors to matrix elements */ 01272 { 01273 enum {DEBUG_LOC=false}; 01274 if( DEBUG_LOC || lgDeBug) 01275 { 01276 fprintf(ioQQQ,"DEBUG H2 matexcit"); 01277 for(ilo=0; ilo<nXLevelsMatrix; ++ilo ) 01278 { 01279 fprintf(ioQQQ,"\t%li",ilo ); 01280 } 01281 fprintf(ioQQQ,"\n"); 01282 for(ihi=0; ihi<nXLevelsMatrix;++ihi) 01283 { 01284 fprintf(ioQQQ,"\t%.2e",excit[ihi] ); 01285 } 01286 fprintf(ioQQQ,"\n"); 01287 for(ihi=0; ihi<nXLevelsMatrix;++ihi) 01288 { 01289 fprintf(ioQQQ,"\t%.2e",stat_levn[ihi] ); 01290 } 01291 fprintf(ioQQQ,"\n"); 01292 01293 fprintf(ioQQQ,"AulEscp[n][]\\[][n] = Aul*Pesc\n"); 01294 for(ilo=0; ilo<nXLevelsMatrix; ++ilo ) 01295 { 01296 fprintf(ioQQQ,"\t%li",ilo ); 01297 } 01298 fprintf(ioQQQ,"\n"); 01299 for(ihi=0; ihi<nXLevelsMatrix;++ihi) 01300 { 01301 fprintf(ioQQQ,"%li", ihi); 01302 for(ilo=0; ilo<nXLevelsMatrix; ++ilo ) 01303 { 01304 fprintf(ioQQQ,"\t%.2e",AulEscp[ilo][ihi] ); 01305 } 01306 fprintf(ioQQQ,"\n"); 01307 } 01308 01309 fprintf(ioQQQ,"AulPump [n][]\\[][n]\n"); 01310 for(ilo=0; ilo<nXLevelsMatrix; ++ilo ) 01311 { 01312 fprintf(ioQQQ,"\t%li",ilo ); 01313 } 01314 fprintf(ioQQQ,"\n"); 01315 for(ihi=0; ihi<nXLevelsMatrix;++ihi) 01316 { 01317 fprintf(ioQQQ,"%li", ihi); 01318 for(ilo=0; ilo<nXLevelsMatrix; ++ilo ) 01319 { 01320 fprintf(ioQQQ,"\t%.2e",AulPump[ihi][ilo] ); 01321 } 01322 fprintf(ioQQQ,"\n"); 01323 } 01324 01325 fprintf(ioQQQ,"CollRate_levn [n][]\\[][n]\n"); 01326 for(ilo=0; ilo<nXLevelsMatrix; ++ilo ) 01327 { 01328 fprintf(ioQQQ,"\t%li",ilo ); 01329 } 01330 fprintf(ioQQQ,"\n"); 01331 for(ihi=0; ihi<nXLevelsMatrix;++ihi) 01332 { 01333 fprintf(ioQQQ,"%li", ihi); 01334 for(ilo=0; ilo<nXLevelsMatrix; ++ilo ) 01335 { 01336 fprintf(ioQQQ,"\t%.2e",CollRate_levn[ihi][ilo] ); 01337 } 01338 fprintf(ioQQQ,"\n"); 01339 } 01340 fprintf(ioQQQ,"SOURCE"); 01341 for(ihi=0; ihi<nXLevelsMatrix;++ihi) 01342 { 01343 fprintf(ioQQQ,"\t%.2e",create[ihi]); 01344 } 01345 fprintf(ioQQQ,"\nSINK"); 01346 for(ihi=0; ihi<nXLevelsMatrix;++ihi) 01347 { 01348 fprintf(ioQQQ,"\t%.2e",destroy[ihi]); 01349 } 01350 fprintf(ioQQQ,"\n"); 01351 } 01352 } 01353 01354 atom_levelN( 01355 /* number of levels */ 01356 nXLevelsMatrix, 01357 abundance, 01358 stat_levn, 01359 excit, 01360 'K', 01361 pops, 01362 depart, 01363 /* net transition rate, A * escape prob, s-1, indices are [upper][lower] */ 01364 &AulEscp, 01365 /* col str from high to low */ 01366 &col_str, 01367 &AulDest, 01368 &AulPump, 01369 &CollRate_levn, 01370 create, 01371 destroy, 01372 /* say that we have evaluated the collision rates already */ 01373 true, 01374 &rot_cooling, 01375 &dCoolDT, 01376 " H2 ", 01377 /* nNegPop positive if negative pops occurred, negative if too cold */ 01378 &nNegPop, 01379 &lgZeroPop, 01380 lgDeBug );/* option to print stuff - set to true for debug printout */ 01381 01382 for( i=0; i< nXLevelsMatrix; ++i ) 01383 { 01384 ip = H2_ipX_ener_sort[i]; 01385 iRot = ipRot_H2_energy_sort[ip]; 01386 iVib = ipVib_H2_energy_sort[ip]; 01387 /* >>chng 05 feb 08, do not update h2_old_populations here, since not done 01388 * like this anywhere else - only update H2_populations here, and let single loop 01389 * in main calling routine handle updating various forms of the population 01390 H2_old_populations[0][iVib][iRot] = H2_populations[0][iVib][iRot]; */ 01391 H2_populations[0][iVib][iRot] = pops[i]; 01392 } 01393 01394 if( 0 && mole.nH2_TRACE >= mole.nH2_trace_full) 01395 { 01396 /*static int nn=0; ++nn; if( nn>5)cdEXIT(1);*/ 01397 /* print pops that came out of matrix */ 01398 fprintf(ioQQQ,"\n DEBUG H2_Level_lowJ hmi.H2_total: %.3e matrix rel pops\n",hmi.H2_total); 01399 fprintf(ioQQQ,"v\tJ\tpop\n"); 01400 for( i=0; i<nXLevelsMatrix; ++i ) 01401 { 01402 ip = H2_ipX_ener_sort[i]; 01403 iRot = ipRot_H2_energy_sort[ip]; 01404 iVib = ipVib_H2_energy_sort[ip]; 01405 fprintf(ioQQQ,"%3li\t%3li\t%.3e\t%.3e\t%.3e\n", 01406 iVib , iRot , H2_populations[0][iVib][iRot]/hmi.H2_total , create[i] , destroy[i]); 01407 } 01408 } 01409 01410 /* nNegPop positive if negative pops occurred, negative if too cold */ 01411 if( nNegPop > 0 ) 01412 { 01413 fprintf(ioQQQ," H2_Level_low_matrix called atom_levelN which returned negative H2_populations.\n"); 01414 ConvFail( "pops" , "H2" ); 01415 } 01416 return; 01417 } 01418 /* do level H2_populations for H2, called by Hydrogenic after ionization and H chemistry 01419 * has been recomputed */ 01420 void H2_LevelPops( void ) 01421 { 01422 static double TeUsedColl=-1.f; 01423 double H2_renorm_conserve=0., 01424 H2_renorm_conserve_init=0. , 01425 sumold, 01426 H2_BigH2_H2s, 01427 H2_BigH2_H2g; 01428 double old_solomon_rate=-1.; 01429 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 01430 long int i; 01431 long int n_pop_oscil = 0; 01432 int kase=0; 01433 bool lgConv_h2_soln, 01434 lgPopsConv_total, 01435 lgPopsConv_relative, 01436 lgHeatConv, 01437 lgSolomonConv, 01438 lgOrthoParaRatioConv; 01439 double quant_old=-1., 01440 quant_new=-1.; 01441 01442 bool lgH2_pops_oscil=false, 01443 lgH2_pops_ever_oscil=false; 01444 long int nEner, 01445 ipHi, ipLo; 01446 long int iElec , iVib , iRot,ip; 01447 double sum_pops_matrix; 01448 realnum collup; 01449 /* old and older ortho - para ratios, used to determine whether soln is converged */ 01450 static double ortho_para_old=0. , ortho_para_older=0. , ortho_para_current=0.; 01451 realnum frac_new_oscil=1.f; 01452 01453 /* keep track of changes in population */ 01454 double PopChgMax_relative=0. , PopChgMaxOld_relative=0., PopChgMax_total=0., PopChgMaxOld_total=0.; 01455 long int iRotMaxChng_relative , iVibMaxChng_relative, 01456 iRotMaxChng_total , iVibMaxChng_total, 01457 nXLevelsMatrix_save; 01458 double popold_relative , popnew_relative , popold_total , popnew_total; 01459 /* reason not converged */ 01460 char chReason[100]; 01461 01462 double flux_accum_photodissoc_BigH2_H2g, flux_accum_photodissoc_BigH2_H2s; 01463 long int ip_H2_level; 01464 01465 /* these are convergence criteria - will be increased during search phase */ 01466 double converge_pops_relative=0.1 , 01467 converge_pops_total=1e-3, 01468 converge_ortho_para=1e-2; 01469 01470 /* H2 not on, so space not allocated and return, 01471 * also return if calculation has been declared a failure */ 01472 if( !h2.lgH2ON || lgAbort ) 01473 return; 01474 01475 DEBUG_ENTRY( "H2_LevelPops()" ); 01476 01477 if(mole.nH2_TRACE >= mole.nH2_trace_full ) 01478 { 01479 fprintf(ioQQQ, 01480 "\n***************H2_LevelPops call %li this iteration, zone is %.2f, H2/H:%.e Te:%e ne:%e\n", 01481 nCallH2_this_iteration, 01482 fnzone, 01483 hmi.H2_total/dense.gas_phase[ipHYDROGEN], 01484 phycon.te, 01485 dense.eden 01486 ); 01487 } 01488 else if( mole.nH2_TRACE >= mole.nH2_trace_final ) 01489 { 01490 static long int nzone_prt=-1; 01491 if( nzone!=nzone_prt ) 01492 { 01493 nzone_prt = nzone; 01494 fprintf(ioQQQ,"DEBUG zone %li H2/H:%.3e Te:%.3e *ne:%.3e n(H2):%.3e\n", 01495 nzone, 01496 hmi.H2_total / dense.gas_phase[ipHYDROGEN], 01497 phycon.te, 01498 dense.eden, 01499 hmi.H2_total ); 01500 } 01501 } 01502 01503 /* evaluate Boltzmann factors and LTE unit population - for trivial abundances 01504 * LTE populations are used in place of full solution */ 01505 mole_H2_LTE(); 01506 01507 /* zero out H2_populations and cooling, and return, if H2 fraction is small 01508 * but, if H2 has ever been done, redo irregardless of abundance - 01509 * if large H2 is ever evaluated then mole.H2_to_H_limit is ignored */ 01510 if( (!hmi.lgBigH2_evaluated && hmi.H2_total/dense.gas_phase[ipHYDROGEN] < mole.H2_to_H_limit ) 01511 || hmi.H2_total < 1e-20 ) 01512 { 01513 /* will not compute the molecule */ 01514 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 01515 fprintf(ioQQQ, 01516 " H2_LevelPops pops too small, not computing, set to LTE and return, H2/H is %.2e and mole.H2_to_H_limit is %.2e.", 01517 hmi.H2_total/dense.gas_phase[ipHYDROGEN] , 01518 mole.H2_to_H_limit); 01519 H2_zero_pops_too_low(); 01520 /* end of zero abundance branch */ 01521 return; 01522 } 01523 01524 /* check whether we need to update the H2_Boltzmann factors, LTE level H2_populations, 01525 * and partition function. LTE level pops normalized by partition function, 01526 * so sum of pops is unity */ 01527 01528 /* say that H2 has been computed, ignore previous limit to abund 01529 * in future - this is to prevent oscillations as model is engaged */ 01530 hmi.lgBigH2_evaluated = true; 01531 /* end loop setting H2_Boltzmann factors, partition function, and LTE H2_populations */ 01532 01533 /* >>chng 05 jun 21, 01534 * during search phase we want to use full matrix - save number of levels so that 01535 * we can restore it */ 01536 nXLevelsMatrix_save = nXLevelsMatrix; 01537 if( conv.lgSearch ) 01538 { 01539 nXLevelsMatrix = nLevels_per_elec[0]; 01540 } 01541 01542 /* 05 oct 27, had only reevaluated collision rates when 5% change in temperature 01543 * caused temp failures in large G0 sims - 01544 * do not check whether we need to update the collision rates but 01545 * reevaluate them always 01546 * >>chng 05 nov 04, above caused a 25% increase in the exec time for constant-T sims 01547 * in test suite- original code had reevaluated if > 0.05 change in T - was too much 01548 * change to 10x smaller, change > 0.005 */ 01549 if( fabs(1. - TeUsedColl / phycon.te ) > 0.005 ) 01550 { 01551 H2_CollidRateEvalAll(); 01552 TeUsedColl = phycon.te; 01553 } 01554 01555 /* set the H2_populations when this is the first call to this routine on 01556 * current iteration- will use LTE H2_populations - populations were set by 01557 * call to mole_H2_LTE before above block */ 01558 if( nCallH2_this_iteration==0 || mole.lgH2_LTE ) 01559 { 01560 /* very first call so use LTE H2_populations */ 01561 if(mole.nH2_TRACE >= mole.nH2_trace_full ) 01562 fprintf(ioQQQ,"H2 1st call - using LTE level pops\n"); 01563 01564 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 01565 { 01566 pops_per_elec[iElec] = 0.; 01567 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 01568 { 01569 pops_per_vib[iElec][iVib] = 0.; 01570 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 01571 { 01572 /* LTE populations are for unit H2 density, so need to multiply 01573 * by total H2 density */ 01574 H2_old_populations[iElec][iVib][iRot] = 01575 (realnum)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total; 01576 H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot]; 01577 } 01578 } 01579 } 01580 /* first guess at ortho and para densities */ 01581 h2.ortho_density = 0.75*hmi.H2_total; 01582 h2.para_density = 0.25*hmi.H2_total; 01583 ortho_para_current = h2.ortho_density / SDIV( h2.para_density ); 01584 ortho_para_older = ortho_para_current; 01585 ortho_para_old = ortho_para_current; 01586 /* this is the fraction of the H2 pops that are within the levels done with a matrix */ 01587 frac_matrix = 1.; 01588 } 01589 01590 /* find sum of all H2_populations in X */ 01591 iElec = 0; 01592 pops_per_elec[0] = 0.; 01593 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 01594 { 01595 pops_per_vib[0][iVib] = 0.; 01596 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 01597 { 01598 pops_per_elec[0] += H2_populations[iElec][iVib][iRot]; 01599 pops_per_vib[0][iVib] += H2_populations[iElec][iVib][iRot]; 01600 } 01601 } 01602 ASSERT( pops_per_elec[0]>SMALLFLOAT ); 01603 /* now renorm the old populations to the correct current H2 density, 01604 * At this point pops_per_elec[0] (pop of X) 01605 * is the result of the previous evaluation of the H2 population, 01606 * following is the ratio of the current chemistry solution H2 to the previous H2 */ 01607 H2_renorm_chemistry = hmi.H2_total/ SDIV(pops_per_elec[0]); 01608 01609 /* >>chng 05 jul 13, TE, 01610 * evaluate ratio of H2g and H2s from chemical network and big molecule model */ 01611 iElec = 0; 01612 hmi.H2_chem_BigH2_H2g = 0.; 01613 hmi.H2_chem_BigH2_H2s = 0.; 01614 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 01615 { 01616 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 01617 { 01618 if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR ) 01619 { 01620 hmi.H2_chem_BigH2_H2s += H2_populations[iElec][iVib][iRot]; 01621 01622 } 01623 else 01624 { 01625 hmi.H2_chem_BigH2_H2g += H2_populations[iElec][iVib][iRot]; 01626 01627 } 01628 } 01629 } 01630 01631 hmi.H2_chem_BigH2_H2g = hmi.Hmolec[ipMH2g]/SDIV(hmi.H2_chem_BigH2_H2g); 01632 hmi.H2_chem_BigH2_H2s = hmi.Hmolec[ipMH2s]/SDIV(hmi.H2_chem_BigH2_H2s); 01633 01634 01635 if(mole.nH2_TRACE >= mole.nH2_trace_full) 01636 fprintf(ioQQQ, 01637 "H2 H2_renorm_chemistry is %.4e, hmi.H2_total is %.4e pops_per_elec[0] is %.4e\n", 01638 H2_renorm_chemistry , 01639 hmi.H2_total, 01640 pops_per_elec[0]); 01641 01642 /* renormalize all level populations for the current chemical solution */ 01643 iElec = 0; 01644 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 01645 { 01646 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 01647 { 01648 H2_populations[iElec][iVib][iRot] *= H2_renorm_chemistry; 01649 H2_old_populations[iElec][iVib][iRot] = H2_populations[iElec][iVib][iRot]; 01650 } 01651 } 01652 ASSERT( fabs(h2.ortho_density+h2.para_density-hmi.H2_total)/hmi.H2_total < 0.001 ); 01653 01654 if(mole.nH2_TRACE >= mole.nH2_trace_full ) 01655 fprintf(ioQQQ, 01656 " H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n", 01657 pops_per_elec[0], 01658 hmi.H2_total); 01659 01660 /* >>chng 05 feb 10, reset convergence criteria if we are in search phase */ 01661 if( conv.lgSearch ) 01662 { 01663 converge_pops_relative *= 2.; /*def is 0.1 */ 01664 converge_pops_total *= 3.; /*def is 1e-3*/ 01665 converge_ortho_para *= 3.; /*def is 1e-2*/ 01666 } 01667 01668 /* update state specific rates in H2_X_formation (cm-3 s-1) that H2 forms from grains and H- */ 01669 mole_H2_form(); 01670 01671 /* evaluate total collision rates */ 01672 H2_X_coll_rate_evaluate(); 01673 01674 /* this flag will say whether H2 H2_populations have converged, 01675 * by comparing old and new values */ 01676 lgConv_h2_soln = false; 01677 /* this will count number of passes around following loop */ 01678 loop_h2_pops = 0; 01679 { 01680 static long int nzoneEval=-1; 01681 if( nzone != nzoneEval ) 01682 { 01683 nzoneEval = nzone; 01684 /* this is number of zones with H2 solution in this iteration */ 01685 ++nH2_zone; 01686 } 01687 } 01688 01689 /* begin - start level population solution 01690 * first do electronic excited states, Lyman, Werner, etc 01691 * using old solution for X 01692 * then do matrix if used, then solve for pops of rest of X 01693 * >>chng 04 apr 06, subtract number of oscillations from limit - don't waste loops 01694 * if solution is unstable */ 01695 while( loop_h2_pops < LIM_H2_POP_LOOP-n_pop_oscil && !lgConv_h2_soln && !mole.lgH2_LTE ) 01696 { 01697 double rate_in , rate_out; 01698 static double old_HeatH2Dexc_BigH2=0., HeatChangeOld=0. , HeatChange=0.; 01699 01700 /* this is number of trips around loop this time */ 01701 ++loop_h2_pops; 01702 /* this is number of times through this loop in entire iteration */ 01703 ++nH2_pops; 01704 01705 /* beginning solution for electronic excited states 01706 * loop over all possible pumping routes to excited electronic states 01707 * to get radiative excitation and dissociation rates */ 01708 hmi.H2_H2g_to_H2s_rate_BigH2 = 0.; 01709 01710 rate_out = 0.; 01711 01712 /* these will store radiative rates between electronic excited states and X */ 01713 iElec = 0; 01714 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 01715 { 01716 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 01717 { 01718 /* radiative rates [cm-3 s-1] from electronic excited states into X vibration and rot */ 01719 H2_X_rate_from_elec_excited[iVib][iRot] = 0.; 01720 /* radiative & cosmic ray rates [s-1] to electronic excited states from X */ 01721 H2_X_rate_to_elec_excited[iVib][iRot] = 0.; 01722 } 01723 } 01724 01725 iElecHi = -INT32_MAX; 01726 /* solve for population of electronic excited states by back-substitution */ 01727 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 01728 { 01729 /* for safety, these loop vars are set to insane value, will crash 01730 * if used without being set */ 01731 iVib = -INT32_MAX; 01732 iRot = -INT32_MAX; 01733 iVibLo = -INT32_MAX; 01734 iRotLo = -INT32_MAX; 01735 iVibHi = -INT32_MAX; 01736 iRotHi = -INT32_MAX; 01737 01738 /* this will be total population in each electronic state */ 01739 pops_per_elec[iElecHi] = 0.; 01740 01741 if(mole.nH2_TRACE >= mole.nH2_trace_full) 01742 fprintf(ioQQQ," Pop(e=%li):",iElecHi); 01743 01744 double factor = ( hmi.lgLeidenCRHack ) ? secondaries.x12tot/6.e-17 : 0.; 01745 01746 /* electronic excited state, loop over all vibration */ 01747 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 01748 { 01749 pops_per_vib[iElecHi][iVibHi] = 0.; 01750 /* ======================= EXCITED ELEC STATE POPS LOOP =====================*/ 01751 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 01752 { 01753 /* Solomon process done here, 01754 * sum of all rates into and out of these upper levels 01755 * all inward rates have units cm-3 s-1 */ 01756 rate_in = 0.; 01757 /* this term is spontaneous dissociation of excited electronic states into 01758 * the X continuum 01759 * all outward rates have units s-1 */ 01760 rate_out = H2_dissprob[iElecHi][iVibHi][iRotHi]; 01761 01762 realnum H2gHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->g; 01763 01764 /* now loop over all levels within X find Solomon rate */ 01765 iElecLo=0; 01766 01767 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo ) 01768 { 01769 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 01770 mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 01771 long nr = h2.nRot_hi[iElecLo][iVibLo]; 01772 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 01773 { 01774 /* consider all radiatively permitted transitions between X and excit states */ 01775 if( *lgH2le++ ) 01776 { 01777 ASSERT( H2L->ipCont > 0 ); 01778 01779 /*>>chng 07 jan 4, GS, add collisional excitation by non-thermal electrons in to 01780 * the Singlet state of H2 (B,C,B',D) from equation 5 of Liu & Dalgarno 1994 ApJ, 428, 769, 01781 * assumed incoming energy and outgoing energies are 20eV and 10eV;Eqn 5 gives downward cross 01782 * section and I multiply it by statistical weights to convert into upward cross section */ 01783 double cr_cross_section = H2L->Coll.cs; 01784 # if 0 01785 /* one-time evaluation of this done in mole_h2_create.cpp */ 01786 cr_cross_section = 01787 pow(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng*1e-8,3)* 01788 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g/ 01789 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g)* 01790 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul* 01791 log(3.)*HPLANCK/(160.f*pow(PI,3)*0.5*1e-8*1.6e-12); 01792 ASSERT( fabs(cr_cross_section-H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cs)/ 01793 SDIV(cr_cross_section) < 1e-7 ); 01794 # endif 01795 01796 /* solve electronic excited state, 01797 * rate lower level in X goes to electronic excited state, s-1 01798 * first term is direct pump, second is cosmic ray excitation */ 01799 double rate_one = H2L->Emis->pump 01800 /*>>chng 05 jun 16, GS, add collisional excitation by non-thermal electrons in to 01801 * the Singlet state of H2 (B,C,B',D) from Dalgarno, Yan, & Liu 1999 ApJs;*/ 01802 /*>>chng 07 jan 4, GS, add collisional excitation by non-thermal electrons in to 01803 * the Singlet state of H2 (B,C,B',D) from equation 5 of 01804 *>>refer H2 cr exc Liu & Dalgarno 1994 ApJ, 428, 769 01805 * in terms of hydrogen ionization cross-section 01806 >>refer H2 elec coll Shemansky, D.E., Hall, D.T.& Ajello, J.M. 1985ApJ, 296, 765 */ 01807 +factor*cr_cross_section; 01808 01809 /* this is a permitted electronic transition, must preserve nuclear spin */ 01810 ASSERT( H2_lgOrtho[iElecHi][iVibHi][iRotHi] == H2_lgOrtho[iElecLo][iVibLo][iRotLo] ); 01811 01812 /* this is the rate [cm-3 s-1] electrons move into the upper level from X */ 01813 rate_in += H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_one; 01814 01815 /* this is total X -> excited electronic state rate, cm-3 s-1 */ 01816 /*rate_tot += H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_one;*/ 01817 01818 /* rate [s-1] from levels within X to electronic excited states, 01819 * includes photoexcitation and cosmic ray excitation */ 01820 H2_X_rate_to_elec_excited[iVibLo][iRotLo] += rate_one; 01821 01822 /* excitation rate for Solomon process - this currently has units 01823 * cm-3 s-1 but will be divided by total pop of X and become s-1 */ 01824 /* >>chng 04 may 25, Gargi Shaw found this bug - had been total sum */ 01825 /* hmi.H2_H2g_to_H2s_rate_BigH2 += rate_one/SDIV(hmi.H2_total);*/ 01826 /* this has unit s-1 and will be used for H2g->H2s */ 01827 rate_one = H2L->Emis->Aul* 01828 /* escape and destruction */ 01829 (H2L->Emis->Pesc + 01830 H2L->Emis->Pelec_esc + 01831 H2L->Emis->Pdest) + 01832 /* induced emission down */ 01833 H2L->Emis->pump * 01834 H2L->Lo->g/H2gHi; 01835 01836 /* this is the rate [s-1] electrons leave the excited electronic upper level 01837 * and decay into X - will be used to get pops of electronic excited states */ 01838 rate_out += rate_one; 01839 01840 ASSERT( rate_in >= 0. && rate_out >= 0. ); 01841 } 01842 ++H2L; 01843 } 01844 } 01845 01846 /* update population [cm-3] of the electronic excited state this only includes 01847 * radiative processes between X and excited electronic states, and cosmic rays - 01848 * thermal collisions are neglected 01849 * X is done below and includes all processes */ 01850 H2_rad_rate_out[iElecHi][iVibHi][iRotHi] = rate_out; 01851 double H2popHi = rate_in / SDIV( rate_out ); 01852 H2_populations[iElecHi][iVibHi][iRotHi] = H2popHi; 01853 if( H2_old_populations[iElecHi][iVibHi][iRotHi]==0. ) 01854 H2_old_populations[iElecHi][iVibHi][iRotHi] = H2popHi; 01855 01856 /* for this population, find rate electronic states decay into X */ 01857 /* now loop over all levels within X find Solomon rate */ 01858 iElecLo=0; 01859 for( iVibLo=0; iVibLo<=h2.nVib_hi[iElecLo]; ++iVibLo ) 01860 { 01861 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 01862 mt6ci H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 01863 01864 long nr = h2.nRot_hi[iElecLo][iVibLo]; 01865 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 01866 { 01867 if( *lgH2le++ ) 01868 { 01869 ASSERT( H2L->ipCont > 0 ); 01870 01871 double rate_one = 01872 /* >>chng 05 may 31, from old to new as in matrix */ 01873 /*H2_old_populations[iElecHi][iVibHi][iRotHi]* */ 01874 H2popHi* 01875 (H2L->Emis->Aul* 01876 /* escape and destruction */ 01877 (H2L->Emis->Pesc + H2L->Emis->Pelec_esc + H2L->Emis->Pdest) + 01878 /* induced emission down */ 01879 H2L->Emis->pump * H2L->Lo->g/H2gHi); 01880 01881 /* radiative rates [cm-3 s-1] from electronic excited states to X */ 01882 H2_X_rate_from_elec_excited[iVibLo][iRotLo] += rate_one; 01883 } 01884 ++H2L; 01885 } 01886 } 01887 01888 ASSERT( H2popHi >= 0. && H2popHi <= hmi.H2_total ); 01889 01890 /* this is total pop in this vibration state */ 01891 pops_per_vib[iElecHi][iVibHi] += H2popHi; 01892 01893 /* ======================= POPS EXCITED ELEC CONVERGE LOOP =====================*/ 01894 }/* end excit electronic state rot pops loop */ 01895 01896 if(mole.nH2_TRACE >= mole.nH2_trace_full) 01897 fprintf(ioQQQ,"\t%.2e",pops_per_vib[iElecHi][iVibHi]/hmi.H2_total); 01898 01899 /* total pop in each electronic state */ 01900 pops_per_elec[iElecHi] += pops_per_vib[iElecHi][iVibHi]; 01901 }/* end excit electronic state all vibration loop */ 01902 01903 /* end excited electronic pops loop */ 01904 if(mole.nH2_TRACE >= mole.nH2_trace_full) 01905 fprintf(ioQQQ,"\n"); 01906 } /* end loop over all electronic excited states */ 01907 /*fprintf(ioQQQ,"DEBUG\t%.3e\t%.3e\n", 01908 H2_X_rate_from_elec_excited[0][0], 01909 pops_per_elec[0] );*/ 01910 /* above set pops of excited electronic levels and found rates between them and X - 01911 * now solve highly excited levels within the X state by back-substitution */ 01912 /* these will do convergence check */ 01913 PopChgMaxOld_relative = PopChgMax_relative; 01914 PopChgMaxOld_total = PopChgMax_total; 01915 PopChgMax_relative = 0.; 01916 PopChgMax_total = 0.; 01917 iElec = 0; 01918 iElecHi = 0; 01919 iRotMaxChng_relative =-1; 01920 iVibMaxChng_relative = -1; 01921 iRotMaxChng_total =-1; 01922 iVibMaxChng_total = -1; 01923 popold_relative = 0.; 01924 popnew_relative = 0.; 01925 popold_total = 0.; 01926 popnew_total = 0.; 01927 01928 /* now evaluate total rates for all levels within X */ 01929 for( nEner=0; nEner<nLevels_per_elec[0]; ++nEner ) 01930 { 01931 /* array of energy sorted indices within X */ 01932 ip = H2_ipX_ener_sort[nEner]; 01933 iVib = ipVib_H2_energy_sort[ip]; 01934 iRot = ipRot_H2_energy_sort[ip]; 01935 01936 realnum H2stat = H2_stat[0][iVib][iRot]; 01937 double H2boltz = H2_Boltzmann[0][iVib][iRot]; 01938 01939 /* these will be total rates into and out of the level */ 01940 double col_rate_in = 0.; 01941 double col_rate_out = 0.; 01942 for( ipLo=0; ipLo<nEner; ++ipLo ) 01943 { 01944 ip = H2_ipX_ener_sort[ipLo]; 01945 iVibLo = ipVib_H2_energy_sort[ip]; 01946 iRotLo = ipRot_H2_energy_sort[ip]; 01947 01948 /* this is rate from this level down to lower level, units s-1 */ 01949 col_rate_out += H2_X_coll_rate[nEner][ipLo]; 01950 /*if( nEner==288 ) fprintf(ioQQQ,"DEBUG %3li %3li %.3e\n", nEner,ipLo,H2_X_coll_rate[nEner][ipLo]);*/ 01951 01952 /* inverse, rate up, cm-3 s-1 */ 01953 collup = (realnum)(H2_old_populations[0][iVibLo][iRotLo] * H2_X_coll_rate[nEner][ipLo] * 01954 H2stat / H2_stat[0][iVibLo][iRotLo] * 01955 H2boltz / SDIV( H2_Boltzmann[0][iVibLo][iRotLo] ) ); 01956 01957 col_rate_in += collup; 01958 } 01959 01960 for( ipHi=nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi ) 01961 { 01962 double colldn; 01963 ip = H2_ipX_ener_sort[ipHi]; 01964 iVibHi = ipVib_H2_energy_sort[ip]; 01965 iRotHi = ipRot_H2_energy_sort[ip]; 01966 01967 /* this is rate from this level up to higher level, units s-1 */ 01968 col_rate_out += H2_X_coll_rate[ipHi][nEner] * 01969 H2_stat[0][iVibHi][iRotHi] / H2stat * 01970 (realnum)(H2_Boltzmann[0][iVibHi][iRotHi] / SDIV( H2boltz ) ); 01971 01972 /* rate down from higher level, cm-3 s-1 */ 01973 colldn = H2_old_populations[0][iVibHi][iRotHi] * H2_X_coll_rate[ipHi][nEner]; 01974 /*if( nEner==288 ) fprintf(ioQQQ,"DEBUG %3li %3li %.3e\n", nEner,ipHi,H2_X_coll_rate[ipHi][nEner] * 01975 H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVib][iRot] * 01976 (realnum)(H2_Boltzmann[0][iVibHi][iRotHi] / 01977 SDIV( H2_Boltzmann[0][iVib][iRot] ) ));*/ 01978 01979 col_rate_in += colldn; 01980 } 01981 01982 H2_col_rate_in[iVib][iRot] = col_rate_in; 01983 H2_col_rate_out[iVib][iRot] = col_rate_out; 01984 } 01985 01986 /* =======================INSIDE X POPULATIONS CONVERGE LOOP =====================*/ 01987 /* begin solving for X by back-substitution 01988 * this is the main loop that determines H2_populations within X 01989 * units of all rates in are cm-3 s-1, all rates out are s-1 01990 * nLevels_per_elec is number of levels within electronic 0 - so nEner is one 01991 * beyond end of array here - but will be decremented at start of loop 01992 * this starts at the highest energy wihtin X and moves down to lower energies */ 01993 nEner = nLevels_per_elec[0]; 01994 while( (--nEner) >= nXLevelsMatrix ) 01995 { 01996 01997 /* array of energy sorted indices within X - we are moving down 01998 * starting from highest level within X */ 01999 ip = H2_ipX_ener_sort[nEner]; 02000 iVib = ipVib_H2_energy_sort[ip]; 02001 iRot = ipRot_H2_energy_sort[ip]; 02002 02003 if( nEner+1 < nLevels_per_elec[0] ) 02004 ASSERT( H2_energies[H2_ipX_ener_sort[nEner]] < H2_energies[H2_ipX_ener_sort[nEner+1]] ); 02005 02006 /* >>chng 05 apr 30,GS, Instead of hmi.H2_total, the specific populations are used because high levels have much less 02007 * populations than ground levels which consists most of the H2 population. 02008 * only do this if working level is not v=0, J=0, 1 */ 02009 if( nEner >1 ) 02010 { 02011 H2_col_rate_out[iVib][iRot] += 02012 /* H2 grain interactions 02013 * rate (s-1) all v,J levels go to 0 or 1 preserving spin */ 02014 (realnum)(hmi.rate_grain_h2_op_conserve); 02015 02016 /* this goes into v=0, and J=0 or 1 depending on whether initial 02017 * state is ortho or para */ 02018 H2_col_rate_in[0][H2_lgOrtho[0][iVib][iRot]] += 02019 /* H2 grain interactions 02020 * rate (cm-3 s-1) all v,J levels go to 0 or 1 preserving spin, 02021 * in above lgOrtho says whether should go to 0 or 1 */ 02022 (realnum)(hmi.rate_grain_h2_op_conserve*H2_old_populations[0][iVib][iRot]); 02023 } 02024 else if( nEner == 1 ) 02025 { 02026 /* this is special J=1 to J=0 collision, which is only fast at 02027 * very low grain temperatures */ 02028 H2_col_rate_out[0][1] += 02029 /* H2 grain interactions 02030 * H2 ortho - para conversion on grain surface, 02031 * rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */ 02032 (realnum)(hmi.rate_grain_h2_J1_to_J0); 02033 02034 H2_col_rate_in[0][0] += 02035 /* H2 grain interactions 02036 * H2 ortho - para conversion on grain surface, 02037 * rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */ 02038 (realnum)(hmi.rate_grain_h2_J1_to_J0 *H2_old_populations[0][0][1]); 02039 } 02040 02041 /* will become rate (cm-3 s-1) other levels have radiative transitions to here */ 02042 H2_rad_rate_in[iVib][iRot] = 0.; 02043 H2_rad_rate_out[0][iVib][iRot] = 0.; 02044 02045 /* the next two account for the Solomon process, 02046 * the first is the sum of decays from electronic excited into X 02047 * second is X going into all excited electronic states 02048 * units cm-3 s-1 */ 02049 H2_rad_rate_in[iVib][iRot] += H2_X_rate_from_elec_excited[iVib][iRot]; 02050 02051 /* radiative & cosmic ray rates [s-1] to electronic excited states from X only vibration and rot */ 02052 H2_rad_rate_out[0][iVib][iRot] += H2_X_rate_to_elec_excited[iVib][iRot]; 02053 02054 /* now sum over states within X which are higher than current state */ 02055 iElecHi = 0; 02056 for( ipHi = nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi ) 02057 { 02058 ip = H2_ipX_ener_sort[ipHi]; 02059 iVibHi = ipVib_H2_energy_sort[ip]; 02060 iRotHi = ipRot_H2_energy_sort[ip]; 02061 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02062 /* the rate we enter this state from more highly excited states within X 02063 * by radiative decays, which have delta J = 0 or 2 */ 02064 /* note test on vibration is needed - iVibHi<iVib, energy order ok and space not allocated */ 02065 /* >>chng 05 feb 07, tried to use use lgH2_line_exists but cant */ 02066 if( ( abs(iRotHi-iRot) ==2 || iRotHi==iRot ) && (iVib <= iVibHi) && 02067 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].ipCont > 0 ) 02068 { 02069 double rateone; 02070 rateone = 02071 H2_old_populations[iElecHi][iVibHi][iRotHi]* 02072 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Aul* 02073 (H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pesc + 02074 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pelec_esc + 02075 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pdest); 02076 ASSERT( rateone >=0 ); 02077 02078 /* units cm-3 s-1 */ 02079 H2_rad_rate_in[iVib][iRot] += rateone; 02080 } 02081 } 02082 02083 /* =======================INSIDE POPULATIONS X CONVERGE LOOP =====================*/ 02084 /* we now have total rate this state is populated from above, now get rate 02085 * this state interacts with levels that are below */ 02086 iElecLo = 0; 02087 for( ipLo = 0; ipLo<nEner; ++ipLo ) 02088 { 02089 ip = H2_ipX_ener_sort[ipLo]; 02090 iVibLo = ipVib_H2_energy_sort[ip]; 02091 iRotLo = ipRot_H2_energy_sort[ip]; 02092 /* radiative interactions between this level and lower levels */ 02093 /* the test on vibration is needed - the energies are ok but the space does not exist */ 02094 /* >>chng 05 feb 07, can't use lgH2_line_exists */ 02095 if( ( abs(iRotLo-iRot) == 2 || iRotLo == iRot ) && (iVibLo <= iVib) && 02096 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].ipCont > 0 ) 02097 { 02098 H2_rad_rate_out[0][iVib][iRot] += 02099 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Aul* 02100 (H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pesc + 02101 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pelec_esc + 02102 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pdest); 02103 } 02104 } 02105 /* =======================INSIDE X POPULATIONS CONVERGE LOOP =====================*/ 02106 02107 /* we now have the total rates into and out of this level, get its population 02108 * units cm-3 */ 02109 H2_populations[iElec][iVib][iRot] = 02110 (H2_col_rate_in[iVib][iRot]+ H2_rad_rate_in[iVib][iRot]+H2_X_source[nEner]) / 02111 SDIV(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[0][iVib][iRot]+H2_X_sink[nEner]); 02112 02113 ASSERT( H2_populations[iElec][iVib][iRot] >= 0. ); 02114 } 02115 /* >>chng 05 may 10, move to following back substitution part within X */ 02116 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02117 /* now do lowest levels H2_populations with matrix, 02118 * these should be collisionally dominated */ 02119 if( nXLevelsMatrix ) 02120 { 02121 H2_Level_low_matrix( 02122 /* the total abundance - frac_matrix is fraction of pop that was in these 02123 * levels the last time this was done */ 02124 hmi.H2_total * (realnum)frac_matrix ); 02125 } 02126 iElecHi = 0; 02127 if(mole.nH2_TRACE >= mole.nH2_trace_full) 02128 { 02129 fprintf(ioQQQ," Rel pop(e=%li)" ,iElecHi); 02130 } 02131 02132 /* find ortho and para densites, sum of pops in each vibration */ 02133 /* this will become total pop is X, which will be renormed to equal hmi.H2_total */ 02134 pops_per_elec[0] = 0.; 02135 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 02136 { 02137 double sumv; 02138 sumv = 0.; 02139 pops_per_vib[0][iVibHi] = 0.; 02140 02141 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 02142 { 02143 pops_per_elec[0] += H2_populations[iElecHi][iVibHi][iRotHi]; 02144 sumv += H2_populations[iElecHi][iVibHi][iRotHi]; 02145 pops_per_vib[0][iVibHi] += H2_populations[iElecHi][iVibHi][iRotHi]; 02146 } 02147 /* print sum of H2_populations in each vibration if trace on */ 02148 if(mole.nH2_TRACE >= mole.nH2_trace_full) 02149 fprintf(ioQQQ,"\t%.2e",sumv/hmi.H2_total); 02150 } 02151 ASSERT( pops_per_elec[0] > SMALLFLOAT ); 02152 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02153 if(mole.nH2_TRACE >= mole.nH2_trace_full) 02154 { 02155 fprintf(ioQQQ,"\n"); 02156 /* print the ground vibration state */ 02157 fprintf(ioQQQ," Rel pop(0,J)"); 02158 iElecHi = 0; 02159 iVibHi = 0; 02160 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 02161 { 02162 fprintf(ioQQQ,"\t%.2e",H2_populations[iElecHi][iVibHi][iRotHi]/hmi.H2_total); 02163 } 02164 fprintf(ioQQQ,"\n"); 02165 } 02166 02167 /* now find population in states done with matrix - this is only used to pass 02168 * to matrix solver */ 02169 iElec = 0; 02170 sum_pops_matrix = 0.; 02171 ip =0; 02172 for( i=0; i<nXLevelsMatrix; ++i ) 02173 { 02174 ip = H2_ipX_ener_sort[i]; 02175 iVib = ipVib_H2_energy_sort[ip]; 02176 iRot = ipRot_H2_energy_sort[ip]; 02177 sum_pops_matrix += H2_populations[iElec][iVib][iRot]; 02178 } 02179 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02180 /* this is self consistent since pops_per_elec[0] came from current soln, 02181 * as did the matrix. pops will be renormalized by results from the chemistry 02182 * a few lines down */ 02183 frac_matrix = sum_pops_matrix / SDIV(pops_per_elec[0]); 02184 02185 /* assuming that all H2 population is in X, this is the 02186 * ratio of H2 that came out of the chemistry network to what we just obtained - 02187 * we need to multiply the pops by renorm to agree with the chemistry, 02188 * this routine does not alter hmi.H2_total, but does change pops_per_elec */ 02189 H2_renorm_conserve = hmi.H2_total/ SDIV(pops_per_elec[0]); 02190 /*pops_per_elec[0] = hmi.H2_total;*/ 02191 02192 /* renormalize H2_populations - the H2_populations were updated by renorm when routine entered, 02193 * before pops determined - but population determinations above do not have a sum rule on total 02194 * population - this renorm is to preserve total population */ 02195 H2_sum_excit_elec_den = 0.; 02196 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 02197 { 02198 pops_per_elec[iElecHi] *= H2_renorm_conserve; 02199 if( iElecHi > 0 ) 02200 H2_sum_excit_elec_den += pops_per_elec[iElecHi]; 02201 02202 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 02203 { 02204 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 02205 { 02206 H2_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve; 02207 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02208 } 02209 } 02210 } 02211 02212 /* this loop first checks for largest changes in populations, to determine whether 02213 * we have converged, then updates the population array with a new value, 02214 * which may be a mean of old and new 02215 * update populations check convergence converged */ 02216 sumold = 0.; 02217 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 02218 { 02219 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 02220 { 02221 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 02222 { 02223 double rel_change; 02224 /* keep track of largest relative change in H2_populations to 02225 * determines convergence */ 02226 if( fabs(H2_populations[iElec][iVib][iRot] - 02227 H2_old_populations[iElec][iVib][iRot])/ 02228 /* on first call some very high J states can have zero pop , 02229 * hence the SDIV, will retain sign for checks on oscilations, 02230 * hence the fabs */ 02231 SDIV(H2_populations[iElec][iVib][iRot]) > fabs(PopChgMax_relative) && 02232 /* >>chng 03 jul 19, this had simply been H2_populations > SMALLFLOAT, 02233 * change to relative pops > 1e-15, spent too much time converging 02234 * levels at pops = 1e-37 */ 02235 /* >>chng 03 dec 27, from rel pop 1e-15 to 1e-6 since converging heating will 02236 * be main convergence criteria check convergence */ 02237 /*H2_populations[iElecHi][iVibHi][iRotHi]/SDIV(hmi.H2_total)>1e-15 )*/ 02238 H2_populations[iElec][iVib][iRot]/SDIV(hmi.H2_total)>1e-6 ) 02239 { 02240 PopChgMax_relative = 02241 (H2_populations[iElec][iVib][iRot] - 02242 H2_old_populations[iElec][iVib][iRot])/ 02243 SDIV(H2_populations[iElec][iVib][iRot]); 02244 iRotMaxChng_relative = iRot; 02245 iVibMaxChng_relative = iVib; 02246 popold_relative = H2_old_populations[iElec][iVib][iRot]; 02247 popnew_relative = H2_populations[iElec][iVib][iRot]; 02248 } 02249 /* >>chng 05 feb 08, add largest rel change in total, this will be converged 02250 * down to higher accuracy than above 02251 * keep track of largest change in H2_populations relative to total H2 to 02252 * determine convergence check convergence */ 02253 rel_change = (H2_populations[iElec][iVib][iRot] - 02254 H2_old_populations[iElec][iVib][iRot])/SDIV(hmi.H2_total); 02255 /* retain sign for checks on oscillations hence the fabs */ 02256 if( fabs(rel_change) > fabs(PopChgMax_total) ) 02257 { 02258 PopChgMax_total = rel_change; 02259 iRotMaxChng_total = iRot; 02260 iVibMaxChng_total = iVib; 02261 popold_total = H2_old_populations[iElec][iVib][iRot]; 02262 popnew_total = H2_populations[iElec][iVib][iRot]; 02263 } 02264 02265 kase = -1; 02266 /* update populations - we used the old populations to update the 02267 * current new populations - will do another iteration if they changed 02268 * by much. here old populations are updated for next sweep through molecule */ 02269 /* pop oscillations have occurred - use small changes */ 02270 /* >>chng 04 may 10, turn this back on - now with min on how small frac new 02271 * can become */ 02272 rel_change = fabs( H2_old_populations[iElec][iVib][iRot] - 02273 H2_populations[iElec][iVib][iRot] )/ 02274 SDIV( H2_populations[iElec][iVib][iRot] ); 02275 02276 /* this branch very large changes, use mean of logs but onlly if both are positive*/ 02277 if( rel_change > 3. && 02278 H2_old_populations[iElec][iVib][iRot]*H2_populations[iElec][iVib][iRot]>0 ) 02279 { 02280 /* large changes or oscillations - take average in the log */ 02281 H2_old_populations[iElec][iVib][iRot] = pow( 10. , 02282 log10(H2_old_populations[iElec][iVib][iRot])/2. + 02283 log10(H2_populations[iElec][iVib][iRot])/2. ); 02284 kase = 2; 02285 } 02286 02287 /* modest change, use means of old and new */ 02288 else if( rel_change> 0.1 ) 02289 { 02290 realnum frac_old=0.25f; 02291 /* large changes or oscillations - take average */ 02292 H2_old_populations[iElec][iVib][iRot] = 02293 frac_old*H2_old_populations[iElec][iVib][iRot] + 02294 (1.f-frac_old)*H2_populations[iElec][iVib][iRot]; 02295 kase = 3; 02296 } 02297 else 02298 { 02299 /* small changes, use new value */ 02300 H2_old_populations[iElec][iVib][iRot] = 02301 H2_populations[iElec][iVib][iRot]; 02302 kase = 4; 02303 } 02304 sumold += H2_old_populations[iElec][iVib][iRot]; 02305 } 02306 } 02307 } 02308 /* will renormalize so that total population is correct */ 02309 H2_renorm_conserve_init = hmi.H2_total/sumold; 02310 02311 /* renormalize H2_populations - the H2_populations were updated by renorm when routine entered, 02312 * before pops determined - but population determinations above do not have a sum rule on total 02313 * population - this renorm is to preserve total population */ 02314 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 02315 { 02316 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 02317 { 02318 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 02319 { 02320 H2_old_populations[iElecHi][iVibHi][iRotHi] *= H2_renorm_conserve_init; 02321 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02322 } 02323 } 02324 } 02325 /* get current ortho-para ratio, will be used as test on convergence */ 02326 iElecHi = 0; 02327 h2.ortho_density = 0.; 02328 h2.para_density = 0.; 02329 H2_den_s = 0.; 02330 H2_den_g = 0.; 02331 02332 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 02333 { 02334 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 02335 { 02336 /* find current population in H2s and H2g */ 02337 if( energy_wn[0][iVibHi][iRotHi]> ENERGY_H2_STAR ) 02338 { 02339 H2_den_s += H2_populations[iElecHi][iVibHi][iRotHi]; 02340 } 02341 else 02342 { 02343 H2_den_g += H2_populations[iElecHi][iVibHi][iRotHi]; 02344 } 02345 if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] ) 02346 { 02347 h2.ortho_density += H2_populations[iElecHi][iVibHi][iRotHi]; 02348 } 02349 else 02350 { 02351 h2.para_density += H2_populations[iElecHi][iVibHi][iRotHi]; 02352 } 02353 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02354 } 02355 } 02356 /* these will be used to determine whether solution has converged */ 02357 ortho_para_older = ortho_para_old; 02358 ortho_para_old = ortho_para_current; 02359 ortho_para_current = h2.ortho_density / SDIV( h2.para_density ); 02360 02361 /* this will be evaluated in call to routine that follows - will check 02362 * whether this has converged */ 02363 old_solomon_rate = hmi.H2_Solomon_dissoc_rate_BigH2_H2g; 02364 02365 /* >>chng 05 jul 24, break code out into separate routine for clarify 02366 * located in mole_h2_etc.c - true says to only do Solomon rate */ 02367 H2_Solomon_rate( ); 02368 02369 /* are changes too large? must decide whether population shave converged, 02370 * will check whether H2_populations themselves have changed by much, 02371 * but also change in heating by collisional deexcitation is stable */ 02372 HeatChangeOld = HeatChange; 02373 HeatChange = old_HeatH2Dexc_BigH2 - hmi.HeatH2Dexc_BigH2; 02374 { 02375 static long int loop_h2_oscil=-1; 02376 /* check whether pops are oscillating, as evidenced by change in 02377 * heating changing sign */ 02378 if( loop_h2_pops>2 && ( 02379 (HeatChangeOld*HeatChange<0. ) || 02380 (PopChgMax_relative*PopChgMaxOld_relative<0. ) ) ) 02381 { 02382 lgH2_pops_oscil = true; 02383 if( loop_h2_pops > 6 ) 02384 { 02385 loop_h2_oscil = loop_h2_pops; 02386 lgH2_pops_ever_oscil = true; 02387 /* make this smaller in attempt to damp out oscillations, 02388 * but don't let get too small*/ 02389 frac_new_oscil *= 0.8f; 02390 frac_new_oscil = MAX2( frac_new_oscil , 0.1f); 02391 ++n_pop_oscil; 02392 } 02393 } 02394 else 02395 { 02396 lgH2_pops_oscil = false; 02397 /* turn off flag if no oscillations for a while */ 02398 if( loop_h2_pops - loop_h2_oscil > 4 ) 02399 { 02400 frac_new_oscil = 1.f; 02401 lgH2_pops_ever_oscil = false; 02402 } 02403 } 02404 } 02405 02406 /* reevaluate heating - cooling if H2 molecule is significant source or either, 02407 * since must have stable heating cooling rate */ 02408 old_HeatH2Dexc_BigH2 = hmi.HeatH2Dexc_BigH2; 02409 if(fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > conv.HeatCoolRelErrorAllowed/10. || 02410 hmi.HeatH2Dexc_BigH2==0. ) 02411 H2_Cooling("H2lup"); 02412 02413 /* begin check on whether solution is converged */ 02414 lgConv_h2_soln = true; 02415 lgPopsConv_total = true; 02416 lgPopsConv_relative = true; 02417 lgHeatConv = true; 02418 lgSolomonConv = true; 02419 lgOrthoParaRatioConv = true; 02420 02421 /* these are all the convergence tests 02422 * check convergence converged */ 02423 if( fabs(PopChgMax_relative)>converge_pops_relative ) 02424 { 02425 /*lgPopsConv = (fabs(PopChgMax_relative)<=0.1);*/ 02426 lgConv_h2_soln = false; 02427 lgPopsConv_relative = false; 02428 /* >>chng 04 sep 08, set quant_new to new chng max gs */ 02429 /*quant_old = PopChgMax_relative;*/ 02430 quant_old = PopChgMaxOld_relative; 02431 /*quant_new = 0.;*/ 02432 quant_new = PopChgMax_relative; 02433 02434 strcpy( chReason , "rel pops changed" ); 02435 } 02436 02437 /* check largest change in a level population relative to total h2 02438 * population convergence converged check */ 02439 else if( fabs(PopChgMax_total)>converge_pops_total) 02440 { 02441 lgConv_h2_soln = false; 02442 lgPopsConv_total = false; 02443 /* >>chng 04 sep 08, set quant_new to new chng max gs */ 02444 /*quant_old = PopChgMax_relative;*/ 02445 quant_old = PopChgMaxOld_total; 02446 /*quant_new = 0.;*/ 02447 quant_new = PopChgMax_total; 02448 02449 strcpy( chReason , "tot pops changed" ); 02450 } 02451 02452 /* >>chng 04 apr 30, look at change in ortho-para ratio, also that is not 02453 * oscillating */ 02454 /* >>chng 04 dec 15, only look at change, and don't make allowed change so tiny - 02455 * these were attempts at fixing problems that were due to shielding not thin*/ 02456 else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> converge_ortho_para ) 02457 /* else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> 1e-3 02458 && (ortho_para_current-ortho_para_old)*(ortho_para_old-ortho_para_older)>0. )*/ 02459 { 02460 lgConv_h2_soln = false; 02461 lgOrthoParaRatioConv = false; 02462 quant_old = ortho_para_old; 02463 quant_new = ortho_para_current; 02464 strcpy( chReason , "ortho/para ratio changed" ); 02465 } 02466 /* >>chng 04 dec 16, reduce error allowed fm /5 to /2, to be similar to 02467 * logic in conv_base */ 02468 else if( !thermal.lgTemperatureConstant && 02469 fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/MAX2(thermal.ctot,thermal.htot) > 02470 conv.HeatCoolRelErrorAllowed/2. 02471 /* >>chng 04 may 09, do not check on error in heating if constant temperature */ 02472 /*&& !(thermal.lgTemperatureConstant || phycon.te <= phycon.TEMP_LIMIT_LOW )*/ ) 02473 { 02474 /* default on HeatCoolRelErrorAllowed is 0.02 */ 02475 /*lgHeatConv = (fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot <= 02476 * conv.HeatCoolRelErrorAllowed/5.);*/ 02477 lgConv_h2_soln = false; 02478 lgHeatConv = false; 02479 quant_old = old_HeatH2Dexc_BigH2/MAX2(thermal.ctot,thermal.htot); 02480 quant_new = hmi.HeatH2Dexc_BigH2/MAX2(thermal.ctot,thermal.htot); 02481 strcpy( chReason , "heating changed" ); 02482 /*fprintf(ioQQQ,"DEBUG old new trip \t%.4e \t %.4e\n", 02483 old_HeatH2Dexc_BigH2, 02484 hmi.HeatH2Dexc_BigH2);*/ 02485 } 02486 02487 /* check on Solomon rate, 02488 * >>chng 04 aug 28, do not do this check if induced processes are disabled, 02489 * since Solomon process is then irrelevant */ 02490 /* >>chng 04 sep 21, GS*/ 02491 else if( rfield.lgInducProcess && 02492 /* this is check that H2 abundance has not been set - if it has been 02493 * then we don't care what the Solomon rate is doing */ 02494 hmi.H2_frac_abund_set==0 && 02495 /*>>chng 05 feb 10, rather than checking change in Solomon relative to Solomon, 02496 * check it relative to total h2 destruction rate */ 02497 fabs( hmi.H2_Solomon_dissoc_rate_BigH2_H2g - old_solomon_rate)/SDIV(hmi.H2_rate_destroy) > 02498 conv.EdenErrorAllowed/5.) 02499 { 02500 lgConv_h2_soln = false; 02501 lgSolomonConv = false; 02502 quant_old = old_solomon_rate; 02503 quant_new = hmi.H2_Solomon_dissoc_rate_BigH2_H2g; 02504 strcpy( chReason , "Solomon rate changed" ); 02505 } 02506 02507 /* did we pass all the convergence test */ 02508 if( !lgConv_h2_soln ) 02509 { 02510 /* this branch H2 H2_populations within X are not converged, 02511 * print diagnostic */ 02512 02513 if( PRT_POPS || mole.nH2_TRACE >=mole.nH2_trace_iterations ) 02514 { 02515 /*fprintf(ioQQQ,"temppp\tnew\t%.4e\tnew\t%.4e\t%.4e\n", 02516 hmi.HeatH2Dexc_BigH2, 02517 old_HeatH2Dexc_BigH2, 02518 fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot );*/ 02519 fprintf(ioQQQ," loop %3li no conv oscl?%c why:%s ", 02520 loop_h2_pops, 02521 TorF(lgH2_pops_ever_oscil), 02522 chReason ); 02523 if( !lgPopsConv_relative ) 02524 fprintf(ioQQQ," PopChgMax_relative:%.4e v:%li J:%li old:%.4e new:%.4e", 02525 PopChgMax_relative, 02526 iVibMaxChng_relative, 02527 iRotMaxChng_relative , 02528 popold_relative , 02529 popnew_relative ); 02530 else if( !lgPopsConv_total ) 02531 fprintf(ioQQQ," PopChgMax_total:%.4e v:%li J:%li old:%.4e new:%.4e", 02532 PopChgMax_total, 02533 iVibMaxChng_total, 02534 iRotMaxChng_total , 02535 popold_total , 02536 popnew_total ); 02537 else if( !lgHeatConv ) 02538 fprintf(ioQQQ," heat:%.4e old:%.4e new:%.4e", 02539 (hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/MAX2(thermal.ctot,thermal.htot), 02540 quant_old , 02541 quant_new); 02542 /* Solomon rate changed */ 02543 else if( !lgSolomonConv ) 02544 fprintf(ioQQQ," d(sol rate)/tot dest\t%2e",(old_solomon_rate - hmi.H2_Solomon_dissoc_rate_BigH2_H2g)/SDIV(hmi.H2_rate_destroy)); 02545 else if( !lgOrthoParaRatioConv ) 02546 fprintf(ioQQQ," current, old, older ratios are %.4e %.4e %.4e", 02547 ortho_para_current , ortho_para_old, ortho_para_older ); 02548 else 02549 TotalInsanity(); 02550 fprintf(ioQQQ,"\n"); 02551 } 02552 } 02553 /* end convergence criteria */ 02554 02555 /*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", 02556 loop_h2_pops, 02557 fnzone, 02558 phycon.te, 02559 dense.eden, 02560 hmi.HeatH2Dexc_BigH2, 02561 hmi.HeatH2Dexc_BigH2/thermal.ctot , 02562 hmi.H2_total, 02563 H2_renorm_chemistry , 02564 H2_renorm_conserve, 02565 hmi.H2_H2g_to_H2s_rate_BigH2);*/ 02566 if( trace.nTrConvg >= 5 ) 02567 { 02568 fprintf( ioQQQ, 02569 " H2 5lev %li Conv?%c", 02570 loop_h2_pops , 02571 TorF(lgConv_h2_soln) ); 02572 02573 if( fabs(PopChgMax_relative)>0.1 ) 02574 fprintf(ioQQQ," pops, rel chng %.3e",PopChgMax_relative); 02575 else 02576 fprintf(ioQQQ," rel heat %.3e rel chng %.3e H2 heat/cool %.2e", 02577 hmi.HeatH2Dexc_BigH2/thermal.ctot , 02578 fabs(hmi.HeatH2Dexc_BigH2-old_HeatH2Dexc_BigH2)/thermal.ctot , 02579 hmi.HeatH2Dexc_BigH2/thermal.ctot); 02580 02581 fprintf( ioQQQ, 02582 " Oscil?%c Ever Oscil?%c", 02583 TorF(lgH2_pops_oscil) , 02584 TorF(lgH2_pops_ever_oscil) ); 02585 if( lgH2_pops_ever_oscil ) 02586 fprintf(ioQQQ," frac_new_oscil %.4f",frac_new_oscil); 02587 fprintf(ioQQQ,"\n"); 02588 } 02589 02590 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 02591 { 02592 fprintf(ioQQQ, 02593 "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", 02594 loop_h2_pops, 02595 kase, 02596 H2_renorm_chemistry, 02597 h2.ortho_density / h2.para_density , 02598 frac_matrix); 02599 02600 /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/ 02601 if( iVibMaxChng_relative>=0 && iRotMaxChng_relative>=0 ) 02602 fprintf(ioQQQ, 02603 "end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n", 02604 loop_h2_pops, 02605 PopChgMax_relative , 02606 H2_old_populations[0][iVibMaxChng_relative][iRotMaxChng_relative], 02607 H2_populations[0][iVibMaxChng_relative][iRotMaxChng_relative], 02608 iVibMaxChng_relative , iRotMaxChng_relative 02609 ); 02610 } 02611 } 02612 /* =======================END POPULATIONS CONVERGE LOOP =====================*/ 02613 02614 /* evaluate H2 rates over H2g and H2s for use in chemistry network */ 02615 H2_gs_rates(); 02616 02617 /* >>chng 05 feb 08, do not print if we are in search phase */ 02618 if( !lgConv_h2_soln && !conv.lgSearch ) 02619 { 02620 conv.lgConvPops = false; 02621 strcpy( conv.chConvIoniz, "H2 pop cnv" ); 02622 fprintf(ioQQQ, 02623 " H2_LevelPops: H2_populations not converged in %li tries; due to %s, old, new are %.4e %.4e, iteration %li zone %.2f.\n", 02624 loop_h2_pops, 02625 chReason, 02626 quant_old, 02627 quant_new , 02628 iteration , 02629 fnzone ); 02630 ConvFail("pops","H2"); 02631 } 02632 02633 /* loop over all possible lines and set H2_populations, 02634 * and quantities that depend on them */ 02635 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 02636 { 02637 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 02638 { 02639 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 02640 { 02641 long int lim_elec_lo = 0; 02642 double H2popHi = H2_populations[iElecHi][iVibHi][iRotHi]; 02643 // this will update Lo->Pop as well as Lo and Hi point to the same set of data structures 02644 // the loops below are OK since Lo->Pop is only accessed after Hi->Pop has been set. 02645 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop = H2popHi; 02646 ASSERT( H2popHi >= 0. ); 02647 realnum H2gHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->g; 02648 /* now the lower levels */ 02649 /* NB - X is the only lower level considered here, since we are only 02650 * concerned with excited electronic levels as a photodissociation process 02651 * code exists to relax this assumption - simply change following to iElecHi */ 02652 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 02653 { 02654 /* want to include all vibration states in lower level if different electronic level, 02655 * but only lower vibration levels if same electronic level */ 02656 long int nv = h2.nVib_hi[iElecLo]; 02657 if( iElecLo==iElecHi ) 02658 nv = iVibHi; 02659 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 02660 { 02661 long nr = h2.nRot_hi[iElecLo][iVibLo]; 02662 if( iElecLo==iElecHi && iVibHi==iVibLo ) 02663 nr = iRotHi-1; 02664 02665 mb6ci lgH2le = lgH2_line_exists.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 02666 mt6i H2L = H2Lines.ptr(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,h2.Jlowest[iElecLo]); 02667 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 02668 { 02669 if( *lgH2le++ ) 02670 { 02671 /* following two heat exchange excitation, deexcitation */ 02672 H2L->Coll.cool = 0.; 02673 H2L->Coll.heat = 0.; 02674 02675 H2L->Emis->PopOpc = 02676 H2L->Lo->Pop - H2popHi * H2L->Lo->g / H2gHi; 02677 02678 /* number of photons in the line */ 02679 H2L->Emis->phots = H2L->Emis->Aul * 02680 (H2L->Emis->Pesc + H2L->Emis->Pelec_esc) * 02681 H2popHi; 02682 02683 /* intensity of line */ 02684 H2L->Emis->xIntensity = 02685 H2L->Emis->phots * 02686 H2L->EnergyErg; 02687 02688 if( iElecHi==0 ) 02689 { 02691 /* the ground electronic state, most excitations are not direct pumping 02692 * (rather indirect, which does not count for ColOvTot) */ 02693 H2L->Emis->ColOvTot = 1.; 02694 } 02695 else 02696 { 02697 /* these are excited electronic states, mostly pumped, except for supras */ 02699 H2L->Emis->ColOvTot = 0.; 02700 } 02701 } 02702 ++H2L; 02703 } 02704 } 02705 } 02706 } 02707 } 02708 } 02709 02710 /* add up H2 + hnu => 2H, continuum photodissociation, 02711 * this is not the Solomon process, true continuum */ 02712 /* >>chng 05 jun 16, GS, add dissociation to triplet states*/ 02713 hmi.H2_photodissoc_BigH2_H2s = 0.; 02714 hmi.H2_photodissoc_BigH2_H2g = 0.; 02715 hmi.H2_tripletdissoc_H2s =0.; 02716 hmi.H2_tripletdissoc_H2g =0.; 02717 hmi.H2_BigH2_H2g_av = 0.; 02718 hmi.H2_BigH2_H2s_av = 0.; 02719 /* >>chng 05 jul 20, GS, add dissociation by H2 g and H2s*/ 02720 hmi.Average_collH2s_dissoc = 0.; 02721 hmi.Average_collH2g_dissoc = 0.; 02722 02723 iElec = 0; 02724 H2_BigH2_H2s = 0.; 02725 H2_BigH2_H2g = 0.; 02726 hmi.H2g_BigH2 =0; 02727 hmi.H2s_BigH2 = 0; 02728 hmi.H2_total_BigH2 =0; 02729 hmi.H2g_LTE_bigH2 =0.; 02730 hmi.H2s_LTE_bigH2 = 0.; 02731 /* >>chng 05 oct 20, no need to reset this var here */ 02732 /*exp_disoc = sexp(H2_DissocEnergies[0]/phycon.te_wn);*/ 02733 02734 /* >>chng 05 sep 12, TE, define a cutoff wavelength of 800 Angstrom 02735 * this is chosen as the cross sections given by 02736 *>>refer H2 photo cs Allison, A.C. & Dalgarno, A. 1969, Atomic Data, 1, 91 02737 * show a sharp decline in the cross section*/ 02738 { 02739 static long ip_cut_off = -1; 02740 if( ip_cut_off < 0 ) 02741 { 02742 /* one-time initialization of this pointer */ 02743 ip_cut_off = ipoint( 1.14 ); 02744 } 02745 02746 /* >>chng 05 sep 12, TE, assume all H2s is at 2.5 eV 02747 * the dissociation threshold is at 1.07896 Rydberg*/ 02748 flux_accum_photodissoc_BigH2_H2s = 0; 02749 ip_H2_level = ipoint( 1.07896 - 2.5 / EVRYD); 02750 for( i= ip_H2_level; i < ip_cut_off; ++i ) 02751 { 02752 flux_accum_photodissoc_BigH2_H2s += ( rfield.flux[i-1] + rfield.ConInterOut[i-1]+ 02753 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1] ); 02754 } 02755 02756 /* sum over all levels to obtain s and g populations and dissociation rates */ 02757 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 02758 { 02759 for( iRot=h2.Jlowest[iElec]; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 02760 { 02761 /* >>chng 05 mar 22, TE, moved H2_photodissoc_BigH2_H2s in this statement and divide by the 02762 * density of H2s not total H2, we consider direct photodissociation only for H2s */ 02763 /* >>chng 05 mar 22, TE, this should be for H2* rather than total */ 02764 /* this is the total rate of direct photo-dissociation of excited electronic states into 02765 * the X continuum - this is continuum photodissociation, not the Solomon process */ 02766 /* >>chng 03 sep 03, make sum of pops of excited states */ 02767 if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR ) 02768 { 02769 double arg_ratio; 02770 hmi.H2_photodissoc_BigH2_H2s += 02771 H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2s; 02772 02773 /* cosmic ray & secondary electron excitation to triplets 02774 * physics described where similar process is done for 02775 * big molecule 02776 *>>chng 07 apr 08, from 3 to 10 to better capture results 02777 * of Dalgarno et al 99 */ 02778 hmi.H2_tripletdissoc_H2s += 02779 H2_populations[iElec][iVib][iRot] * 10.f*secondaries.x12tot; 02780 02781 /* sum of pops in levels in H2* for use in chemistry network */ 02782 H2_BigH2_H2s += H2_populations[iElec][iVib][iRot]; 02783 02784 /* >>chng 05 jun 28, TE, determine average energy level in H2s */ 02785 hmi.H2_BigH2_H2s_av += (H2_populations[iElec][iVib][iRot] * energy_wn[0][iVib][iRot]); 02786 02787 /* >>chng 05 july 20, GS, collisional dissociation by H2s, unit s-1*/ 02788 hmi.Average_collH2s_dissoc += H2_populations[iElec][iVib][iRot] * H2_coll_dissoc_rate_coef_H2[iVib][iRot]; 02789 02790 /* >>chng 05 oct 17, GS, LTE populations of H2s*/ 02791 arg_ratio = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]); 02792 if( arg_ratio > 0. ) 02793 { 02794 /* >>chng 05 oct 21, GS, only add ratio if Boltzmann factor > 0 */ 02795 hmi.H2s_LTE_bigH2 += H2_populations[0][iVib][iRot]*SAHA/SDIV(phycon.te32*arg_ratio)* 02796 (H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5; 02797 } 02798 } 02799 else 02800 { 02801 double arg_ratio; 02802 /* >>chng 05 sep 12, TE, for H2g do the sum explicitly for every level*/ 02803 flux_accum_photodissoc_BigH2_H2g = 0; 02804 /* this is the dissociation energy needed for the level*/ 02805 ip_H2_level = ipoint( 1.07896 - energy_wn[0][iVib][iRot] * WAVNRYD); 02806 02807 for( i= ip_H2_level; i < ip_cut_off; ++i ) 02808 { 02809 flux_accum_photodissoc_BigH2_H2g += ( rfield.flux[i-1] + rfield.ConInterOut[i-1]+ 02810 rfield.outlin[i-1]+ rfield.outlin_noplot[i-1] ); 02811 } 02812 02813 hmi.H2_photodissoc_BigH2_H2g += 02814 H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2g; 02815 02816 02817 /* cosmic ray & secondary electron excitation to triplets 02818 * physics described where similar process is done for 02819 * big molecule 02820 *>>chng 07 apr 08, from 3 to 10 to better capture results 02821 * of Dalgarno et al 99 */ 02822 hmi.H2_tripletdissoc_H2g += 02823 H2_populations[iElec][iVib][iRot] * 10.f*secondaries.x12tot; 02824 02825 /* sum of pops in levels in H2g for use in chemistry network */ 02826 H2_BigH2_H2g += H2_populations[iElec][iVib][iRot]; 02827 02828 /* >>chng 05 jun 28, TE, determine average energy level in H2g */ 02829 hmi.H2_BigH2_H2g_av += (H2_populations[iElec][iVib][iRot] * energy_wn[0][iVib][iRot]); 02830 02831 /* >>chng 05 jul 20, GS, collisional dissociation by H2s, unit s-1*/ 02832 hmi.Average_collH2g_dissoc += H2_populations[iElec][iVib][iRot] * H2_coll_dissoc_rate_coef_H2[iVib][iRot]; 02833 02834 /* >>chng 05 oct 17, GS, LTE populations of H2g*/ 02835 arg_ratio = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]); 02836 if( arg_ratio > 0. ) 02837 { 02838 hmi.H2g_LTE_bigH2 += H2_populations[0][iVib][iRot]*(SAHA/SDIV(phycon.te32*arg_ratio)* 02839 (H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5); 02840 } 02841 } 02842 } 02843 } 02844 } 02845 hmi.H2g_BigH2 = (realnum)H2_BigH2_H2g; 02846 hmi.H2s_BigH2 = (realnum)H2_BigH2_H2s; 02847 hmi.H2_total_BigH2 =hmi.H2g_BigH2+hmi.H2s_BigH2; 02848 02849 /* average energy in H2s */ 02850 hmi.H2_BigH2_H2s_av = hmi.H2_BigH2_H2s_av / H2_BigH2_H2s; 02851 /* average energy in H2g */ 02852 hmi.H2_BigH2_H2g_av = hmi.H2_BigH2_H2g_av / H2_BigH2_H2g; 02853 02854 /* above sum was rate per unit vol since mult by H2 density, now div by H2* density to get rate s-1 */ 02855 /* 0.25e-18 is wild guess of typical photodissociation cross section, from 02856 * >>refer H2 dissoc Allison, A.C. & Dalgarno, A. 1969, Atomic Data, 1, 91 02857 * this is based on an average of the highest v values they gave. unfortunately, we want 02858 * the highest J values - 02859 * final units are s-1*/ 02860 hmi.H2_photodissoc_BigH2_H2s = hmi.H2_photodissoc_BigH2_H2s / SDIV(H2_BigH2_H2s) * H2_DISS_ALLISON_DALGARNO; 02861 hmi.H2_photodissoc_BigH2_H2g = hmi.H2_photodissoc_BigH2_H2g / SDIV(H2_BigH2_H2g) * H2_DISS_ALLISON_DALGARNO; 02862 hmi.H2_tripletdissoc_H2g = hmi.H2_tripletdissoc_H2g/SDIV(H2_BigH2_H2g); 02863 hmi.H2_tripletdissoc_H2s = hmi.H2_tripletdissoc_H2s/SDIV(H2_BigH2_H2s); 02864 hmi.Average_collH2g_dissoc = hmi.Average_collH2g_dissoc /SDIV(H2_BigH2_H2g);/* unit cm3s-1*/ 02865 hmi.Average_collH2s_dissoc = hmi.Average_collH2s_dissoc /SDIV(H2_BigH2_H2s);/* unit cm3s-1*/ 02866 hmi.H2s_LTE_bigH2 = hmi.H2s_LTE_bigH2/SDIV(H2_BigH2_H2s); 02867 hmi.H2g_LTE_bigH2 = hmi.H2g_LTE_bigH2/SDIV(H2_BigH2_H2g); 02868 02869 /* >>chng 05 jul 09, GS*/ 02870 /* average Einstein value for H2* to H2g, GS*/ 02871 double sumpop1 = 0.; 02872 double sumpopA1 = 0.; 02873 double sumpopcollH2O_deexcit = 0.; 02874 double sumpopcollH2p_deexcit = 0.; 02875 double sumpopcollH_deexcit = 0.; 02876 double popH2s = 0.; 02877 double sumpopcollH2O_excit = 0.; 02878 double sumpopcollH2p_excit = 0.; 02879 double sumpopcollH_excit = 0.; 02880 double popH2g = 0.; 02881 02882 02883 iElecLo = 0; 02884 for( iVibHi=0; iVibHi<=h2.nVib_hi[0]; ++iVibHi ) 02885 { 02886 long nr1 = h2.nRot_hi[0][iVibHi]; 02887 for( iRotHi=h2.Jlowest[0]; iRotHi<=nr1; ++iRotHi ) 02888 { 02889 double ewnHi = energy_wn[0][iVibHi][iRotHi]; 02890 for( iVibLo=0; iVibLo<=h2.nVib_hi[0]; ++iVibLo ) 02891 { 02892 long nr2 = h2.nRot_hi[0][iVibLo]; 02893 md3ci ewnLo = energy_wn.ptr(0,iVibLo,h2.Jlowest[0]); 02894 for( iRotLo=h2.Jlowest[0]; iRotLo<=nr2; ++iRotLo ) 02895 { 02896 double ewnLo2 = energy_wn[0][iVibLo][iRotLo]; 02897 /*if( ewnHi > ENERGY_H2_STAR && *ewnLo++ < ENERGY_H2_STAR )*/ 02898 if( ewnHi > ENERGY_H2_STAR && ewnLo2 < ENERGY_H2_STAR ) 02899 { 02900 /* >>chng 05 jul 10, GS*/ 02901 /* average collisional rate for H2* to H2g, GS*/ 02902 if( H2_lgOrtho[0][iVibHi][iRotHi] == H2_lgOrtho[0][iVibLo][iRotLo] ) 02903 { 02904 /* sums of populations */ 02905 popH2s += H2_populations[0][iVibHi][iRotHi]; 02906 popH2g += H2_populations[0][iVibLo][iRotLo]; 02907 02908 /* sums of deexcitation rates - H2* to H2g */ 02909 sumpopcollH_deexcit += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0]; 02910 sumpopcollH2O_deexcit += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][2]; 02911 sumpopcollH2p_deexcit += H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][3]; 02912 02913 /* sums of excitation rates - H2g to H2* */ 02914 sumpopcollH_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] * 02915 H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] ); 02916 sumpopcollH2O_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][2]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] * 02917 H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] ); 02918 sumpopcollH2p_excit += H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][3]*H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] * 02919 H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] ); 02920 02921 02922 02923 /* if( (abs((iRotHi-iRotLo) ))==2 || (iRotHi==iRotLo ) ) */ 02924 if( lgH2_line_exists[0][iVibHi][iRotHi][0][iVibLo][iRotLo] ) 02925 { 02926 sumpop1 += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Hi->Pop; 02927 sumpopA1 += H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Hi->Pop* 02928 H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Aul; 02929 } 02930 } 02931 } 02932 } 02933 } 02934 } 02935 } 02936 hmi.Average_A = sumpopA1/SDIV(sumpop1); 02937 02938 /* collisional excitation and deexcitation of H2g and H2s */ 02939 hmi.Average_collH2_deexcit = (sumpopcollH2O_deexcit+sumpopcollH2p_deexcit)/SDIV(popH2s); 02940 hmi.Average_collH2_excit = (sumpopcollH2O_excit+sumpopcollH2p_excit)/SDIV(popH2g); 02941 hmi.Average_collH_excit = sumpopcollH_excit/SDIV(popH2g); 02942 hmi.Average_collH_deexcit = sumpopcollH_deexcit/SDIV(popH2s); 02943 /* populations are badly off during search phase */ 02944 if( conv.lgSearch ) 02945 { 02946 hmi.Average_collH_excit /=10.; 02947 hmi.Average_collH_deexcit/=10.; 02948 } 02949 02950 /*fprintf(ioQQQ, 02951 "DEBUG Average_collH_excit sumpop = %.2e %.2e %.2e %.2e %.2e %.2e \n", 02952 popH2g,popH2s,sumpopcollH_deexcit ,sumpopcollH_excit , 02953 sumpopcollH_deexcit/SDIV(popH2s) ,sumpopcollH_excit/SDIV(popH2g));*/ 02954 /*fprintf(ioQQQ,"sumpop = %le sumpopA = %le Av= %le\n", 02955 sumpop1,sumpopA1 , hmi.Average_A );*/ 02956 02957 if( mole.nH2_TRACE >= mole.nH2_trace_full|| (trace.lgTrace && trace.lgTr_H2_Mole) ) 02958 { 02959 fprintf(ioQQQ," H2_LevelPops exit2 Sol dissoc %.2e (TH85 %.2e)", 02960 hmi.H2_Solomon_dissoc_rate_BigH2_H2g + 02961 hmi.H2_Solomon_dissoc_rate_BigH2_H2s , 02962 hmi.H2_Solomon_dissoc_rate_TH85_H2g); 02963 02964 /* Solomon process rate from X into the X continuum with units s-1 02965 * rates are total rate, and rates from H2g and H2s */ 02966 fprintf(ioQQQ," H2g Sol %.2e H2s Sol %.2e", 02967 hmi.H2_Solomon_dissoc_rate_used_H2g , 02968 hmi.H2_Solomon_dissoc_rate_BigH2_H2s ); 02969 02970 /* photoexcitation from H2g to H2s */ 02971 fprintf(ioQQQ," H2g->H2s %.2e (TH85 %.2e)", 02972 hmi.H2_H2g_to_H2s_rate_BigH2 , 02973 hmi.H2_H2g_to_H2s_rate_TH85); 02974 02975 /* add up H2s + hnu => 2H, continuum photodissociation, 02976 * this is not the Solomon process, true continuum, units s-1 */ 02977 fprintf(ioQQQ," H2 con diss %.2e (TH85 %.2e)\n", 02978 hmi.H2_photodissoc_BigH2_H2s , 02979 hmi.H2_photodissoc_TH85); 02980 } 02981 else if( mole.nH2_TRACE ) 02982 { 02983 fprintf(ioQQQ," H2_LevelPops exit1 %8.2f loops:%3li H2/H:%.3e Sol dis old %.3e new %.3e", 02984 fnzone , 02985 loop_h2_pops , 02986 hmi.H2_total / dense.gas_phase[ipHYDROGEN], 02987 old_solomon_rate, 02988 hmi.H2_Solomon_dissoc_rate_BigH2_H2g ); 02989 fprintf(ioQQQ,"\n"); 02990 } 02991 02992 /* >>chng 03 sep 01, add this population - before had just used H2star from chem network */ 02993 /* if big H2 molecule is turned on and used for this zone, use its 02994 * value of H2* (pops of all states with v > 0 ) rather than simple network */ 02995 02996 /* update number of times we have been called */ 02997 ++nCallH2_this_iteration; 02998 02999 /* this will say how many times the large H2 molecule has been called in this zone - 03000 * if not called (due to low H2 abundance) then not need to update its line arrays */ 03001 ++h2.nCallH2_this_zone; 03002 03003 /* >>chng 05 jun 21, 03004 * during search phase we want to use full matrix - save number of levels so that 03005 * we can restore it */ 03006 nXLevelsMatrix = nXLevelsMatrix_save; 03007 03008 /* >>chng 05 jan 19, check how many levels should be in the matrix if first call on 03009 * new zone, and we have a solution */ 03010 /* end loop setting very first LTE H2_populations */ 03011 if( nCallH2_this_iteration && nzone != nzone_nlevel_set ) 03012 { 03013 /* this is fraction of populations to include in matrix */ 03014 # define FRAC 0.99999 03015 /* this loop is over increasing energy */ 03016 double sum_pop = 0.; 03017 nEner = 0; 03018 iElec = 0; 03019 # define PRT false 03020 if( PRT ) fprintf(ioQQQ,"DEBUG pops "); 03021 while( nEner < nLevels_per_elec[0] && sum_pop/hmi.H2_total < FRAC ) 03022 { 03023 03024 /* array of energy sorted indices within X */ 03025 ip = H2_ipX_ener_sort[nEner]; 03026 iVib = ipVib_H2_energy_sort[ip]; 03027 iRot = ipRot_H2_energy_sort[ip]; 03028 sum_pop += H2_old_populations[iElec][iVib][iRot]; 03029 if( PRT ) fprintf(ioQQQ,"\t%.3e ", H2_old_populations[iElec][iVib][iRot]); 03030 ++nEner; 03031 } 03032 if( PRT ) fprintf(ioQQQ,"\n "); 03033 nzone_nlevel_set = nzone; 03034 /*fprintf(ioQQQ,"DEBUG zone %.2f old nmatrix %li proposed nmatrix %li sum_pop %.4e H2_total %.4e\n", 03035 fnzone , nXLevelsMatrix ,nEner , sum_pop,hmi.H2_total); 03036 nXLevelsMatrix = nEner;*/ 03037 # undef FRAC 03038 # undef PRT 03039 } 03040 return; 03041 } 03042 /*lint -e802 possible bad pointer */ 03043 03044 /*H2_cooling evaluate cooling and heating due to H2 molecule, called by 03045 * H2_LevelPops in convergence loop when h2 heating is important, also 03046 * called by CoolEvaluate to get final heating - argument is name of 03047 * routine that called it */ 03048 #if defined(__ICC) && defined(__i386) 03049 #pragma optimization_level 1 03050 #endif 03051 void H2_Cooling( 03052 /* string saying who called this routine, 03053 * "H2lup" call within H2 level populations solver 03054 * "CoolEvaluate" call from main cooling routine */ 03055 const char *chRoutine) 03056 { 03057 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 03058 double heatone , 03059 rate_dn_heat, 03060 rate_up_cool; 03061 long int nColl, 03062 ipHi, ipLo; 03063 double Big1_heat , Big1_cool, 03064 H2_X_col_cool , H2_X_col_heat; 03065 long int ipVib_big_heat_hi,ipVib_big_heat_lo ,ipRot_big_heat_hi , 03066 ipRot_big_heat_lo ,ipVib_big_cool_hi,ipVib_big_cool_lo , 03067 ipRot_big_cool_hi , ipRot_big_cool_lo; 03068 /*# define DEBUG_DIS_DEAT*/ 03069 # ifdef DEBUG_DIS_DEAT 03070 double heatbig; 03071 long int iElecBig , iVibBig , iRotBig; 03072 # endif 03073 03074 /* option to keep track of strongest single heating agent due to collisions 03075 * within X */ 03076 enum {DEBUG_H2_COLL_X_HEAT=false }; 03077 03078 DEBUG_ENTRY( "H2_Cooling()" ); 03079 03080 /* possible debug counters, counter itself, counter to turn on 03081 * debug output, and counter for stopping */ 03082 static long int nCount=-1; 03083 long int nCountDebug = 930, 03084 nCountStop = 940; 03085 ++nCount; 03086 03087 /* nCallH2_this_iteration is not incremented until after the level 03088 * populations have converged the first time. so for the first n calls 03089 * this will return zero, a good idea since populations will be wildly 03090 * incorrect during search for first valid pops */ 03091 if( !h2.lgH2ON || !nCallH2_this_iteration ) 03092 { 03093 hmi.HeatH2Dexc_BigH2 = 0.; 03094 hmi.HeatH2Dish_BigH2 = 0.; 03095 hmi.deriv_HeatH2Dexc_BigH2 = 0.; 03096 return; 03097 } 03098 03099 hmi.HeatH2Dish_BigH2 = 0.; 03100 heatone = 0.; 03101 # ifdef DEBUG_DIS_DEAT 03102 heatbig = 0.; 03103 iElecBig = -1; 03104 iVibBig = -1; 03105 iRotBig = -1; 03106 # endif 03107 /* heating due to dissociation of electronic excited states */ 03108 for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 03109 { 03110 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 03111 { 03112 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 03113 { 03114 /* population, cm-3, in excited state */ 03115 heatone = H2_populations[iElecHi][iVibHi][iRotHi] * 03116 H2_dissprob[iElecHi][iVibHi][iRotHi] * 03117 H2_disske[iElecHi][iVibHi][iRotHi]; 03118 hmi.HeatH2Dish_BigH2 += heatone; 03119 # ifdef DEBUG_DIS_DEAT 03120 if( heatone > heatbig ) 03121 { 03122 heatbig = heatone; 03123 iElecBig = iElecHi; 03124 iVibBig = iVibHi; 03125 iRotBig = iRotHi; 03126 } 03127 # endif 03128 } 03129 } 03130 } 03131 # ifdef DEBUG_DIS_DEAT 03132 fprintf(ioQQQ,"DEBUG H2 dis heat\t%.2f\t%.3f\t%li\t\t%li\t%li\n", 03133 fnzone , 03134 heatbig / SDIV( hmi.HeatH2Dish_BigH2 ) , 03135 iElecBig , 03136 iVibBig , 03137 iRotBig ); 03138 # endif 03139 /* dissociation heating HeatH2Dish_BigH2 was in eV - 03140 * convert to ergs */ 03141 hmi.HeatH2Dish_BigH2 *= EN1EV; 03142 03143 /*fprintf(ioQQQ,"DEBUG H2 heat/dissoc %.3e\n", 03144 hmi.HeatH2Dish_BigH2/ SDIV(hmi.H2_rate_destroy*hmi.H2_total) );*/ 03145 /* now work on collisional heating due to bound-bound 03146 * collisional transitions within X */ 03147 hmi.HeatH2Dexc_BigH2 = 0.; 03148 H2_X_col_cool = 0.; 03149 H2_X_col_heat = 0.; 03150 /* these are the colliders that will be considered as depopulating agents */ 03151 /* the colliders are H, He, H2 ortho, H2 para, H+ */ 03152 /* atomic hydrogen */ 03153 long int nBug1 = 200 , nBug2 = 201; 03154 # if 0 03155 /* fudge(-1) returns number of fudge factors entered on fudge command */ 03156 if( fudge(-1) > 0 ) 03157 { 03158 nBug1 = (long)fudge(0); 03159 nBug2 = (long)fudge(1); 03160 } 03161 # endif 03162 if( DEBUG_H2_COLL_X_HEAT && (nCount == nBug1 || nCount==nBug2) ) 03163 { 03164 FILE *ioBAD=NULL; 03165 if( nCount==nBug1 ) 03166 { 03167 ioBAD = fopen("firstpop.txt" , "w" ); 03168 } 03169 else if( nCount==nBug2 ) 03170 { 03171 ioBAD = fopen("secondpop.txt" , "w" ); 03172 } 03173 for( ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi ) 03174 { 03175 long int ip = H2_ipX_ener_sort[ipHi]; 03176 iVibHi = ipVib_H2_energy_sort[ip]; 03177 iRotHi = ipRot_H2_energy_sort[ip]; 03178 fprintf(ioBAD , "%li\t%li\t%.2e\n", 03179 iVibHi , iRotHi , 03180 H2_populations[0][iVibHi][iRotHi] ); 03181 } 03182 fclose(ioBAD); 03183 } 03184 03185 /* now make sum of all collisions within X itself */ 03186 iElecHi = 0; 03187 iElecLo = 0; 03188 Big1_heat = 0.; 03189 Big1_cool = 0.; 03190 ipVib_big_heat_hi = -1; 03191 ipVib_big_heat_lo = -1; 03192 ipRot_big_heat_hi = -1; 03193 ipRot_big_heat_lo = -1; 03194 ipVib_big_cool_hi = -1; 03195 ipVib_big_cool_lo = -1; 03196 ipRot_big_cool_hi = -1; 03197 ipRot_big_cool_lo = -1; 03198 /* this will be derivative */ 03199 hmi.deriv_HeatH2Dexc_BigH2 = 0.; 03200 for( ipHi=1; ipHi<nLevels_per_elec[iElecHi]; ++ipHi ) 03201 { 03202 long int ip = H2_ipX_ener_sort[ipHi]; 03203 iVibHi = ipVib_H2_energy_sort[ip]; 03204 iRotHi = ipRot_H2_energy_sort[ip]; 03205 if( iVibHi > VIB_COLLID ) 03206 continue; 03207 03208 realnum H2statHi = H2_stat[iElecHi][iVibHi][iRotHi]; 03209 double H2boltzHi = H2_Boltzmann[iElecHi][iVibHi][iRotHi]; 03210 double H2popHi = H2_populations[iElecHi][iVibHi][iRotHi]; 03211 double ewnHi = energy_wn[iElecHi][iVibHi][iRotHi]; 03212 03213 for( ipLo=0; ipLo<ipHi; ++ipLo ) 03214 { 03215 double coolone , oneline; 03216 ip = H2_ipX_ener_sort[ipLo]; 03217 iVibLo = ipVib_H2_energy_sort[ip]; 03218 iRotLo = ipRot_H2_energy_sort[ip]; 03219 if( iVibLo > VIB_COLLID) 03220 continue; 03221 03222 rate_dn_heat = 0.; 03223 03224 /* this sum is total downward heating summed over all colliders */ 03225 mr5ci H2cr = H2_CollRate.begin(iVibHi,iRotHi,iVibLo,iRotLo); 03226 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 03227 /* downward collision rate */ 03228 rate_dn_heat += H2cr[nColl]*collider_density[nColl]; 03229 03230 /* now get upward collisional cooling by detailed balance */ 03231 rate_up_cool = rate_dn_heat * H2_populations[iElecLo][iVibLo][iRotLo] * 03232 /* rest converts into upward collision rate */ 03233 H2statHi / H2_stat[iElecLo][iVibLo][iRotLo] * 03234 H2boltzHi / SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] ); 03235 03236 rate_dn_heat *= H2popHi; 03237 03238 /* net heating due to collisions within X - 03239 * positive if heating, negative is cooling 03240 * this will usually be heating if X is photo pumped 03241 * in printout and in punch heating this is called "H2cX" */ 03242 double conversion = (ewnHi - energy_wn[iElecLo][iVibLo][iRotLo]) * ERG1CM; 03243 heatone = rate_dn_heat * conversion; 03244 coolone = rate_up_cool * conversion; 03245 /* this is net heating, negative if cooling */ 03246 oneline = heatone - coolone; 03247 hmi.HeatH2Dexc_BigH2 += oneline; 03248 03249 /* keep track of heating and cooling separately */ 03250 H2_X_col_cool += coolone; 03251 H2_X_col_heat += heatone; 03252 03253 if( 0 && DEBUG_H2_COLL_X_HEAT && (nCount == 692 || nCount==693) ) 03254 { 03255 static FILE *ioBAD=NULL; 03256 if(ipHi == 1 && ipLo == 0) 03257 { 03258 if( nCount==692 ) 03259 { 03260 ioBAD = fopen("firstheat.txt" , "w" ); 03261 } 03262 else if( nCount==693 ) 03263 { 03264 ioBAD = fopen("secondheat.txt" , "w" ); 03265 } 03266 fprintf(ioBAD,"DEBUG start \n"); 03267 } 03268 fprintf(ioBAD,"DEBUG BAD DAY %li %li %li %li %.3e %.3e\n", 03269 iVibHi , iRotHi, iVibLo , iRotLo , heatone , coolone ); 03270 if( ipHi==nLevels_per_elec[iElecHi]-1 && 03271 ipLo==ipHi-1 ) 03272 fclose(ioBAD); 03273 } 03274 03275 /* derivative wrt temperature - assume exp wrt ground - 03276 * this needs to be divided by square of temperature in wn - 03277 * done at end of loop */ 03278 hmi.deriv_HeatH2Dexc_BigH2 += (realnum)(oneline * ewnHi); 03279 03280 /* oneline is net heating, positive for heating, negative 03281 * for cooling */ 03282 if( DEBUG_H2_COLL_X_HEAT ) 03283 { 03284 if( oneline < Big1_cool ) 03285 { 03286 Big1_cool = oneline; 03287 ipVib_big_cool_hi = iVibHi; 03288 ipVib_big_cool_lo = iVibLo; 03289 ipRot_big_cool_hi = iRotHi; 03290 ipRot_big_cool_lo = iRotLo; 03291 } 03292 else if( oneline > Big1_heat ) 03293 { 03294 Big1_heat = oneline; 03295 ipVib_big_heat_hi = iVibHi; 03296 ipVib_big_heat_lo = iVibLo; 03297 ipRot_big_heat_hi = iRotHi; 03298 ipRot_big_heat_lo = iRotLo; 03299 } 03300 } 03301 03302 /* this would be a major logical error */ 03303 ASSERT( 03304 (rate_up_cool==0 && rate_dn_heat==0) || 03305 (energy_wn[iElecHi][iVibHi][iRotHi] > energy_wn[iElecLo][iVibLo][iRotLo]) ); 03306 }/* end loop over lower levels, all collisions within X */ 03307 }/* end loop over upper levels, all collisions within X */ 03308 03309 /* this debug statement will identify the single strongest heating or 03310 * cooling agent within X and give the lowest 7 level populations */ 03311 if( DEBUG_H2_COLL_X_HEAT ) 03312 { 03313 /* nCount counts number of calls through here, 03314 * nCountDebug is count to turn on debug prints, 03315 * nCountStop is count to abort code */ 03316 if(nCount>nCountDebug ) 03317 { 03318 fprintf(ioQQQ, 03319 "DEBUG H2_Cooling A %li %15s, Te %.3e net Heat(Xcol) %.2e " 03320 "heat %.2e cool %.2e H+/0 " 03321 "%.2e n(H2)%.3e Sol rat %.3e grn J1->0%.2e frac heat 1 line " 03322 "%.2e Hi(v,j)%li %li Lo(v,J)%li %li frac cool 1 line %.2e " 03323 "Hi(v,j)%li %li Lo(v,J)%li %li POP(J=1,13)", 03324 nCount , chRoutine, 03325 phycon.te, 03326 hmi.HeatH2Dexc_BigH2 , 03327 H2_X_col_cool , 03328 H2_X_col_heat , 03329 dense.xIonDense[ipHYDROGEN][1]/SDIV(dense.xIonDense[ipHYDROGEN][0]), 03330 hmi.H2_total, 03331 hmi.H2_Solomon_dissoc_rate_BigH2_H2g, 03332 hmi.rate_grain_h2_J1_to_J0 , 03333 Big1_heat/hmi.HeatH2Dexc_BigH2 , 03334 ipVib_big_heat_hi , ipRot_big_heat_hi , ipVib_big_heat_lo , ipRot_big_heat_lo , 03335 Big1_cool/hmi.HeatH2Dexc_BigH2 , 03336 ipVib_big_cool_hi , ipRot_big_cool_hi , ipVib_big_cool_lo , ipRot_big_cool_lo ); 03337 03338 for( iRotLo=0; iRotLo<14; ++iRotLo ) 03339 { 03340 fprintf(ioQQQ,"\t%.2e" , 03341 H2_populations[0][0][iRotLo]/hmi.H2_total ); 03342 } 03343 fprintf(ioQQQ,"\t%li\n",nCount); 03344 /* now give collision rates for strongest heat/cool level */ 03345 fprintf(ioQQQ,"DEBUG H2_Cooling B heat Coll Rate (lrg col) dn,up" ); 03346 double HeatNet = 0.; 03347 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 03348 { 03349 fprintf(ioQQQ,"\t%.2e" , 03350 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]* 03351 collider_density[nColl] ); 03352 fprintf(ioQQQ,"\t%.2e" , 03353 /* downward collision rate */ 03354 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]* 03355 collider_density[nColl]* 03356 /* rest converts into upward collision rate */ 03357 H2_stat[0][ipVib_big_heat_hi][ipRot_big_heat_hi] / H2_stat[0][ipVib_big_heat_lo][ipRot_big_heat_lo] * 03358 H2_Boltzmann[0][ipVib_big_heat_hi][ipRot_big_heat_hi] / 03359 SDIV( H2_Boltzmann[0][ipVib_big_heat_lo][ipRot_big_heat_lo] ) ); 03360 HeatNet += 03361 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]* 03362 collider_density[nColl]*H2_populations[iElecLo][ipRot_big_heat_hi][ipVib_big_heat_lo] 03363 - 03364 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]* 03365 collider_density[nColl]* 03366 /* rest converts into upward collision rate */ 03367 H2_stat[0][ipVib_big_heat_hi][ipRot_big_heat_hi] / H2_stat[0][ipVib_big_heat_lo][ipRot_big_heat_lo] * 03368 H2_Boltzmann[0][ipVib_big_heat_hi][ipRot_big_heat_hi] / 03369 SDIV( H2_Boltzmann[0][ipVib_big_heat_lo][ipRot_big_heat_lo] ) * 03370 H2_populations[iElecLo][ipVib_big_heat_lo][ipRot_big_heat_lo] ; 03371 } 03372 /* HeatNet includes the level populations, rates do not */ 03373 fprintf( ioQQQ , " HeatNet %.2e",HeatNet); 03374 fprintf(ioQQQ,"\n"); 03375 fprintf(ioQQQ,"DEBUG H2_Cooling C cool Coll Rate (lrg col) dn,up" ); 03376 double CoolNet = 0.; 03377 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 03378 { 03379 fprintf(ioQQQ,"\t%.2e" , 03380 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]* 03381 collider_density[nColl] ); 03382 fprintf(ioQQQ,"\t%.2e" , 03383 /* downward collision rate */ 03384 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]* 03385 collider_density[nColl]* 03386 /* rest converts into upward collision rate */ 03387 H2_stat[0][ipVib_big_cool_hi][ipRot_big_cool_hi] / H2_stat[0][ipVib_big_cool_lo][ipRot_big_cool_lo] * 03388 H2_Boltzmann[0][ipVib_big_cool_hi][ipRot_big_cool_hi] / 03389 SDIV( H2_Boltzmann[0][ipVib_big_cool_lo][ipRot_big_cool_lo] ) ); 03390 CoolNet += 03391 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]* 03392 collider_density[nColl]*H2_populations[iElecLo][ipRot_big_cool_hi][ipVib_big_cool_lo] 03393 - 03394 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]* 03395 collider_density[nColl]* 03396 /* rest converts into upward collision rate */ 03397 H2_stat[0][ipVib_big_cool_hi][ipRot_big_cool_hi] / H2_stat[0][ipVib_big_cool_lo][ipRot_big_cool_lo] * 03398 H2_Boltzmann[0][ipVib_big_cool_hi][ipRot_big_cool_hi] / 03399 SDIV( H2_Boltzmann[0][ipVib_big_cool_lo][ipRot_big_cool_lo] ) * 03400 H2_populations[iElecLo][ipVib_big_cool_lo][ipRot_big_cool_lo] ; 03401 } 03402 /* CoolNet includes the level populations, rates do not */ 03403 fprintf( ioQQQ , " CoolNet %.2e",CoolNet); 03404 fprintf(ioQQQ,"\n"); 03405 } 03406 } 03407 03408 /* this is inside h2 cooling, and is called extra times when H2 heating is important */ 03409 if( PRT_POPS ) 03410 fprintf(ioQQQ, 03411 " DEBUG H2 heat fnzone\t%.2f\trenrom\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n", 03412 fnzone , 03413 H2_renorm_chemistry , 03414 phycon.te , 03415 hmi.HeatH2Dexc_BigH2, 03416 hmi.HeatH2Dexc_BigH2/thermal.ctot); 03417 03418 /* this is derivative of collisional heating wrt temperature - needs 03419 * to be divided by square of temperature in wn */ 03420 hmi.deriv_HeatH2Dexc_BigH2 /= (realnum)POW2(phycon.te_wn); 03421 03422 { 03423 enum {DEBUG_LOC=false }; 03424 if( DEBUG_H2_COLL_X_HEAT && DEBUG_LOC && 03425 (fabs(hmi.HeatH2Dexc_BigH2) > SMALLFLOAT) ) 03426 { 03427 int iVib = 0; 03428 03429 /*fprintf(ioQQQ," H2_cooling pops\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 03430 H2_populations[0][iVib][0]/hmi.H2_total, 03431 H2_populations[0][iVib][1]/hmi.H2_total, 03432 H2_populations[0][iVib][2]/hmi.H2_total, 03433 H2_populations[0][iVib][3]/hmi.H2_total, 03434 H2_populations[0][iVib][4]/hmi.H2_total, 03435 H2_populations[0][iVib][5]/hmi.H2_total);*/ 03436 03437 iElecHi = iElecLo = 0; 03438 iVibHi = iVibLo = 0; 03439 iRotHi = 7; 03440 iRotLo = 5; 03441 rate_dn_heat = rate_up_cool = 0.; 03442 /* this sum is total downward heating */ 03443 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 03444 { 03445 rate_dn_heat += 03446 H2_populations[iElecHi][iVibHi][iRotHi] * 03447 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl]* 03448 collider_density[nColl]; 03449 03450 /* now get upward collisional cooling by detailed balance */ 03451 rate_up_cool += 03452 H2_populations[iElecLo][iVibLo][iRotLo] * 03453 /* downward collision rate */ 03454 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl]* 03455 collider_density[nColl]* 03456 /* rest converts into upward collision rate */ 03457 H2_stat[iElecHi][iVibHi][iRotHi] / H2_stat[iElecLo][iVibLo][iRotLo] * 03458 H2_Boltzmann[iElecHi][iVibHi][iRotHi] / 03459 SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] ); 03460 } 03461 03462 fprintf(ioQQQ,"DEBUG H2_cooling D pop %li ov %li\t%.3e\tdn up 31\t%.3e\t%.3e\n", 03463 iRotHi , iRotLo , 03464 H2_populations[0][iVib][iRotHi]/H2_populations[0][iVib][iRotLo], 03465 rate_dn_heat, 03466 rate_up_cool); 03467 if( nCount>= nCountStop ) 03468 cdEXIT(EXIT_FAILURE); 03469 } 03470 } 03471 { 03472 enum {DEBUG_LOC=false}; 03473 if( DEBUG_LOC ) 03474 { 03475 static long nzdone=-1 , nzincre; 03476 if( nzone!=nzdone ) 03477 { 03478 nzdone = nzone; 03479 nzincre = -1; 03480 } 03481 ++nzincre; 03482 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", 03483 fnzone, 03484 nzincre, 03485 phycon.te, 03486 hmi.H2_total , 03487 hmi.HeatH2Dexc_BigH2, 03488 hmi.deriv_HeatH2Dexc_BigH2 , 03489 hmi.HeatH2Dish_BigH2); 03490 03491 } 03492 } 03493 03494 # if 0 03495 /* this can be noisy due to finite accuracy of solution, so take average with 03496 * previous value */ 03497 /*>>chng 04 mar 01, do not take average */ 03498 if( 1 || nzone <1 || old_HeatH2Dexc==0. || nCallH2_this_iteration <2) 03499 { 03500 old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2; 03501 } 03502 else 03503 { 03504 hmi.HeatH2Dexc_BigH2 = (hmi.HeatH2Dexc_BigH2+old_HeatH2Dexc)/2.f; 03505 old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2; 03506 } 03507 # endif 03508 { 03509 enum {DEBUG_LOC=false}; 03510 if( DEBUG_LOC /*&& DEBUG_H2_COLL_X_HEAT*/ ) 03511 { 03512 fprintf(ioQQQ,"DEBUG H2_cooling E %15s %c vib deex %li Te %.3e net heat %.3e cool %.3e heat %.3e\n", 03513 chRoutine , 03514 TorF(conv.lgSearch), 03515 nCount, 03516 phycon.te, hmi.HeatH2Dexc_BigH2, 03517 H2_X_col_cool , 03518 H2_X_col_heat /*, 03519 H2_populations[0][0][7]/SDIV(H2_populations[0][0][5]) , 03520 H2_populations[0][0][13]/SDIV(H2_populations[0][0][11])*/ ); 03521 if( 0 && nCount > nCountStop ) 03522 { 03523 cdEXIT( EXIT_FAILURE ); 03524 } 03525 } 03526 } 03527 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 03528 fprintf(ioQQQ, 03529 " H2_Cooling Ctot\t%.4e\t HeatH2Dish_BigH2 \t%.4e\t HeatH2Dexc_BigH2 \t%.4e\n" , 03530 thermal.ctot , 03531 hmi.HeatH2Dish_BigH2 , 03532 hmi.HeatH2Dexc_BigH2 ); 03533 03534 /* when we are very far from solution, during search phase, collisions within 03535 * X can be overwhelmingly large heating and cooling terms, which nearly 03536 * cancel out. Some dense cosmic ray heated clouds could not find correct 03537 * initial solution due to noise introduced by large net heating which was 03538 * the very noisy tiny difference between very large heating and cooling 03539 * terms. Do not include collisions with x as heat/cool during the 03540 * initial search phase */ 03541 if( conv.lgSearch ) 03542 hmi.HeatH2Dexc_BigH2 = 0.; 03543 return; 03544 } 03545 03546 03547 /*cdH2_colden return column density in H2, negative -1 if cannot find state, 03548 * header is cdDrive */ 03549 double cdH2_colden( long iVib , long iRot ) 03550 { 03551 03552 /*if iVib is negative, return 03553 * total column density - iRot=0 03554 * ortho column density - iRot 1 03555 * para column density - iRot 2 03556 * else return column density in iVib, iRot */ 03557 if( iVib < 0 ) 03558 { 03559 if( iRot==0 ) 03560 { 03561 /* return total H2 column density */ 03562 return( h2.ortho_colden + h2.para_colden ); 03563 } 03564 else if( iRot==1 ) 03565 { 03566 /* return ortho H2 column density */ 03567 return h2.ortho_colden; 03568 } 03569 else if( iRot==2 ) 03570 { 03571 /* return para H2 column density */ 03572 return h2.para_colden; 03573 } 03574 else 03575 { 03576 fprintf(ioQQQ," iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n"); 03577 return -1.; 03578 } 03579 } 03580 else if( h2.lgH2ON ) 03581 { 03582 /* this branch want state specific column density, which can only result from 03583 * evaluation of big molecule */ 03584 int iElec = 0; 03585 if( iRot <0 || iVib >h2.nVib_hi[iElec] || iRot > h2.nRot_hi[iElec][iVib]) 03586 { 03587 fprintf(ioQQQ," iVib and iRot must lie within X, returning -2.\n"); 03588 fprintf(ioQQQ," iVib must be <= %li and iRot must be <= %li.\n", 03589 h2.nVib_hi[iElec],h2.nRot_hi[iElec][iVib]); 03590 return -2.; 03591 } 03592 else 03593 { 03594 return H2_X_colden[iVib][iRot]; 03595 } 03596 } 03597 /* error condition - no valid parameter */ 03598 else 03599 return -1; 03600 } 03601 03602 /*H2_Colden maintain H2 column densities within X */ 03603 void H2_Colden( const char *chLabel ) 03604 { 03605 long int iVib , iRot; 03606 03607 /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */ 03608 if( !h2.lgH2ON /*|| !h2.nCallH2_this_zone*/ ) 03609 return; 03610 03611 DEBUG_ENTRY( "H2_Colden()" ); 03612 03613 if( strcmp(chLabel,"ZERO") == 0 ) 03614 { 03615 /* zero out formation rates and column densites */ 03616 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 03617 { 03618 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 03619 { 03620 /* space for the rotation quantum number */ 03621 H2_X_colden[iVib][iRot] = 0.; 03622 H2_X_colden_LTE[iVib][iRot] = 0.; 03623 } 03624 } 03625 } 03626 03627 else if( strcmp(chLabel,"ADD ") == 0 ) 03628 { 03629 /* add together column densities */ 03630 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 03631 { 03632 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 03633 { 03634 /* state specific H2 column density */ 03635 H2_X_colden[iVib][iRot] += (realnum)(H2_populations[0][iVib][iRot]*radius.drad_x_fillfac); 03636 /* LTE state specific H2 column density - H2_populations_LTE is normed to unity 03637 * so must be multiplied by total H2 density */ 03638 H2_X_colden_LTE[iVib][iRot] += (realnum)(H2_populations_LTE[0][iVib][iRot]* 03639 hmi.H2_total*radius.drad_x_fillfac); 03640 } 03641 } 03642 } 03643 03644 /* we will not print column densities so skip that - if not print then we have a problem */ 03645 else if( strcmp(chLabel,"PRIN") != 0 ) 03646 { 03647 fprintf( ioQQQ, " H2_Colden does not understand the label %s\n", 03648 chLabel ); 03649 cdEXIT(EXIT_FAILURE); 03650 } 03651 03652 return; 03653 } 03654 03655 /*H2_DR choose next zone thickness based on H2 big molecule */ 03656 double H2_DR(void) 03657 { 03658 return BIGFLOAT; 03659 } 03660 03661 /*H2_RT_OTS - add H2 ots fields */ 03662 void H2_RT_OTS( void ) 03663 { 03664 03665 long int iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo; 03666 03667 /* do not compute if H2 not turned on, or not computed for these conditions */ 03668 if( !h2.lgH2ON || !h2.nCallH2_this_zone ) 03669 return; 03670 03671 DEBUG_ENTRY( "H2_RT_OTS()" ); 03672 03673 /* loop over all possible lines and set H2_populations, and quantities that depend on escape prob, dest, etc */ 03674 long int lim_elec_hi = 0; 03675 for( iElecHi=0; iElecHi<=lim_elec_hi; ++iElecHi ) 03676 { 03677 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 03678 { 03679 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 03680 { 03681 double H2popHi = H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop; 03682 long int lim_elec_lo = 0; 03683 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 03684 { 03685 /* want to include all vibration states in lower level if different electronic level, 03686 * but only lower vibration levels if same electronic level */ 03687 long int nv = h2.nVib_hi[iElecLo]; 03688 if( iElecLo==iElecHi ) 03689 nv = iVibHi; 03690 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 03691 { 03692 long nr = h2.nRot_hi[iElecLo][iVibLo]; 03693 if( iElecLo==iElecHi && iVibHi==iVibLo ) 03694 nr = iRotHi-1; 03695 03696 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 03697 { 03698 /* >>chng 05 feb 07, use lgH2_line_exists */ 03699 if( iElecHi==0 && lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 03700 { 03701 /* ots destruction rate */ 03702 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots = 03703 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul * H2popHi * 03704 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Pdest; 03705 03706 /* dump the ots rate into the stack - but only for ground electronic state*/ 03707 RT_OTS_AddLine( 03708 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots, 03709 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont ); 03710 } 03711 } 03712 } 03713 } 03714 } 03715 } 03716 } 03717 03718 return; 03719 } 03720 /*lint +e802 possible bad pointer */ 03721