47 static long int nhit = 0,
50 static double TeEvalCS = 0., TeEvalCS_21cm=0.;
51 static double TeUsedBrems=-1.f;
52 static int nzoneUsedBrems=-1;
54 static double electron_rate_21cm,
67 static double oltcool=0.,
70 long int coolnum, coolcal;
87 fprintf(
ioQQQ,
" COOLR TE:%.4e zone %li %li Cool:%.4e Heat:%.4e eden:%.4e edenTrue:%.4e\n",
100 if(
gv.lgGrainPhysicsOn )
112 if(
gv.lgDustOn() &&
gv.lgDHetOn )
115 thermal.heating[0][13] =
gv.GasHeatPhotoEl;
123 thermal.heating[0][25] =
gv.GasHeatTherm;
132 else if(
gv.lgBakesPAH_heat )
136 thermal.heating[0][13] =
gv.GasHeatPhotoEl;
160 hmi.HeatH2Dish_used = 0.;
161 hmi.HeatH2Dexc_used = 0.;
162 hmi.deriv_HeatH2Dexc_used = 0.;
171 if(
h2.lgEnabled &&
hmi.lgH2_Thermal_BigH2 )
179 hmi.HeatH2Dish_used =
h2.HeatDiss;
180 hmi.HeatH2Dexc_used =
h2.HeatDexc;
182 fprintf(
ioQQQ,
"DEBUG big %.2f\t%.5e\t%.2e\t%.2e\t%.2e\n",
187 hmi.deriv_HeatH2Dexc_used = -
h2.HeatDexc_deriv;
191 hmi.HeatH2Dish_used = 0;
192 hmi.HeatH2Dexc_used = 0;
193 hmi.deriv_HeatH2Dexc_used = 0;
197 else if(
hmi.chH2_small_model_type ==
'T' )
201 hmi.HeatH2Dish_used =
hmi.HeatH2Dish_TH85;
202 hmi.HeatH2Dexc_used =
hmi.HeatH2Dexc_TH85;
203 hmi.deriv_HeatH2Dexc_used =
hmi.deriv_HeatH2Dexc_TH85;
205 else if(
hmi.chH2_small_model_type ==
'H' )
208 hmi.HeatH2Dish_used =
hmi.HeatH2Dish_BHT90;
209 hmi.HeatH2Dexc_used =
hmi.HeatH2Dexc_BHT90;
210 hmi.deriv_HeatH2Dexc_used =
hmi.deriv_HeatH2Dexc_BHT90;
212 else if(
hmi.chH2_small_model_type ==
'B')
215 hmi.HeatH2Dish_used =
hmi.HeatH2Dish_BD96;
216 hmi.HeatH2Dexc_used =
hmi.HeatH2Dexc_BD96;
217 hmi.deriv_HeatH2Dexc_used =
hmi.deriv_HeatH2Dexc_BD96;
219 else if(
hmi.chH2_small_model_type ==
'E')
222 hmi.HeatH2Dish_used =
hmi.HeatH2Dish_ELWERT;
223 hmi.HeatH2Dexc_used =
hmi.HeatH2Dexc_ELWERT;
224 hmi.deriv_HeatH2Dexc_used =
hmi.deriv_HeatH2Dexc_ELWERT;
238 if( (*diatom)->lgEnabled &&
mole_global.lgStancil )
240 heat += (*diatom)->Cont_Diss_Heat_Rate();
260 if(
hmi.HeatH2Dexc_used < 0. )
284 rothi = pow(10.,-19.24 + 0.474*x - 1.247*x*x);
297 rotlow = pow(10.,-22.90 - 0.553*x - 1.148*x*x)*qn;
302 CoolHeavy.h2line +=
hmi.H2_total*rothi/(1. + rothi/rotlow);
309 enum {DEBUG_LOC=
false};
312 fprintf(
ioQQQ,
"h2coolbug\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
317 hmi.HMinus_photo_rate,
333 factor/(1416.+
phycon.sqrte*
hmi.H2_total * (1. + 3.*factor));
339 double chemical_heating =
mole.chem_heat();
340 thermal.heating[0][29] =
MAX2(0.,chemical_heating);
372 enum {DEBUG_LOC=
false};
375 fprintf(
ioQQQ,
"CoolEvaluate debuggg\t%.2f\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n",
401 for(
long int nelem=ipISO; nelem <
LIMELM; nelem++ )
419 nzone != nzoneUsedBrems )
421 double BremsThisEnergy;
427 nzoneUsedBrems =
nzone;
438 double bhfac, bhMinusfac;
440 long int ion_lo , ion_hi;
454 for(
long int i=
rfield.ipEnergyBremsThin; i < limit; i++ )
462 CoolHeavy.brems_cool_h += BremsThisEnergy;
466 CoolHeavy.brems_cool_hminus += BremsThisEnergy *
opac.OpacStack[i-1+
opac.iphmra];
471 CoolHeavy.brems_cool_hminus *= bhMinusfac;
475 for(
long int i=
rfield.ipEnergyBremsThin; i < limit; i++ )
482 CoolHeavy.brems_cool_he += BremsThisEnergy;
491 for(
long int ion=1; ion<=
LIMELM; ++ion )
496 if(
dense.lgElmtOn[nelem] && ion<=nelem+1 )
498 sumion[ion] +=
dense.xIonDense[nelem][ion];
506 for(
long ipMol = 0; ipMol<
mole_global.num_calc; ipMol++ )
527 while( sumion[ion_lo]==0 && ion_lo<
LIMELM-1 )
530 while( sumion[ion_hi]==0 && ion_hi>0 )
536 for(
long int i=
rfield.ipEnergyBremsThin; i < limit; i++ )
538 BremsThisEnergy = 0.;
539 for(
long int ion=ion_lo; ion<=ion_hi; ++ion )
540 BremsThisEnergy += sumion[ion]*
rfield.gff[ion][i];
550 enum {DEBUG_LOC=
false};
551 if( DEBUG_LOC &&
nzone>60 )
553 double sumfield = 0., sumtot=0., sum1=0., sum2=0.;
554 for(
long int i=
rfield.ipEnergyBremsThin; i<limit; i++ )
559 sum2 +=
opac.FreeFreeOpacity[i]*
rfield.flux[0][i];
561 fprintf(
ioQQQ,
"DEBUG brems heat\t%.2f\t%.3e\t%.3e\t%.3e\t%e\t%.3e\t%.3e\n",
564 sumtot/
SDIV(sumfield) ,
568 opac.FreeFreeOpacity[1218]);
612 if(
dense.lgElmtOn[nelem] )
615 long limit_lo =
MAX2( 1 ,
dense.IonLow[nelem] );
617 for(
long int ion=limit_lo; ion<=limit_hi; ++ion )
625 double one =
dense.xIonDense[nelem][ion] *
ionbal.RR_rate_coef_used[nelem][ion-1]*
680 if(
wind.lgBallistic() )
686 else if(
dynamics.lgTimeDependentStatic &&
721 enum {DEBUG_LOC=
false};
726 double teval[
N21CM_TE]={2.,5.,10.,20.,50.,100.,200.,500.,1000.,
727 2000.,3000.,5000.,7000.,10000.,15000.,20000.};
731 ioQQQ,
"DEBUG 21 cm deex Te=\t%.2e\tH0=\t%.2e\tp=\t%.2e\te=\t%.2e\n",
746 TeEvalCS_21cm =
phycon.te;
750 cs = (electron_rate_21cm *
dense.eden +
760 for(
long int i=1; i <
nHFLines; i++ )
779 for(
long int i=1; i <
nHFLines; i++ )
803 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
812 for(
int nelem=0; nelem <
LIMELM; nelem++ )
814 for(
int ion=0; ion<=nelem+1; ++ion )
816 xIonDenseSave[nelem][ion] =
dense.xIonDense[nelem][ion];
818 if(
dense.lgIonChiantiOn[nelem][ion] ||
dense.lgIonStoutOn[nelem][ion] )
819 dense.xIonDense[nelem][ion] = 0.;
827 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
835 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
843 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
852 fprintf(
ioQQQ,
"DEBUG dCdT Ne %.3e dHdT %.3e\n",
thermal.dCooldT
854 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
861 fprintf(
ioQQQ,
"DEBUG dCdT 8 %.3e dHdT %.3e\n",
thermal.dCooldT
863 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
869 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
877 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
885 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
893 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
899 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
905 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
911 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
919 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
925 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
931 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
937 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
944 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
952 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
958 for( coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
967 for(
int nelem=0; nelem <
LIMELM; nelem++ )
969 for(
int ion=0; ion<=nelem+1; ++ion )
971 dense.xIonDense[nelem][ion] = xIonDenseSave[nelem][ion];
979 for(
int coolcal = coolnum; coolcal <
thermal.ncltot; coolcal++ )
987 enum {DEBUG_LOC=
false};
995 for(
long ipSpecies=1; ipSpecies<
nSpecies; ipSpecies++ )
1000 printf(
"DEBUG Max\t%li" ,
dBaseSpecies[0].numLevels_max );
1001 for(
long ipSpecies=1; ipSpecies<
nSpecies; ipSpecies++ )
1003 printf(
"\t%li" ,
dBaseSpecies[ipSpecies].numLevels_max );
1007 printf(
"DEBUG Local\t%li" ,
dBaseSpecies[0].numLevels_local );
1008 for(
long ipSpecies=1; ipSpecies<
nSpecies; ipSpecies++ )
1010 printf(
"\t%li" ,
dBaseSpecies[ipSpecies].numLevels_local );
1024 fprintf(
ioQQQ,
" COOLR; cooling is <=0, this is impossible.\n" );
1032 fprintf(
ioQQQ,
" COOLR; cooling slope <=0, this is impossible.\n" );
1035 fprintf(
ioQQQ,
" Probably due to very low density.\n" );
1047 if( (((((!
thermal.lgTemperatureConstant) && *tot < 0.) &&
called.lgTalk) &&
1051 " NOTE Negative cooling, zone %4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n",
1068 deriv = (oltcool - *tot)/(oldtemp -
phycon.te);
1075 if(
nzone != nzSave )