Changeset 1092

Show
Ignore:
Timestamp:
05/03/07 18:54:55 (21 months ago)
Author:
gary
Message:

dynamics.cpp - time option on hextra command
hextra.h - vars for time option for hextra command
zerologic.cpp - int vars for time option on hextra
parse_commands.cpp - time option on hextra.h
iso_cool.cpp - clean out dead code & comments,
prt_zone.cpp - fix comment
punch_do.cpp - add transmission coefficient on punch continuum

*.in fix trailing comments in class field

Location:
trunk
Files:
24 modified

Legend:

Unmodified
Added
Removed

  • trunk/source/date.h

    r1074 r1092  
    88 * in the morning so no date conflict happens */ 
    99#define YEAR    107 
    10 #define MONTH   3 
    11 #define DAY     28 
     10#define MONTH   4 
     11#define DAY     3 

  • trunk/source/dynamics.cpp

    r812 r1092  
    2222#include "radius.h" 
    2323#include "thirdparty.h" 
     24#include "hextra.h" 
    2425#include "rfield.h" 
    2526#include "trace.h" 
    … …  
    13401341 
    13411342                        /* this is set true on CORONAL XXX TIME INIT command, says to use constant 
    1342                          * temperature for forst N_INITIAL_RELAX iterations, then let run free */ 
     1343                         * temperature for first N_INITIAL_RELAX iterations, then let run free */ 
    13431344                        if( dynamics.lg_coronal_time_init  ) 
    13441345                        { 
    … …  
    14031404                                        ++nTimeVary; 
    14041405                        } 
    1405                         if( !nTimeVary ) 
    1406                         { 
    1407                                 fprintf(ioQQQ," DISASTER - there were no variable continua - " 
    1408                                         "put TIME option on at least one luminosity command.\n"); 
     1406 
     1407                        if( hextra.lgTurbHeatVaryTime ) 
     1408                        { 
     1409                                /* vary extra heating */ 
     1410                                hextra.TurbHeat = hextra.TurbHeatSave * time_continuum_scale; 
     1411                                fprintf(ioQQQ,"DEBUG TurbHeat vary new heat %.2e\n", 
     1412                                        hextra.TurbHeat); 
     1413                        } 
     1414                        else if( nTimeVary ) 
     1415                        { 
     1416                                for( i=0; i<rfield.nupper; ++i ) 
     1417                                { 
     1418                                        /* >>chng 06 mar 22, break into constant and time dependent parts */ 
     1419                                        rfield.flux[i] = rfield.flux_const[i] + 
     1420                                                rfield.flux_time[i]*time_continuum_scale; 
     1421                                        rfield.FluxSave[i] = rfield.flux[i]; 
     1422                                } 
     1423                        } 
     1424                        else 
     1425                        { 
     1426                                fprintf(ioQQQ," DISASTER - there were no variable continua " 
     1427                                        "or heat sources - put TIME option on at least one " 
     1428                                        "luminosity or hextra command.\n"); 
    14091429                                puts( "[Stop in DynaEndIter]" ); 
    14101430                                cdEXIT(EXIT_FAILURE); 
    1411                         } 
    1412  
    1413                         for( i=0; i<rfield.nupper; ++i ) 
    1414                         { 
    1415                                 /* >>chng 06 mar 22, break into constant and time dependent parts */ 
    1416                                 rfield.flux[i] = rfield.flux_const[i] + 
    1417                                         rfield.flux_time[i]*time_continuum_scale; 
    1418                                 rfield.FluxSave[i] = rfield.flux[i]; 
    14191431                        } 
    14201432                        /* this is heat radiated, after correction for change of H density in constant 
    … …  
    14251437                else 
    14261438                { 
    1427                         /* this branch, during initial stabilazation of soln */ 
     1439                        /* this branch, during initial stabilization of solution */ 
    14281440                        HeatInitial = 1.5 * pressure.PresGasCurr; 
    14291441                        HeatRadiated = 0.; 
    … …  
    14501462                                dynamics.lgRecom = true; 
    14511463                                /* set largest possible zone thickness to value on previous 
    1452                                         * iteratin when light was on - during recombination conditions 
     1464                                        * iteration when light was on - during recombination conditions 
    14531465                                        * become much more homogeneous and dr can get too large, 
    14541466                                        * crashing into H i-front */ 
    … …  
    14581470                        if( lgtime_dt_specified ) 
    14591471                        { 
    1460  
    14611472                                /* this is the new timestep */ 
    14621473                                fprintf(ioQQQ,"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element, 
    … …  
    19791990        timestep_init = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 
    19801991        /* this is log of timestep */ 
     1992        if( timestep_init > 30. ) 
     1993        { 
     1994                /* this number is large and will almost certainly crash */ 
     1995                fprintf(ioQQQ,"WARNING - the log of the initial time step is too " 
     1996                        "large, I shall probably crash.  The value was log t = %.2e\n", 
     1997                        timestep_init ); 
     1998                fflush(ioQQQ); 
     1999        } 
    19812000        timestep_init = pow( 10. , timestep_init ); 
    19822001        timestep = timestep_init; 

  • trunk/source/hextra.h

    r740 r1092  
    2121        float cryden_ov_background; 
    2222 
    23         /** this is default cosmic ray background  density and rate */ 
     23        /** this is default cosmic ray background density and rate */ 
    2424        float background_density; 
    2525        float background_rate; 
    … …  
    2727        /** these three set with hextra command, first the heating rate */ 
    2828        float TurbHeat,  
     29 
     30          /** save the initial value in case TurbHeat varies with time */ 
     31          TurbHeatSave , 
    2932 
    3033          /** the scale radius for the heating */ 
    … …  
    3437          turback; 
    3538 
     39        /** set true if extra heat varies with time in time dependent sims */ 
     40        bool lgTurbHeatVaryTime; 
     41 
    3642        }       hextra; 
    37  
    38  

  • trunk/source/hot_issues.txt

    r1075 r1092  
    11HOT ISSUES 
    2  
    3 cdColm has option to get level col den for large H2.  this is also repeated in cdH2_colden.  Doc for first is removed from hazy, and code for it should be removed too.  there is no reason to have two ways of doing exactly the same thing. 
    42 
    53ming gu dr data received 2007 feb 8, incorporate into cloudy 
    64 
    75table starburst - write claus 
     6 
     7table stars in hazy is wrong 
    88 
    99check e2 function at line 300 of service.cpp 
    … …  
    1313I am also not completely sure about what happens further down between lines 1145 and 1198. 
    1414 
     15NI pumping 
     16 
    1517what happened with e2 
    1618 
    1719uta code query, get out of conv_base 
     20 
     21rename lines.in 
    1822 
    1923hmi.H2_H2g_to_H2s_rate_used - make multiple occurances like rest 
    … …  
    2327 
    2428but make sure all "used" are saved and rest 
     29 
    2530 
    2631purify crash in func_trans_read after trans write 

  • trunk/source/iso_cool.cpp

    r1075 r1092  
    11/* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland 
    22 * For conditions of distribution and use see copyright notice in license.txt */ 
    3 /*iso_cool compute net cooling due to species in iso-sequences */ 
     3/*iso_cool compute net cooling due to hydrogenc atom species, ground state  
     4 * photoionization of hydrogenic species done in sumheat */ 
    45#include "cddefines.h" 
    56#include "physconst.h" 
    … …  
    1213#include "cooling.h" 
    1314#include "iso.h" 
     15/** \todo       2       - if pc lint ever fixes this bug in their product, remove this -e */ 
    1416/*lint -e661 Possible access of out-of-bounds pointer*/ 
    1517/*lint -e662  (Warning -- Possible creation of out-of-bounds pointer  */ 
    … …  
    3537          edenHCorr_IonAbund,  
    3638          edenIonAbund,  
    37           coll_ion_cool_one,  
     39          hex,  
    3840          HeatExcited, 
    39           *coll_ion_cool_array,  
    40           *hlone , 
     41          *hexx,  /* will become save arrays dimd nhlvl+1*/ 
     42          *hlone ,/* will become save arrays dimd nhlvl+1*/ 
    4143          ThinCoolingCaseB,  
    4244          ThinCoolingSum, 
    … …  
    8587        DEBUG_ENTRY( "iso_cool()" ); 
    8688 
     89        /** \todo       2       this routine dumps all heating and cooling into only a very few lables, 
     90         * it would be best to break these out into individual labels that show element, 
     91         * iso sequence, and agent */ 
     92 
    8793        /* validate the incoming data */ 
    8894        ASSERT( nelem >= ipISO ); 
    8995        ASSERT( ipISO <NISO ); 
    9096        ASSERT( nelem < LIMELM ); 
    91         /* local number of levels may be less than malloced number if continuum 
    92          * has been lowered due to high density */ 
    93         ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] ); 
    94  
    95         /* for zero abundance or if element has been turned off we need to  
    96          * set some produced variables to zero */ 
     97 
     98        if( ipISO==0 ) 
     99        { 
     100                /* >>chng 06 aug 28, continuum lowering now on, allow less than here. */ 
     101                ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] ); 
     102        } 
     103        else if( ipISO==1 ) 
     104        { 
     105                ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] ); 
     106        } 
     107        else 
     108                TotalInsanity(); 
     109 
     110        /* for zero abundance we need to set some produced variables to zero */ 
     111        /* >>chng 04 nov 13, add check on element off */ 
    97112        if( dense.xIonDense[nelem][nelem+1-ipISO] <= 0. ||  
    98113                /* >>chng 04 dec 07, add test for recombined species having zero abundance, 
    … …  
    131146        } 
    132147 
    133         /* create some space, these go to numLevels_local instead of _max 
    134          * since continuum may have been lowered by density */ 
    135         coll_ion_cool_array = (double *)MALLOC( (unsigned)(iso.numLevels_local[ipISO][nelem])*sizeof(double) ); 
     148        /* create some space */ 
     149        /* >>chng 06 aug 17, all these should go to numLevels_local instead of _max. */ 
     150        hexx = (double *)MALLOC( (unsigned)(iso.numLevels_local[ipISO][nelem])*sizeof(double) ); 
    136151        hlone = (double *)MALLOC( (unsigned)(iso.numLevels_local[ipISO][nelem])*sizeof(double) ); 
    137152        SavePhotoHeat = (double *)MALLOC( (unsigned)(iso.numLevels_local[ipISO][nelem])*sizeof(double) ); 
    … …  
    139154        SaveRadRecCool = (double *)MALLOC( (unsigned)(iso.numLevels_local[ipISO][nelem])*sizeof(double) ); 
    140155 
     156        /* compute the total cooling due to the hydrogen atom, and its derivative 
     157         * halfte is 1/(2T) */ 
     158 
    141159        /*********************************************************************** 
    142160         *                                                                     * 
    … …  
    145163         ***********************************************************************/ 
    146164 
    147         /* will be net collisional ionization cooling, units erg/cm^3/s */ 
    148         coll_ion_cool_one = 0.; 
     165        /* hex will be net collisional ionization cooling, units erg/cm^3/s */ 
     166        hex = 0.; 
    149167        dhexdt = 0.; 
    150168 
    151         /* collisional ionization cooling minus three body heating  
    152          * depending on how topoff is done, highest level can have large population 
    153          * and its coupling to continuum can be large, at various times code  
    154          * had to ignore effects of very highest level, but starting in mid 
    155          * 20006 all levels have been included */ 
    156         for( n=0; n < iso.numLevels_local[ipISO][nelem]; ++n ) 
     169        /* collisional ionization cooling, three body heating  
     170         * highest level is fastest collision by far, its coupling to 
     171         * continuum can overwhelm heating, so only count nearly up to the top */ 
     172        /* >>chng 02 apr 25, from numLevels to numLevels-1 - departure coef for highest can be large */ 
     173        /* >>chng 06 aug 17, upper limit from numLevels_max to _local with continuum lowering 
     174         * treated consistently */ 
     175        /* >>chng 06 sep 1, take "-1" take out, negate 02 apr 25 change. */ 
     176        for( n=0; n < iso.numLevels_local[ipISO][nelem]/*-1*/; ++n ) 
    157177        { 
    158178                /* total collisional ionization cooling less three body heating */ 
    159                 coll_ion_cool_array[n] =  
     179                hexx[n] =  
    160180                  iso.xIsoLevNIonRyd[ipISO][nelem][n]*iso.ColIoniz[ipISO][nelem][n]*collider* 
    161181                  (iso.Pop2Ion[ipISO][nelem][n] -iso.PopLTE[ipISO][nelem][n]*dense.eden); 
    162                 coll_ion_cool_one += coll_ion_cool_array[n]; 
    163                 /*fprintf(ioQQQ,"DEBUG coll_ion_cool_one %.2e %2li %2li %2li %.2e %.4e %.4e %.3f\n", 
    164                         coll_ion_cool_one, 
     182                hex += hexx[n]; 
     183                /*fprintf(ioQQQ,"DEBUG hex %.2e %2li %2li %2li %.2e %.4e %.4e %.3f\n", 
     184                        hex, 
    165185                        ipISO , nelem , n , 
    166                         coll_ion_cool_array[n] ,  
     186                        hexx[n] ,  
    167187                        iso.Pop2Ion[ipISO][nelem][n] , 
    168188                        iso.PopLTE[ipISO][nelem][n]*dense.eden , 
    … …  
    171191 
    172192                /* the derivative of the cooling  
    173                  * need extra factor of temp of 1 Ryd since div by square of temp in Ryd */ 
    174                 dhexdt += coll_ion_cool_array[n]*(iso.xIsoLevNIonRyd[ipISO][nelem][n]/TE1RYD/POW2(phycon.te_ryd) -  
     193                 * need extra factor of temp of 1 ryd since div by square of temp in ryd */ 
     194                dhexdt += hexx[n]*(iso.xIsoLevNIonRyd[ipISO][nelem][n]/TE1RYD/POW2(phycon.te_ryd) -  
    175195                  thermal.halfte); 
    176196        } 
     197        /*if( ipISO==1 && nelem==1 ) 
     198        { 
     199                n = 85; 
     200                fprintf(ioQQQ,"DEBUG hexx %.2e %.2e\n", 
     201                iso.Pop2Ion[ipISO][nelem][n] , 
     202                iso.PopLTE[ipISO][nelem][n]*dense.eden ); 
     203        }*/ 
     204 
     205        /* convert to physical units, need to convert ryd to ergs,  
     206         * and bring to density per vol not per ion */ 
     207        /* >>chng 05 dec 28, leave hex in hexx units for simple debugging */ 
     208        /*hex *= EN1RYD*dense.xIonDense[nelem][nelem+1-ipISO]; 
     209        dhexdt *= EN1RYD*dense.xIonDense[nelem][nelem+1-ipISO];*/ 
    177210 
    178211        /* save net collisional ionization cooling less H-3 body heating 
    179          * convert to physical units, need to convert Ryd to ergs,  
     212         * for inclusion in printout  
     213         * convert to physical units, need to convert ryd to ergs,  
    180214         * and bring to density per vol not per ion */ 
    181         iso.coll_ion[ipISO][nelem] = coll_ion_cool_one * EN1RYD*dense.xIonDense[nelem][nelem+1-ipISO]; 
     215        iso.coll_ion[ipISO][nelem] = hex * EN1RYD*dense.xIonDense[nelem][nelem+1-ipISO]; 
    182216        /* create a meaningful label */ 
    183217        sprintf(chLabel , "IScion%2s%2s" , elementnames.chElementSym[ipISO] ,  
    … …  
    187221        CoolAdd(chLabel,(float)nelem,MAX2(0.,iso.coll_ion[ipISO][nelem])); 
    188222 
    189         /* derivative wrt temp  
    190          * convert to physical units, need to convert Ryd to ergs,  
     223        /* deriv wrt temp  
     224         * convert to physical units, need to convert ryd to ergs,  
    191225         * and bring to density per vol not per ion */ 
    192226        thermal.dCooldT += dhexdt*EN1RYD*dense.xIonDense[nelem][nelem+1-ipISO]; 
    193227 
    194         /* heating[0][3] is all three body heating, opposite of collisional  
    195          * ionization cooling, 
     228        /* heating[0][3] is all three body heating, opposite of coll ion cooling, 
    196229         * would be unusual for this to be non-zero since would require excited 
    197230         * state departure coefficients to be greater than unity */ 
    198231        thermal.heating[0][3] += MAX2(0.,-iso.coll_ion[ipISO][nelem]); 
    199         { 
    200                 /* debug block printing individual contributors to this process */ 
     232        /*if( ipISO==1 && nelem==1 ) 
     233                fprintf(ioQQQ,"DEBUG %li %2li %.2e %.4e\n",  
     234                        ipISO , nelem ,  
     235                        -iso.coll_ion[ipISO][nelem]/ MAX2( SMALLFLOAT , thermal.htot ) , 
     236                        phycon.te );*/ 
     237        { 
     238                /*@-redef@*/ 
    201239                enum {DEBUG_LOC=false}; 
     240                /*@+redef@*/ 
    202241                if( DEBUG_LOC  && nelem==1 && ipISO==1 ) 
    203242                { 
    … …  
    209248                                for( n=0; n < (iso.numLevels_local[ipISO][nelem] - 1); n++ ) 
    210249                                { 
    211                                         if( coll_ion_cool_array[n] / SDIV( coll_ion_cool_one ) > 0.1 ) 
     250                                        if( hexx[n] / SDIV( hex ) > 0.1 ) 
    212251                                                fprintf(ioQQQ," %2li %.1e", 
    213252                                                n, 
    214                                                 coll_ion_cool_array[n]/ SDIV( coll_ion_cool_one ) ); 
     253                                                hexx[n]/ SDIV( hex ) ); 
    215254                                } 
    216255                                fprintf(ioQQQ,"\n"); 
    … …  
    236275        edenHCorr_IonAbund = collider*dense.xIonDense[nelem][nelem+1-ipISO]; 
    237276 
     277        /* >>chng 02 apr 23 */ 
     278        /* new logic, for any iso seq */ 
    238279        /* now do case b sum to compare with exact value below */ 
    239280        iso.RadRecCool[ipISO][nelem] = 0.; 
    … …  
    254295                 thin = iso.RadRecomb[ipISO][nelem][0][ipRecRad]* 
    255296                        /* arg is the scaled temperature, T * n^2 / Z^2,  
    256                          * n is principle quantum number, Z is charge, 1 for H */ 
     297                         * n is prin quant number, Z is charge, 1 for H */ 
    257298                        HCoolRatio(  
    258299                        phycon.te * POW2( (double)iso.quant_desig[ipISO][nelem][0].n / (double)(nelem+1-ipISO) ))* 
    … …  
    260301                        phycon.te * BOLTZMANN; 
    261302        } 
    262         /* the cooling, corrected for optical depth */ 
     303        /* the cooling, corrected for optical dephs */ 
    263304        SaveRadRecCool[0] = iso.RadRecomb[ipISO][nelem][0][ipRecNetEsc] * thin; 
    264305        /* this is now total free-bound cooling */ 
    … …  
    266307 
    267308        /* radiative recombination cooling for all excited states */ 
     309        /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 
    268310        for( n=1; n < iso.numLevels_local[ipISO][nelem]; n++ ) 
    269311        { 
    … …  
    271313                 thin = iso.RadRecomb[ipISO][nelem][n][ipRecRad]* 
    272314                        /* arg is the scaled temperature, T * n^2 / Z^2,  
    273                          * n is principle quantum number, Z is charge, 1 for H */ 
     315                         * n is prin quant number, Z is charge, 1 for H */ 
    274316                        HCoolRatio(  
    275317                        phycon.te * POW2( (double)iso.quant_desig[ipISO][nelem][n].n / (double)(nelem+1-ipISO) ))* 
    … …  
    277319                        phycon.te * BOLTZMANN; 
    278320 
    279                 /* the cooling, corrected for optical depth */ 
     321                /* the cooling, corrected for optical dephs */ 
    280322                SaveRadRecCool[n] = iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc] * thin * edenIonAbund; 
    281323                /* this is now total free-bound cooling */ 
    … …  
    286328        } 
    287329        { 
    288                 /* debug block for state specific recombination cooling */ 
     330                /*@-redef@*/ 
    289331                enum {DEBUG_LOC=false}; 
     332                /*@+redef@*/ 
    290333                if( DEBUG_LOC  ) 
    291334                { 
    … …  
    302345        } 
    303346 
    304         /* Case b sum of optically thin radiative recombination cooling.  
    305          * add any remainder to the sum from above - high precision is needed  
    306          * to get STE result to converge close to equilibrium - only done for 
    307          * H-like ions where exact result is known */ 
     347        /* following is expression for case b sum of  
     348         * optically thin radiative recombination cooling - we will need to add any remainder to 
     349         * the sum from above - high precision is needed to get STE result to converge 
     350         * close to equilibrium - as in lte.in and induchhe.in */ 
     351        /* >>chng 02 apr 26, do not print for he-like - they are always negative at low t */ 
    308352        if( ipISO == ipH_LIKE ) 
    309353        { 
    … …  
    350394        RecCoolExtra = MAX2(0., RecCoolExtra ); 
    351395 
    352         /* add error onto total - this is significant for approach to STE */ 
     396        /* add error onto total - this is significant for approach to ste */ 
     397        /* >>chng 06 aug 17, should have numLevels_local instead of _max. */ 
    353398        iso.RadRecCool[ipISO][nelem] += RecCoolExtra* edenIonAbund *iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecNetEsc]; 
    354399 
    … …  
    360405         ***********************************************************************/ 
    361406 
    362         /* photoionization of excited levels */ 
     407        /* photoionization of excited levels 
     408         * HPhotHeat(n,nelem) is photoionization heating due to level n, 
     409         * evaluated in iso_photo */ 
    363410        HeatExcited = 0.; 
    364411        ipbig = -1000; 
     412        /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 
    365413        for( n=1; n < (iso.numLevels_local[ipISO][nelem] - 1); ++n ) 
    366414        { 
    … …  
    375423        } 
    376424        { 
    377                 /* debug block for heating due to photo of each n */ 
     425                /* this block to debug heating due to photo of var n */ 
     426                /*@-redef@*/ 
    378427                enum {DEBUG_LOC=false}; 
     428                /*@+redef@*/ 
    379429                if( DEBUG_LOC  && ipISO==0 && nelem==0  && nzone > 700) 
    380430                { 
    … …  
    402452        } 
    403453 
    404         /* FreeBnd_net_Cool_Rate is net cooling due to recombination  
    405          * RadRecCool is total radiative recombination cooling sum to all levels, 
    406          * with n>=2 photoionization heating subtracted */ 
     454        /* HFreBndNet is net cooling due to  
     455         * HRadRecCool is total radiative recombination cooling sum to all levels, 
     456         * net is this less n>=2 photoionization heating */ 
    407457        iso.FreeBnd_net_Cool_Rate[ipISO][nelem] = iso.RadRecCool[ipISO][nelem] - HeatExcited; 
    408458        /*fprintf(ioQQQ,"DEBUG Hn2\t%.3e\t%.3e\n", 
    … …  
    410460                HeatExcited/SDIV(thermal.htot));*/ 
    411461 
    412         /* heating[0][1] is all excited state photoionization heating from ALL  
    413          * species, this is set to zero in CoolEvaluate before loop where this  
    414          * is called. */ 
     462        /* heating[1][0] is all excited state photoionization heating from ALL species, 
     463         * this is set to zero in CoolEvaluate before loop where this is called. 
     464         * >>chng 96 jan 25, was -HeatExcited */ 
    415465        thermal.heating[0][1] += MAX2(0.,-iso.FreeBnd_net_Cool_Rate[ipISO][nelem]); 
    416466 
    … …  
    441491        } 
    442492        { 
    443                 /* print rec cool, induced rec cool, photo heat */ 
     493                /* print rec cool, induc rec cool, photo heat */ 
     494                /*@-redef@*/ 
    444495                enum {DEBUG_LOC=false}; 
     496                /*@+redef@*/ 
    445497                if( DEBUG_LOC && ipISO==0 && nelem==5 /**/ ) 
    446498                { 
    … …  
    476528        CoolAdd(chLabel,(float)nelem,iso.RecomInducCool_Rate[ipISO][nelem]); 
    477529 
    478         /* find total collisional energy exchange due to bound-bound */ 
     530        /* lines are remainder of this routine */ 
     531        /* this will be the only loop in the final code 
     532         * find total collisional energy exchange */ 
    479533        iso.xLineTotCool[ipISO][nelem] = 0.; 
    480534        dCdT_all = 0.; 
    … …  
    483537        nhi_heat_max = -1; 
    484538 
    485         /* loop does not include highest levels - their population may be 
    486          * affected by topoff */ 
    487539        for( ipLo=0; ipLo < iso.numLevels_local[ipISO][nelem]-2; ipLo++ ) 
    488540        { 
    489541                for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipISO][nelem]-1; ipHi++ ) 
    490542                { 
    491                         /* collisional energy exchange between ipHi and ipLo - net cool */ 
    492543              &nb