Changeset 1175

Show
Ignore:
Timestamp:
06/08/07 10:50:41 (1 year ago)
Author:
gary
Message:

readme_data.htm - correct name of FeII energy levels file
cool_eval.cpp - minor changes in comments, no change in executable
init_defaults.cpp rename to init_defaults_preparse.cpp - this is the initialization routine that sets variables that are possibly reset by the parser

mole_h2.cpp - many changes in debug print statements, no change in calculations

mole_h2_etc.cpp - add SDIV in denominator one time, clean up comments

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved
  • trunk/data/readme_data.htm

    r919 r1175  
    161161<h3>fe2nn.dat</h3> 
    162162<p>nn(142) - actual level numbers for each of Zhang and Pradhan levels</p> 
    163 <h3>fe2levels.dat</h3> 
    164 <p>This gives the spectroscopic levels sorted in increasing order within the 
    165 model atom.</p> 
     163<h3> 
     164        fe2energies.dat</h3> 
     165<p>This gives the energy levels in increasing order.</p> 
    166166<hr> 
    167167<h2>Molecular hydrogen files</h2> 
  • trunk/source/cool_eval.cpp

    r1122 r1175  
    139139        if( hmi.lgNoH2Mole ) 
    140140        { 
     141                /* this branch - do not include molecules */ 
    141142                hmi.hmicol = 0.; 
    142143                CoolHeavy.brems_cool_hminus = 0.; 
     
    159160        else 
    160161        { 
    161  
    162                 /* save it */ 
     162                /* save various molecular heating/cooling agent */ 
    163163                thermal.heating[0][15] = hmi.hmihet; 
    164164                thermal.heating[0][16] = hmi.h2plus_heat; 
     
    227227                /* add to heating is net heating is positive */ 
    228228                thermal.heating[0][8] = MAX2(0.,hmi.HeatH2Dexc_used); 
     229 
    229230                /* add to cooling if net heating is negative */ 
    230231                CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used)); 
     
    234235                /*thermal.dCooldT += MAX2(0.,-hmi.HeatH2Dexc_used)* ( 30172. * thermal.tsq1 - thermal.halfte );*/ 
    235236                /* >>chng 04 jan 25, check sign to prevent cooling from entering here, 
    236                  * also enter neg sign since going into cooling stack (bug), in headsum 
     237                 * also enter neg sign since going into cooling stack (bug), in heatsum 
    237238                 * same term adds to deriv of heating */ 
    238239                if( hmi.HeatH2Dexc_used < 0. ) 
  • trunk/source/date.h

    r1168 r1175  
    99#define YEAR    107 
    1010#define MONTH   5 
    11 #define DAY     4 
     11#define DAY     7 
  • trunk/source/mole_h2.cpp

    r1167 r1175  
    33943394        return; 
    33953395} 
     3396/*lint -e668 possible bad pointer */ 
     3397/*lint -e802 possible bad pointer */ 
    33963398 
    33973399/*H2_cooling evaluate cooling and heating due to H2 molecule, called by  
     
    34013403void H2_Cooling( 
    34023404        /* string saying who called this routine,  
    3403          * H2lup call within H2 level populations solver 
    3404          * CoolEvaluate call from main cooling routine */ 
     3405         * "H2lup" call within H2 level populations solver 
     3406         * "CoolEvaluate" call from main cooling routine */ 
    34053407         const char *chRoutine) 
    34063408{ 
    3407         static long int nCount=0; 
    34083409        long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo; 
    34093410        double heatone , 
     
    34123413        long int nColl, 
    34133414                ipHi, ipLo; 
    3414         double Big1_heat , Big1_cool; 
     3415        double Big1_heat , Big1_cool, 
     3416                H2_X_col_cool , H2_X_col_heat; 
    34153417        long int ipVib_big_heat_hi,ipVib_big_heat_lo ,ipRot_big_heat_hi , 
    34163418                ipRot_big_heat_lo ,ipVib_big_cool_hi,ipVib_big_cool_lo , 
    34173419                ipRot_big_cool_hi ,     ipRot_big_cool_lo; 
    3418         static double old_HeatH2Dexc=0.; 
    34193420/*#     define DEBUG_DIS_DEAT*/ 
    34203421#       ifdef DEBUG_DIS_DEAT 
     
    34223423        long int iElecBig , iVibBig , iRotBig; 
    34233424#       endif 
    3424         /* option to keep track of strongest single heating agent */ 
    3425 #       define DEBUG_H2_COLL_X_HEAT     false 
     3425 
     3426        /* option to keep track of strongest single heating agent due to collisions 
     3427         * within X */ 
     3428        enum {DEBUG_H2_COLL_X_HEAT=false }; 
    34263429 
    34273430        DEBUG_ENTRY( "H2_Cooling()" ); 
     3431 
     3432        /* possible debug counters, counter itself, counter to turn on 
     3433         * debug output, and counter for stopping */ 
     3434        static long int nCount=-1; 
     3435        long int nCountDebug = 930, 
     3436                nCountStop = 940; 
     3437        ++nCount; 
    34283438 
    34293439        /* nCallH2_this_iteration is not incremented until after the level 
     
    34393449                return; 
    34403450        } 
    3441  
    3442         old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2; 
    34433451 
    34443452        hmi.HeatH2Dish_BigH2 = 0.; 
     
    34913499         * collisional transitions within X */ 
    34923500        hmi.HeatH2Dexc_BigH2 = 0.; 
     3501        H2_X_col_cool = 0.; 
     3502        H2_X_col_heat = 0.; 
    34933503        /* these are the colliders that will be considered as depopulating agents */ 
    34943504        /* the colliders are H, He, H2 ortho, H2 para, H+ */ 
    34953505        /* atomic hydrogen */ 
    3496         if( DEBUG_H2_COLL_X_HEAT && nCount == 692 || nCount==693 ) 
    3497         { 
    3498                 /*lint -e668 possible bad pointer */ 
    3499                 /*lint -e802 possible bad pointer */ 
     3506        long int nBug1 = 200 , nBug2 = 201;  
     3507        /* fudge(-1) returns number of fudge factors entered on fudge command */ 
     3508        if( fudge(-1) > 0 ) 
     3509        { 
     3510                nBug1 = (long)fudge(0); 
     3511                nBug2 = (long)fudge(1); 
     3512        } 
     3513        if( DEBUG_H2_COLL_X_HEAT && nCount == nBug1 || nCount==nBug2 ) 
     3514        { 
    35003515                FILE *ioBAD=NULL; 
    3501                 if( nCount==692
     3516                if( nCount==nBug1
    35023517                { 
    35033518                        ioBAD = fopen("firstpop.txt" , "w" ); 
    35043519                } 
    3505                 else if( nCount==693
     3520                else if( nCount==nBug2
    35063521                { 
    35073522                        ioBAD = fopen("secondpop.txt" , "w" ); 
     
    35173532                } 
    35183533                fclose(ioBAD); 
    3519                 /*lint +e668 possible bad pointer */ 
    3520                 /*lint +e802 possible bad pointer */ 
    35213534        } 
    35223535 
     
    35463559                for( ipLo=0; ipLo<ipHi; ++ipLo ) 
    35473560                { 
    3548                         double oneline; 
     3561                        double coolone , oneline; 
    35493562                        ip = H2_ipX_ener_sort[ipLo]; 
    35503563                        iVibLo = ipVib_H2_energy_sort[ip]; 
     
    35763589                        }/* end loop over all colliders */ 
    35773590 
    3578                         /* net heating due to collisions within X - this will usually be heating - 
     3591                        /* net heating due to collisions within X -  
     3592                         * positive if heating, negative is cooling 
     3593                         * this will usually be heating if X is photo pumped 
    35793594                         * in printout and in punch heating this is called "H2cX" */ 
    3580                         oneline = (rate_dn_heat - rate_up_cool)* 
     3595                        heatone = rate_dn_heat *  
    35813596                                (energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]) * 
    35823597                                ERG1CM; 
     3598                        coolone = rate_up_cool* 
     3599                                (energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]) * 
     3600                                ERG1CM; 
     3601                        /* this is net heating, negative if cooling */ 
     3602                        oneline = heatone - coolone; 
    35833603                        hmi.HeatH2Dexc_BigH2 += oneline; 
    35843604 
    3585                         if( DEBUG_H2_COLL_X_HEAT && (nCount == 692 || nCount==693) ) 
     3605                        /* keep track of heating and cooling separately */ 
     3606                        H2_X_col_cool += coolone; 
     3607                        H2_X_col_heat += heatone; 
     3608 
     3609                        if( 0 && DEBUG_H2_COLL_X_HEAT && (nCount == 692 || nCount==693) ) 
    35863610                        { 
    35873611                                static FILE *ioBAD=NULL; 
     
    35983622                                        fprintf(ioBAD,"DEBUG start \n"); 
    35993623                                } 
    3600                                 fprintf(ioBAD,"DEBUG BAD DAY %li %li %li %li %.3e\n", 
    3601                                         iVibHi , iRotHi, iVibLo , iRotLo , oneline );/**/ 
     3624                                fprintf(ioBAD,"DEBUG BAD DAY %li %li %li %li %.3e %.3e\n", 
     3625                                        iVibHi , iRotHi, iVibLo , iRotLo , heatone , coolone );/**/ 
    36023626                                if( ipHi==nLevels_per_elec[iElecHi]-1 && 
    36033627                                        ipLo==ipHi-1 ) 
     
    36053629                        } 
    36063630 
    3607                         /* derivative wrt temperature - assume exp wrt ground - this needs to be 
    3608                          * divided by square of temperature in wn -  
     3631                        /* derivative wrt temperature - assume exp wrt ground -  
     3632                         * this needs to be divided by square of temperature in wn -  
    36093633                         * done at end of loop */ 
    3610                         hmi.deriv_HeatH2Dexc_BigH2 +=  (float)(oneline *energy_wn[iElecHi][iVibHi][iRotHi]); 
     3634                        hmi.deriv_HeatH2Dexc_BigH2 +=  (float)(oneline * 
     3635                                energy_wn[iElecHi][iVibHi][iRotHi]); 
    36113636 
    36123637                        /* oneline is net heating, positive for heating, negative 
     
    36463671                 * nCountDebug is count to turn on debug prints,  
    36473672                 * nCountStop is count to abort code */ 
    3648                 long int nCountDebug = 630, 
    3649                         nCountStop = 700; 
    36503673                if(nCount>nCountDebug ) 
    36513674                { 
    36523675                        fprintf(ioQQQ, 
    3653                                 "DEBUG H2_Cooling A %li %15s, Te %.3e Heat(Xcol) %.2e H+/0 " 
     3676                                "DEBUG H2_Cooling A %li %15s, Te %.3e net Heat(Xcol) %.2e " 
     3677                                "heat %.2e cool %.2e H+/0 " 
    36543678                                "%.2e n(H2)%.3e Sol rat %.3e grn J1->0%.2e frac heat 1 line " 
    36553679                                "%.2e Hi(v,j)%li %li Lo(v,J)%li %li frac cool 1 line %.2e " 
     
    36583682                                phycon.te,  
    36593683                                hmi.HeatH2Dexc_BigH2 ,  
     3684                                H2_X_col_cool , 
     3685                                H2_X_col_heat , 
    36603686                                dense.xIonDense[ipHYDROGEN][1]/SDIV(dense.xIonDense[ipHYDROGEN][0]), 
    36613687                                hmi.H2_total, 
     
    37353761                                fprintf(ioQQQ,"\n");  
    37363762                } 
    3737                 ++nCount; 
    3738                 if( nCount > nCountStop ) 
    3739                 { 
    3740                         cdEXIT( EXIT_FAILURE ); 
    3741                 } 
    3742 #               if 0 
    3743                 else if( nCount > 630 ) 
    3744                 { 
    3745                         mole.nH2_TRACE = mole.nH2_trace_full; 
    3746                         trace.npsbug = 1; 
    3747                         trace.lgTrace = true; 
    3748                         geometry.nprint = 1; 
    3749                         iterations.IterPrnt[0] = 1; 
    3750                 } 
    3751 #               endif 
    37523763        } 
    37533764 
     
    37673778 
    37683779        { 
    3769                 enum {DEBUG_LOC=true }; 
     3780                enum {DEBUG_LOC=false }; 
    37703781                if( DEBUG_H2_COLL_X_HEAT && DEBUG_LOC &&  
    37713782                        (fabs(hmi.HeatH2Dexc_BigH2) > SMALLFLOAT) ) 
     
    38113822                                rate_dn_heat, 
    38123823                                rate_up_cool); 
     3824                        if( nCount>= nCountStop ) 
     3825                                cdEXIT(EXIT_FAILURE); 
    38133826                } 
    38143827        } 
     
    38363849        } 
    38373850 
     3851#       if 0 
    38383852        /* this can be noisy due to finite accuracy of solution, so take average with 
    38393853         * previous value */ 
     
    38483862                old_HeatH2Dexc = hmi.HeatH2Dexc_BigH2; 
    38493863        } 
    3850         { 
    3851                 enum {DEBUG_LOC=true}; 
    3852                 if( DEBUG_LOC && DEBUG_H2_COLL_X_HEAT ) 
    3853                 { 
    3854                         static long int kount = 0; 
    3855                         ++kount; 
    3856                         fprintf(ioQQQ,"DEBUG H2_cooling E %15s vib deex %li Te %.3e heat %.2e J7/5 %.4e J13/11 %.4e\n",  
     3864#       endif 
     3865        { 
     3866                enum {DEBUG_LOC=false}; 
     3867                if( DEBUG_LOC /*&& DEBUG_H2_COLL_X_HEAT*/ ) 
     3868                { 
     3869                        fprintf(ioQQQ,"DEBUG H2_cooling E %15s %c vib deex %li Te %.3e net heat %.3e cool %.3e heat %.3e\n",  
    38573870                                chRoutine , 
    3858                                 kount, 
     3871                                TorF(conv.lgSearch), 
     3872                                nCount, 
    38593873                                phycon.te, hmi.HeatH2Dexc_BigH2, 
     3874                                H2_X_col_cool , 
     3875                                H2_X_col_heat /*, 
    38603876                                H2_populations[0][0][7]/SDIV(H2_populations[0][0][5]) , 
    3861                                 H2_populations[0][0][13]/SDIV(H2_populations[0][0][11]) );/**/ 
    3862                         /*if( kount > 635 ) 
    3863                                 cdEXIT(EXIT_FAILURE);*/ 
     3877                                H2_populations[0][0][13]/SDIV(H2_populations[0][0][11])*/ ); 
     3878                        if( 0 && nCount > nCountStop ) 
     3879                        { 
     3880                                cdEXIT( EXIT_FAILURE ); 
     3881                        } 
    38643882                } 
    38653883        } 
     
    38743892        return; 
    38753893} 
    3876 #undef  DEBUG_H2_COLL_X_HEAT 
    38773894 
    38783895 
     
    40584075        return; 
    40594076} 
    4060  
    4061  
     4077/*lint +e668 possible bad pointer */ 
     4078/*lint +e802 possible bad pointer */ 
     4079 
     4080 
  • trunk/source/mole_h2_etc.cpp

    r1163 r1175  
    3939        iElecLo = 0; 
    4040 
    41         /* find rate (s-1) h2 dissoc into X continuum by Solomon process and 
     41        /* find rate (s-1) h2 dissociation into X continuum by Solomon process and 
    4242         * assign to the TH85 g and s states  
    4343         * these will go back into the chemistry network */ 
    4444 
    45         /* rates [s-1] for dissociation from s or g, into elec excited states   
     45        /* rates [s-1] for dissociation from s or g, into electronic excited states   
    4646         * followed by dissociation */ 
    4747        hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.; 
     
    5353 
    5454        /* at this point we have already evaluated the sum of the radiative rates out 
    55          * of the electronic excited states - this is H2_rad_rate_out[elec][vib][rot]  
     55         * of the electronic excited states - this is H2_rad_rate_out[electronic][vib][rot]  
    5656         * and this includes both decays into the continuum and bound states of X */ 
    5757 
     
    6464                        { 
    6565 
    66                                 /* loop over all lower levels within ground of X to find decay rates from H2g to continuum  
    67                                  * nEner_H2_ground is number of levels in H2g as determined by ENERGY_H2_STAR */ 
     66                                /* loop over all lower levels within ground of X to find  
     67                                 * decay rates from H2g to continuum  
     68                                 * nEner_H2_ground is number of levels in H2g as determined  
     69                                 * by ENERGY_H2_STAR */ 
    6870                                for( iplo=0; iplo < nEner_H2_ground; ++iplo ) 
    6971                                { 
     
    7375                                        if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo] ) 
    7476                                        { 
    75                                                 /* this is the rate {cm-3 s-1] that mole goes from lower level into electronic excited 
    76                                                  * states then into continuum */ 
     77                                                /* this is the rate {cm-3 s-1] that mole goes from  
     78                                                 * lower level into electronic excited states then  
     79                                                 * into continuum */ 
    7780                                                double rate_up_cont =  
    7881                                                        H2_populations[iElecLo][iVibLo][iRotLo]* 
     
    8083                                                        H2_dissprob[iElecHi][iVibHi][iRotHi] / H2_rad_rate_out[iElecHi][iVibHi][iRotHi]; 
    8184 
    82                                                 /* this is H2g up to exec excit then to continuum - cm-3 s-1 at this point */ 
     85                                                /* this is H2g up to excited then to continuum -  
     86                                                 * cm-3 s-1 at this point */ 
    8387                                                hmi.H2_Solomon_dissoc_rate_BigH2_H2g += rate_up_cont; 
    8488 
    85                                                 /* rate elec state decays into H2g */ 
     89                                                /* rate electronic state decays into H2g */ 
    8690                                                hmi.H2_Solomon_elec_decay_H2g +=  
    8791                                                        H2_populations[iElecHi][iVibHi][iRotHi]* 
     
    100104                                        if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo] ) 
    101105                                        { 
    102                                                 /* this is the rate {cm-3 s-1] that mole goes from lower level into elec excited 
    103                                                         * states then into continuum */ 
     106                                                /* this is the rate {cm-3 s-1] that mole goes from  
     107                                                 * lower level into electronic excited states then  
     108                                                 * into continuum */ 
    104109                                                double rate_up_cont =  
    105110                                                        H2_populations[iElecLo][iVibLo][iRotLo]* 
     
    107112                                                        H2_dissprob[iElecHi][iVibHi][iRotHi] / H2_rad_rate_out[iElecHi][iVibHi][iRotHi]; 
    108113 
    109                                                 /* this is H2s up to exec excit then to continuum - cm-3 s-1 at this point */ 
     114                                                /* this is H2s up to exec excited then to continuum  
     115                                                 * cm-3 s-1 at this point */ 
    110116                                                hmi.H2_Solomon_dissoc_rate_BigH2_H2s += rate_up_cont; 
    111117 
    112                                                 /* rate elec state decays into H2s */ 
     118                                                /* rate electronic state decays into H2s */ 
    113119                                                hmi.H2_Solomon_elec_decay_H2s +=  
    114120                                                        H2_populations[iElecHi][iVibHi][iRotHi]* 
     
    174180                iRotLoX = ipRot_H2_energy_sort[ip]; 
    175181                iVibLoX = ipVib_H2_energy_sort[ip]; 
    176                 /* now find all pumps up to elec excited states */ 
     182                /* now find all pumps up to electronic excited states */ 
    177183                /* sum over all electronic states, finding dissociation rate */ 
    178184                for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 
     
    184190                                        if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX] ) 
    185191                                        { 
    186                                                 /* this is the rate {cm-3 s-1] that mole goes from lower level into elec excited 
    187                                                  * states then into continuum */ 
     192                                                /* this is the rate {cm-3 s-1] that mole goes from  
     193                                                 * lower level into electronic excited states then  
     194                                                 * into continuum */ 
    188195                                                double rate_up_cont =  
    189196                                                        H2_populations[0][iVibLoX][iRotLoX]* 
     
    191198 
    192199                                                double decay_star = H2_rad_rate_out[iElecHi][iVibHi][iRotHi] - H2_dissprob[iElecHi][iVibHi][iRotHi]; 
    193                                                 /* loop over all other levels in H2g, subtracting their rate - remainder is rate into star  
    194                                                  * this is usually only a few levels */ 
     200                                                /* loop over all other levels in H2g, subtracting  
     201                                                 * their rate - remainder is rate into star, this is  
     202                                                 * usually only a few levels */ 
    195203                                                for( ipOther=0; ipOther < nEner_H2_ground; ++ipOther ) 
    196204                                                { 
     
    209217                                                /* MAX because may underflow to negative numbers is rates very large  
    210218                                                 * this is fraction that returns to H2s */ 
    211                                                 decay_star = MAX2(0., decay_star)/H2_rad_rate_out[iElecHi][iVibHi][iRotHi]
     219                                                decay_star = MAX2(0., decay_star)/SDIV(H2_rad_rate_out[iElecHi][iVibHi][iRotHi])
    212220                                                hmi.H2_H2g_to_H2s_rate_BigH2 += rate_up_cont*decay_star; 
    213221 
    214222                                        }/* end if line exists */ 
    215                                 }/* end loop rot elec excit */ 
    216                         }/* end loop vib elec excit */ 
    217                 }/* end loop elec elec excit */ 
     223                                }/* end loop rot electronic excited */ 
     224                        }/* end loop vib electronic excited */ 
     225                }/* end loop electronic electronic excited */ 
    218226        } 
    219227 
     
    222230 
    223231        DEBUG_EXIT( "H2_gs_rates()" ); 
    224  
    225232        return; 
    226233} 
     
    236243        DEBUG_ENTRY( "H2_zero_pops_too_low()" ); 
    237244 
    238         /* >>chng 05 jan 26, add this block to set populations to lte value */ 
     245        /* >>chng 05 jan 26, add this block to set populations to LTE value */ 
    239246        for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 
    240247        { 
     
    252259        } 
    253260        /* zero everything out - loop over all possible lines */ 
    254         /* >>chng 05 jan 26, set to lte values, since we still need to accumulate Lyman line 
     261        /* >>chng 05 jan 26, set to LTE values, since we still need to accumulate Lyman line 
    255262                * optical depths to have correct self-shielding when large h2 does come on */ 
    256263        for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 
     
    263270                        { 
    264271                                long int lim_elec_lo = 0; 
    265                                 /* now the lower levels */ 
    266                                 /* NB - iElecLo the lower electronic level is only X  
    267                                        * we don't consider excited elec to excit elec trans here, since we are only  
    268                                        * concerned with excited electronic levels as a photodissociation process 
    269                                        * code exists to relax this assumption - simply change following to iElecHi */ 
     272                                /* now the lower levels  
     273                                * NB - iElecLo the lower electronic level is only X  
     274                                * we don't consider excited electronic to excited electronic trans here, since we are only  
     275                                * concerned with excited electronic levels as a photodissociation process 
     276                                * code exists to relax this assumption - simply change following to iElecHi */ 
    270277                                for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 
    271278                                { 
    272                                         /* want to include all vib states in lower level if different elec level, 
    273                                         * but only lower vib levels if same elec level */ 
     279                                        /* want to include all vib states in lower level if different electronic level, 
     280                                        * but only lower vib levels if same electronic level */ 
    274281                                        long int nv = h2.nVib_hi[iElecLo]; 
    275282                                        if( iElecLo==iElecHi ) 
     
    287294 
    288295                                                                /* population of lower level */ 
    289                                                                 /* >>chng 05 jan 26, set to lte value */ 
     296                                                                /* >>chng 05 jan 26, set to LTE value */ 
    290297                                                                H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo.Pop =  
    291298                                                                        H2_populations[iElecLo][iVibLo][iRotLo]; 
     
    329336 
    330337        DEBUG_EXIT( "H2_zero_pops_too_low()" ); 
    331  
    332338        return; 
    333339} 
     
    380386                                                H2_stat[iElec][iVib][iRot] / part_fun; 
    381387                                        /*if( iElec==0 && iVib < 2) 
    382                                                 fprintf(ioQQQ,"DEBUG lte pop\t%i\t%i\t%e\n", 
     388                                                fprintf(ioQQQ,"DEBUG LTE pop\t%i\t%i\t%e\n", 
    383389                                                iVib,iRot,H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total);*/ 
    384390                                } 
     
    394400 
    395401        DEBUG_EXIT( "mole_H2_LTE()" ); 
    396  
    397402        return; 
    398  
    399403} 
    400404 
     
    404408        /* the order of the electronic states is 
    405409         * X, B, C+, C-, B', D+, and D- */ 
    406         /* this will be the number of vibration levels within each elec */ 
     410        /* this will be the number of vibration levels within each electronic */ 
    407411        /* number of vib states within electronic states from 
    408412         * >>refer      H2      energies        Abgrall, */ 
     
    413417        long int Jlowest_init[N_H2_ELEC] = {0 , 0 , 1  , 1 , 0 , 1 , 1 }; 
    414418 
    415         /* number of rotation levels within each elec - vib */ 
     419        /* number of rotation levels within each electronic - vib */ 
    416420        /*lint -e785 too few init for aggregate */ 
    417421        long int nRot_hi_init[N_H2_ELEC][50]= 
     
    450454        /* the order of the electronic states is 
    451455         * X, B, C+, C-, B', D+, and D- */ 
    452         /* this will be the number of vibration levels within each elec */ 
     456        /* this will be the number of vibration levels within each electronic */ 
    453457        /* number of vib states within electronic states from 
    454458         * >>refer      H2      energies        Abgrall, */ 
     
    486490 
    487491        /* number of times large molecules evaluated in this iteration, 
    488          * is false if never evaluated, on next evaluation will start with lte populations */ 
     492         * is false if never evaluated, on next evaluation will start with LTE populations */ 
    489493        nCallH2_this_iteration = 0; 
    490494 
     
    504508 
    505509        DEBUG_EXIT( "H2_Reset()" ); 
    506  
    507510        return; 
    508  
    509511} 
    510512 
     
    519521 
    520522        /* this is the smallest ratio of H2/H where we will bother with the large H2 molecule 
    521          * this value was chosen so that large mole used at very start of TH85 standard pdr
     523         * this value was chosen so that large mole used at very start of TH85 standard PDR
    522524         * NB - this appears in headinfo and must be updated there if changed here */ 
    523525        /* >>chng 03 jun 02, from 1e-6 to 1e-8 - in orion veil large H2 turned on half way 
     
    626628 
    627629        DEBUG_EXIT( "H2_Zero()" ); 
    628  
    629630        return; 
    630631} 
     
    834835                } 
    835836 
    836                 /*k = pow(10.,H2_He_coll_fit_par[init][final][0]+b2+c2) + pow(10.,H2_He_coll_fit_par[init][final][4]+f2);*/ 
    837837                sum1 = H2_He_coll_fit_par[init][final][0]+b2+c2; 
    838838                sum2 = H2_He_coll_fit_par[init][final][4]+f2; 
     
    846846                { 
    847847                        /* rate of 1e-6 is largest possible rate according to Phillip 
    848                         * Stancil email 07 may 27 */ 
     848                        * Stancil email 07 may 27 */ 
    849849                        fprintf(ioQQQ,"PROBLEM H2-He rate coefficient hit cap, upper index of %li" 
    850850                                " lower index of %li rate was %.2e resetting to 1e-6, Te=%e\n",