mole_drive.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*mole_drive - public routine, calls mole_solve to converge molecules */
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "hmi.h"
#include "conv.h"
#include "doppvel.h"
#include "thermal.h"
#include "gammas.h"
#include "iso.h"
#include "opacity.h"
#include "phycon.h"
#include "physconst.h"
#include "radius.h"
#include "rfield.h"
#include "secondaries.h"
#include "timesc.h"
#include "trace.h"
#include "thermal.h"
#include "called.h"
#include "hcmap.h"
#include "coolheavy.h"
#include "mole.h"
#include "mole_priv.h"
#include "dynamics.h"
#include "h2.h"
#include "co.h"
#include "ionbal.h"
 
/* this is the error tolerance for reporting non-convergence */
static const double MOLETOLER = 0.10;
 
STATIC void mole_effects(void);
STATIC void mole_h_rate_diagnostics(void);
STATIC void mole_ion_trim(void);
STATIC void mole_update_limiting_reactants();
 
/*mole_drive main driver for heavy molecular equilibrium routines      */
void mole_drive(void)
{
        DEBUG_ENTRY( "mole_drive()" );
 
        mole_update_species_cache();  /* Update densities of species controlled outside the chemical network */
 
        mole_update_limiting_reactants();
        
        mole_update_rks();
        
        double error = mole_solve();    
        
        bool lgConverged = (error < MOLETOLER); 
 
        if( !lgConverged )
        {
                // should something be done with this?
        }
 
        mole_ion_trim();
 
        mole_effects();
        
        return;
}
 
STATIC void mole_ion_trim(void)
{
        for (ChemAtomList::iterator atom = unresolved_atom_list.begin();
                  atom != unresolved_atom_list.end(); ++atom)
        {
                const long int nelem=(*atom)->el->Z-1;
                if( !dense.lgElmtOn[nelem] )
                        continue;
 
                for (long int ion=0;ion<nelem+2;ion++) 
                {
                        if ((*atom)->ipMl[ion] != -1)
                        {
                                if (dense.IonLow[nelem] > ion )
                                {
                                        if( dense.xIonDense[nelem][ion] > ( ionbal.trimlo * dense.gas_phase[nelem] ) )
                                                dense.IonLow[nelem] = ion;
                                        else
                                        {
                                                dense.xIonDense[nelem][dense.IonLow[nelem]] += dense.xIonDense[nelem][ion];
                                                dense.xIonDense[nelem][ion] = 0.;
                                        }
                                }
 
                                if (dense.IonHigh[nelem] < ion )
                                {
                                        if( dense.xIonDense[nelem][ion] > ( ionbal.trimhi * dense.gas_phase[nelem] ) )
                                                dense.IonHigh[nelem] = ion;
                                        else
                                        {
                                                dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion];
                                                dense.xIonDense[nelem][ion] = 0.;
                                        }
                                }
                        }
                }
        }
}
 
void mole_update_sources(void)
{
        DEBUG_ENTRY( "mole_update_sources()" );
 
        mole_update_species_cache();
        mole_eval_sources(mole_global.num_total);
}
 
STATIC void mole_effects()
{
        double 
                plte;
 
        DEBUG_ENTRY( "mole_effects()" );
 
        mole_eval_sources(mole_global.num_total);
 
        co.CODissHeat = (realnum)(mole.findrate("PHOTON,CO=>C,O")*1e-12);
 
        thermal.heating[0][9] = co.CODissHeat;
 
        /* total H2 - all forms */
        hmi.H2_total = findspecieslocal("H2*")->den + findspecieslocal("H2")->den;
        hmi.HD_total = findspecieslocal("HD")->den;
        /* first guess at ortho and para densities */
        h2.ortho_density = 0.75*hmi.H2_total;
        h2.para_density = 0.25*hmi.H2_total;
        {
                hmi.H2_total_f = (realnum)hmi.H2_total;
                h2.ortho_density_f = (realnum)h2.ortho_density;
                h2.para_density_f = (realnum)h2.para_density;
        }
 
        /* save rate H2 is destroyed units s-1 */
        /* >>chng 05 mar 18, TE, add terms - 
                total destruction rate is: dest_tot = n_H2g/n_H2tot * dest_H2g + n_H2s/n_H2tot * dest_H2s */
        /* as reactions that change H2s to H2g and vice versa are not counted destruction processes, the terms c[ipMH2g][ipMH2s] *
           and c[ipMH2s][ipMH2g], which have a different sign than [ipMH2g][ipMH2g] and [ipMH2s][ipMH2s], have to be added      */
        hmi.H2_rate_destroy = mole.sink_rate_tot("H2");
 
        /* establish local timescale for H2 molecule destruction */
        if(hmi.H2_rate_destroy > SMALLFLOAT )
        {
                /* units are now seconds */
                timesc.time_H2_Dest_here = 1./hmi.H2_rate_destroy;
        }
        else
        {
                timesc.time_H2_Dest_here = 0.;
        }
        
        /* local timescale for H2 formation 
         * both grains and radiative attachment 
         * >>chng 08 mar 15 rjrw -- ratio formation rate to density of _H2_ so comparable with destruction rate,
         * use full rates not rate coefficients so have correct units */
        timesc.time_H2_Form_here = (mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn"));
        /* timescale is inverse of this rate */
        if( timesc.time_H2_Form_here > SMALLFLOAT )
        {
                /* units are now seconds */
                timesc.time_H2_Form_here =  hmi.H2_total/timesc.time_H2_Form_here;
        }
        else
        {
                timesc.time_H2_Form_here = 0.;
        }
 
 
        /* heating due to H2 dissociation */
 
        /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
        /*hmi.HeatH2Dish_TH85 = (float)(1.36e-23*findspecieslocal("H2")->den*h2esc*hmi.UV_Cont_rel2_Habing_TH85_depth);*/
        /* >>chng 04 feb 07, more general to express heating in terms of the assumed
         * photo rates - the 0.25 was obtained by inverting A8 & A9 of TH85 to find that
         * there are 0.25 eV per dissociative pumping, ie, 10% of total 
         * this includes both H2g and H2s - TH85 say just ground but they include
         * process for both H2 and H2s - as we did above - both must be in
         * heating term */
        /* >>chng 05 mar 11, TE, old had used H2_Solomon_dissoc_rate_used, which was only
         * H2g.  in regions where Solomon process is fast, H2s has a large population
         * and the heating rate was underestimated. */
        /* >>chng 05 jun 23, 
         * >>chng 05 dec 05, TE, modified to approximate the heating better for the
         * new approximation */
        /* >>chng 00 nov 25, explicitly break out this heat agent */
        /* 2.6 eV of heat per deexcitation, consider difference
         * between deexcitation (heat) and excitation (cooling) */
        /* >>chng 04 jan 27, code moved here and simplified */
        /* >>chng 05 jul 10, GS*/ 
        /*  average collisional rate for H2* to H2g calculated from big H2, GS*/
        
        /* TH85 dissociation heating - this is ALWAYS defined for reference
         * may be output for comparison with other rates*/
        hmi.HeatH2Dish_TH85 = 0.25 * EN1EV *
                (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
                 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
        
        /* TH85 deexcitation heating */
        hmi.HeatH2Dexc_TH85 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
                                                                                                 - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
        /* this is derivative wrt temperature, only if counted as a cooling source */
        hmi.deriv_HeatH2Dexc_TH85 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_TH85)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
        
        if( hmi.chH2_small_model_type == 'H' )
        {
                /* Burton et al. 1990 */
                hmi.HeatH2Dish_BHT90 = 0.25 * EN1EV *
                        (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
                         hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
                
                /* Burton et al. 1990 heating */
                hmi.HeatH2Dexc_BHT90 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H")) 
                                                                                                                - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
                /* this is derivative wrt temperature, only if counted as a cooling source */
                hmi.deriv_HeatH2Dexc_BHT90 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BHT90)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
        }
        else if( hmi.chH2_small_model_type == 'B')
        {
                /* Bertoldi & Draine */
                hmi.HeatH2Dish_BD96 = 0.25 * EN1EV *
                        (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
                         hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den);
                /* Bertoldi & Draine heating, same as TH85 */
                hmi.HeatH2Dexc_BD96 = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
                                                                                                         - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12;
                /* this is derivative wrt temperature, only if counted as a cooling source */
                hmi.deriv_HeatH2Dexc_BD96 = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_BD96)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
        }
        else if(hmi.chH2_small_model_type == 'E')
        {
                /* heating due to dissociation of H2
                 * >>chng 05 oct 19, TE, define new approximation for the heating due to the destruction of H2
                 *      use this approximation for the specified cloud parameters, otherwise
                 * use limiting cases for 1 <= G0, G0 >= 1e7, n >= 1e7, n <= 1 */
                
                double log_density, 
                        f1, f2,f3, f4, f5;
                static double log_G0_face = -1;
                static double k_f4;
                
                
                /* test for G0 
                 * this is a constant so only do it in zone 0 */
                if( !nzone )
                {
                        if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.) 
                        { 
                                log_G0_face = 0.;
                        }
                        else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7) 
                        { 
                                log_G0_face = 7.;
                        }
                        else 
                        { 
                                log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face); 
                        }
                        /*>>chng 06 oct 24 TE change Go face for spherical geometry*/
                        log_G0_face /= radius.r1r0sq;
                }
 
                /* test for density */
                if(dense.gas_phase[ipHYDROGEN] <= 1.) 
                { 
                        log_density = 0.; 
                }
                else if(dense.gas_phase[ipHYDROGEN] >= 1e7) 
                { 
                        log_density = 7.; 
                }
                else 
                { 
                        log_density = log10(dense.gas_phase[ipHYDROGEN]); 
                }
                
                f1 = 0.15 * log_density + 0.75;
                f2 = -0.5 * log_density + 10.;
                
                hmi.HeatH2Dish_ELWERT = 0.25 * EN1EV *  f1 * 
                        (hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den +
                         hmi.H2_Solomon_dissoc_rate_used_H2s * findspecieslocal("H2*")->den ) + 
                        f2 * secondaries.x12tot * EN1EV * hmi.H2_total;
                
                /*fprintf( ioQQQ, "f1: %.2e, f2: %.2e,heat Solomon: %.2e",f1,f2,hmi.HeatH2Dish_TH85);*/
                
                
                /* heating due to collisional deexcitation within X of H2
                 * >>chng 05 oct 19, TE, define new approximation for the heating due to the destruction of H2
                 *      use this approximation for the specified cloud parameters, otherwise
                 * use limiting cases for 1 <= G0, G0 >= 1e7, n >= 1e7, n <= 1 */
                
                /* set value of k_f4 by testing on value of G0 */
                if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.) 
                { 
                        log_G0_face = 0.;
                }
                else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7) 
                { 
                        log_G0_face = 7.;
                }
                else 
                { 
                        log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face); 
                }
                /* 06 oct 24, TE introduce effects of spherical geometry */
                log_G0_face /= radius.r1r0sq;
                
                /* terms only dependent on G0_face */
                k_f4 = -0.25 * log_G0_face + 1.25;
 
                /* test for density */
                if(dense.gas_phase[ipHYDROGEN] <= 1.) 
                { 
                        log_density = 0.; 
                        f4 = 0.; 
                }
                else if(dense.gas_phase[ipHYDROGEN] >= 1e7) 
                { 
                        log_density = 7.; 
                        f4 = pow2(k_f4) * pow( 10. , 2.2211 * log_density - 29.8506);
                }
                else 
                { 
                        log_density = log10(dense.gas_phase[ipHYDROGEN]); 
                        f4 = pow2(k_f4) * pow( 10., 2.2211 * log_density - 29.8506);
                }
                
                f3 = MAX2(0.1, -4.75 * log_density + 24.25);
                f5 = MAX2(1.,0.95 * log_density - 1.45) * 0.2 * log_G0_face;
                
                hmi.HeatH2Dexc_ELWERT = ((mole.findrate("H2*,H2=>H2,H2") + mole.findrate("H2*,H=>H2,H"))
                        - (mole.findrate("H2,H2=>H2*,H2") + mole.findrate("H2,H=>H2*,H"))) * 4.17e-12 * f3 + 
                        f4 * (mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[0]].den/dense.gas_phase[ipHYDROGEN]) +
                        f5 * secondaries.x12tot * EN1EV * hmi.H2_total;
                
                if(log_G0_face == 0.&& dense.gas_phase[ipHYDROGEN] > 1.) 
                        hmi.HeatH2Dexc_ELWERT *= 0.001 / dense.gas_phase[ipHYDROGEN];
                
                /* >>chng 06 oct 24, TE introduce effects of spherical geometry */
                /*if(radius.depth/radius.rinner >= 1.0) */
                hmi.HeatH2Dexc_ELWERT /= POW2(radius.r1r0sq);
                
                /* this is derivative wrt temperature, only if counted as a cooling source */
                hmi.deriv_HeatH2Dexc_ELWERT = (realnum)(MAX2(0.,-hmi.HeatH2Dexc_ELWERT)* ( 30172. * thermal.tsq1 - thermal.halfte ) );
                
                /*fprintf( ioQQQ, "\tf3: %.2e, f4: %.2e, f5: %.2e, heat coll dissoc: %.2e\n",f3,f4,f5,hmi.HeatH2Dexc_TH85);*/
        }
        /* end Elwert branch for photo rates */
        else
                TotalInsanity();
 
 
        if( findspecieslocal("H-")->den > 0. && hmi.rel_pop_LTE_Hmin > 0. )
        {
                hmi.hmidep = (double)findspecieslocal("H-")->den/ SDIV( 
                        ((double)dense.xIonDense[ipHYDROGEN][0]*dense.eden*hmi.rel_pop_LTE_Hmin));
        }
        else
        {
                hmi.hmidep = 1.;
        }
 
        /* this will be net volume heating rate, photo heat - induc cool */
        hmi.hmihet = hmi.HMinus_photo_heat*findspecieslocal("H-")->den - hmi.HMinus_induc_rec_cooling;
 
        /* H2+  +  HNU  =>  H+  + H */
        t_phoHeat photoHeat;
        (void) GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.,&photoHeat); /* Now used entirely for side-effect.  Yuk */
        hmi.h2plus_heat = photoHeat.HeatNet*findspecieslocal("H2+")->den;
 
        /* departure coefficient for H2 defined rel to n(1s)**2
         * (see Littes and Mihalas Solar Phys 93, 23) */
        plte = (double)dense.xIonDense[ipHYDROGEN][0] * hmi.rel_pop_LTE_H2g * (double)dense.xIonDense[ipHYDROGEN][0];
        if( plte > 0. )
        {
                hmi.h2dep = findspecieslocal("H2")->den/plte;
        }
        else
        {
                hmi.h2dep = 1.;
        }
 
        /* departure coefficient of H2+ defined rel to n(1s) n(p)
         * sec den was HI before 85.34 */
        plte = (double)dense.xIonDense[ipHYDROGEN][0]*hmi.rel_pop_LTE_H2p*(double)dense.xIonDense[ipHYDROGEN][1];
        if( plte > 0. )
        {
                hmi.h2pdep = findspecieslocal("H2+")->den/plte;
        }
        else
        {
                hmi.h2pdep = 1.;
        }
 
        /* departure coefficient of H3+ defined rel to N(H2+) N(p) */
        if( hmi.rel_pop_LTE_H3p > 0. )
        {
                hmi.h3pdep = findspecieslocal("H3+")->den/hmi.rel_pop_LTE_H3p;
        }
        else
        {
                hmi.h3pdep = 1.;
        }
 
        mole_h_rate_diagnostics();
 
        return;
}
#if defined(__HP_aCC)
#pragma OPTIMIZE OFF
#pragma OPTIMIZE ON
#endif
STATIC void mole_h_rate_diagnostics(void)
{
        int i;
        double 
                rate; 
 
        /* identify dominant H2 formation process */
        {
                /* following should be set true to identify H2 formation and destruction processes */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && (nzone>50) )
                {
                        double createsum ,create_from_Hn2 , create_3body_Ho, create_h2p, 
                                create_h3p, create_h3pe, create_grains, create_hminus;
                        double destroysum, destroy_hm ,destroy_solomon ,destroy_2h ,destroy_hp,
                                destroy_h,destroy_hp2,destroy_h3p;
 
                        create_from_Hn2 = mole.findrate("H,H=>H2");                     /* H(n=2) + H(n=1) -> H2 */
                        create_3body_Ho = mole.findrate("H,H,H=>H2,H");
                        create_h2p = mole.findrate("H,H2+=>H+,H2*");
                        create_h3p = mole.findrate("H,H3+=>H2,H2+");
                        create_h3pe = mole.findrate("e-,H3+=>H2,H"); 
                        create_grains = mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn");
                        create_hminus = (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"));
 
                        createsum = create_from_Hn2 + create_3body_Ho + create_h2p +
                                create_h3p + create_h3pe + create_grains + create_hminus;
 
                        fprintf(ioQQQ,"H2 create zone\t%.2f \tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                fnzone,
                                createsum,
                                create_hminus   / createsum,
                                create_from_Hn2 / createsum,
                                create_3body_Ho / createsum,
                                create_h2p      / createsum,
                                create_h3p      / createsum,
                                create_h3pe     / createsum,
                                create_grains   / createsum     );
 
                        /* all h2 molecular hydrogen destruction processes */
                        /* >>chng 04 jan 28, had wrong Boltzmann factor,
                         * caught by gargi Shaw */
                        destroy_hm = (mole.findrate("H2,e-=>H,H-")+mole.findrate("H2*,e-=>H,H-"));
                        /*destroy_hm2 = eh2hhm;*/
                        destroy_solomon = hmi.H2_Solomon_dissoc_rate_used_H2g * findspecieslocal("H2")->den;
                        destroy_2h = (mole.findrate("H2,e-=>H,H,e-")+mole.findrate("H2*,e-=>H,H,e-"));
                        destroy_hp = mole.findrate("H2,H+=>H3+");
                        destroy_h = mole.findrate("H,H2=>H,H,H");
                        destroy_hp2 = mole.findrate("H2,H+=>H2+,H");
                        destroy_h3p = mole.findrate("H2,H3+=>H2,H2+,H");
                        destroysum = destroy_hm + /*destroy_hm2 +*/ destroy_solomon + destroy_2h + 
                                destroy_hp+ destroy_h+ destroy_hp2+ destroy_h3p;
 
                        fprintf(ioQQQ,"H2 destroy\t%.3f \t%.2e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                fnzone,
                                destroysum,
                                destroy_hm / destroysum ,
                                destroy_solomon / destroysum ,
                                destroy_2h / destroysum ,
                                destroy_hp / destroysum ,
                                destroy_h / destroysum ,
                                destroy_hp2 / destroysum ,
                                destroy_h3p / destroysum );
 
                }
        }
 
        {
                /* following should be set true to identify H- formation and destruction processes */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && (nzone>140) )
                {
                        double create_from_Ho,create_3body_Ho,create_batach,destroy_photo,
                                destroy_coll_heavies,destroy_coll_electrons,destroy_Hattach,destroy_fhneut,
                                destsum , createsum;
 
                        create_from_Ho = mole.findrate("H,e-=>H-,PHOTON");
                        create_3body_Ho = mole.findrate("H,e-,e-=>H-,e-");
                        /* total formation is sum of g and s attachment */
                        create_batach = (mole.findrate("H2,e-=>H,H-") + mole.findrate("H2*,e-=>H,H-"));
                        destroy_photo = mole.findrate("H-,PHOTON=>H,e-");
                        destroy_coll_heavies = 0.;
                        for( i=ipLITHIUM; i < LIMELM; i++ )
                        {
                                char react[32];
                                if (dense.lgElmtOn[i]) 
                                {
                                        sprintf(react,"H-,%s=>H,%s",mole_global.list[unresolved_atom_list[i]->ipMl[1]]->label.c_str(),unresolved_atom_list[i]->label().c_str() );
                                        destroy_coll_heavies += mole.findrate(react);
                                }
                        }
                        destroy_coll_electrons = mole.findrate("H-,e-=>H-,e-,e-");
                        destroy_Hattach = (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"));
                        destroy_fhneut = mole.findrate("H-,H+=>H,H");
 
                        destsum = destroy_photo + destroy_coll_heavies + destroy_coll_electrons + 
                                destroy_Hattach + destroy_fhneut;
                        fprintf(ioQQQ,"H- destroy zone\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", 
                                fnzone,
                                phycon.te,
                                destsum,
                                destroy_photo/destsum , 
                                destroy_coll_heavies/destsum,
                                destroy_coll_electrons/destsum, 
                                destroy_Hattach/destsum,
                                destroy_fhneut/destsum );
 
                        createsum = create_from_Ho+create_3body_Ho+create_batach;
                        fprintf(ioQQQ,"H- create\t%.2f\tTe\t%.4e\tsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                fnzone,
                                phycon.te,
                                createsum,
                                dense.eden,
                                create_from_Ho/createsum,
                                create_3body_Ho/createsum,
                                create_batach/createsum);
                }
        }
        {
                /* following should be set true to print populations */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC )
                {
                        /* these are the raw results */
                        fprintf( ioQQQ, " MOLE raw; hi\t%.2e" , dense.xIonDense[ipHYDROGEN][0]);
                        for( i=0; i < mole_global.num_calc; i++ )
                        {
                                fprintf( ioQQQ, "\t%-4.4s\t%.2e", mole_global.list[i]->label.c_str(), mole.species[i].den );
                        }
                        fprintf( ioQQQ, " \n" );
                }
        }
 
        if( trace.lgTrace && trace.lgTraceMole )
        {
                /* these are the raw results */
                fprintf( ioQQQ, " raw; " );
                for( i=0; i < mole_global.num_calc; i++ )
                {
                        fprintf( ioQQQ, " %-4.4s:%.2e", mole_global.list[i]->label.c_str(), mole.species[i].den );
                }
                fprintf( ioQQQ, " \n" );
 
        }
 
        /* option to print rate H2 forms */
        /* trace.lgTraceMole is trace molecules option,
         * punch htwo */
        if( (trace.lgTrace && trace.lgTraceMole) )
        {
                double H2_rate_create = mole.source_rate_tot("H2") + mole.source_rate_tot("H2*");
 
                if( H2_rate_create > SMALLFLOAT )
                {
                        fprintf( ioQQQ, " Create H2, rate=%10.2e grain;%5.3f hmin;%5.3f bhedis;%5.3f h2+;%5.3f hmi.radasc:%5.3f hmi.h3ph2p:%5.3f hmi.h3petc:%5.3f\n", 
                                H2_rate_create, 
                                (mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn"))/H2_rate_create, 
                                (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-"))/H2_rate_create, 
                                mole.findrate("H,H,H=>H2,H")/H2_rate_create, 
                                mole.findrate("H,H2+=>H+,H2*")/H2_rate_create, 
                                mole.findrate("H,H=>H2")/H2_rate_create, 
                                mole.findrate("H,H3+=>H2,H2+")/H2_rate_create, 
                                mole.findrate("H2,H3+=>H2,H2+,H")/H2_rate_create );
                }
                else
                {
                        fprintf( ioQQQ, " Create H2, rate=0\n" );
                }
        }
 
        /* this is H2+ */
        if( trace.lgTrace && trace.lgTraceMole )
        {
                rate = mole.findrate("H2,H+=>H2+,H") + mole.findrate("H,H+,e-=>H2+,e-") + 
                  mole.findrate("H,H3+=>H2,H2+") + mole.findrate("H2,H3+=>H2,H2+,H");
                if( rate > 1e-25 )
                {
                        fprintf( ioQQQ, " Create H2+, rate=%10.2e hmi.rh2h2p;%5.3f b2pcin;%5.3f hmi.h3ph2p;%5.3f hmi.h3petc+;%5.3f\n", 
                          rate, mole.findrate("H2,H+=>H2+,H")/rate, 
                                                         mole.findrate("H,H+,e-=>H2+,e-")/rate, mole.findrate("H,H3+=>H2,H2+")/rate, mole.findrate("H2,H3+=>H2,H2+,H")/rate );
                }
                else
                {
                        fprintf( ioQQQ, " Create H2+, rate=0\n" );
                }
        }
 
        if( trace.lgTrace && trace.lgTraceMole )
        {
 
                double destroy_coll_heavies = 0.;
 
                for( i=ipLITHIUM; i < LIMELM; i++ )
                {
                        char react[32];
                        if (dense.lgElmtOn[i]) 
                        {
                                sprintf(react,"H-,%s=>H,%s",mole_global.list[unresolved_atom_list[i]->ipMl[1]]->label.c_str(),unresolved_atom_list[i]->label().c_str() );
                                destroy_coll_heavies += mole.findrate(react);
                        }
                }
                fprintf( ioQQQ, " MOLE, Dep Coef, H-:%10.2e H2:%10.2e H2+:%10.2e\n", 
                  hmi.hmidep, hmi.h2dep, hmi.h2pdep );
                fprintf( ioQQQ, "     H- creat: Rad atch%10.3e Induc%10.3e bHneut%10.2e 3bod%10.2e b=>H2%10.2e N(H-);%10.2e b(H-);%10.2e\n", 
                                                 mole.findrate("H,e-=>H-,PHOTON"), hmi.HMinus_induc_rec_rate*findspecieslocal("H-")->den, mole.findrate("H,H=>H-,H+"), 
                                                 mole.findrate("H,e-,e-=>H-,e-"), 
                                                 mole.findrate("H2,e-=>H,H-"), findspecieslocal("H-")->den, hmi.hmidep );
 
                fprintf( ioQQQ, "     H- destr: Photo;%10.3e mut neut%10.2e e- coll ion%10.2e =>H2%10.2e x-ray%10.2e p+H-%10.2e\n", 
                                                 mole.findrate("H-,PHOTON=>H,e-"), destroy_coll_heavies, 
                                                 mole.findrate("H-,e-=>H-,e-,e-"), 
                                                 (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")), iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*findspecieslocal("H-")->den, 
                                                 mole.findrate("H-,H+=>H,H") );
                fprintf( ioQQQ, "     H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n", 
                                                 hmi.hmihet, hmi.HMinus_induc_rec_cooling, hmi.hmicol );
        }
 
        /* identify creation and destruction processes for H2+ */
        if( trace.lgTrace && trace.lgTraceMole )
        {
                rate = findspecieslocal("H2+")->snk;
                if( rate != 0. )
                {
                        rate *= findspecieslocal("H2+")->den;
                        if( rate > 0. )
                        {
                                fprintf( ioQQQ, 
                                                                 " Destroy H2+: rate=%10.2e e-;%5.3f phot;%5.3f hard gam;%5.3f H2col;%5.3f h2phhp;%5.3f pion;%5.3f bh2h2p:%5.3f\n", 
                                                                 rate, mole.findrate("H2+,e-=>H,H+,e-")/rate, mole.findrate("H2+,PHOTON=>H,H+")/rate, 
                                                                 mole.findrate("H2+,CRPHOT=>H,H+")/rate, mole.findrate("H2,H2+=>H,H3+")/rate, mole.findrate("H2+,H2=>H,H+,H2")/rate, 
                                                                 mole.findrate("H2+,H+=>H,H+,H+")/rate, mole.findrate("H,H2+=>H+,H2*")/rate );
 
                                fprintf( ioQQQ, 
                                                                 " Create  H2+: rate=%.2e HII HI;%.3f Col H2;%.3f HII H2;%.3f HI HI;%.3f\n", 
                                                                 rate, 
                                                                 mole.findrate("H+,H=>H2+,PHOTON")/rate, 
                                                                 mole.findrate("H2,CRPHOT=>H2+,e-")/rate, 
                                                                 mole.findrate("H2,H+=>H2+,H")/rate, 
                                                                 mole.findrate("H,H+,e-=>H2+,e-")/rate );
                        }
                        else
                        {
                                fprintf( ioQQQ, " Create  H2+: rate= is zero\n" );
                        }
                }
        }
 
        {
                /* following should be set true to print populations */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC )
                {
                        fprintf(ioQQQ,"mole bugg\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                fnzone,
                                iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].gamnc,
                                hmi.HMinus_photo_rate,
                                findspecieslocal("H2")->den , 
                                findspecieslocal("H-")->den ,
                                findspecieslocal("H+")->den);
                }
        }
 
        {
                /*@-redef@*/
                /* this debug print statement compares H2 formation through grn vs H- */
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && nzone>140 )
                {
                        fprintf(ioQQQ," debuggggrn grn\t%.2f\t%.3e\t%.3e\tfrac\t%.3e\tH-\t%.3e\t%.3e\tfrac\t%.3e\t%.3e\t%.3e\t%.3e\n",
                                fnzone ,
                                mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn") , 
                                hmi.H2_forms_grains+hmi.H2star_forms_grains ,
                                                        mole.findrate("H,H,grn=>H2*,grn")/SDIV(mole.findrate("H,H,grn=>H2*,grn")+mole.findrate("H,H,grn=>H2,grn")),
                                                        (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),
                                hmi.H2star_forms_hminus+hmi.H2_forms_hminus,
                                                        frac_H2star_hminus(),
                                                        (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),findspecieslocal("H")->den,findspecieslocal("H-")->den
                                );
                }
        }
 
        {
                /*@-redef@*/
                /* often the H- route is the most efficient formation mechanism for H2,
                 * will be through rate called mole_global.list[unresolved_atom_list[ipHYDROGEN]->ipMl]->den*hmi.assoc_detach
                 * this debug print statement is to trace h2 oscillations */
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && nzone>140/*&& iteration > 1*/)
                {
                        /* rapid increase in H2 density caused by rapid increase in hmi.rel_pop_LTE_H2g */
                        fprintf(ioQQQ,
                                                        "hmi.assoc_detach_backwards_grnd\t%.2f\t%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 
                                                        /* total forward rate */
                                                        fnzone,
                                                        phycon.te, 
                                                        dense.eden,
                                                        (mole.findrate("H,H-=>H2,e-")+mole.findrate("H,H-=>H2*,e-")),
                                                        mole.findrate("H2,e-=>H,H-"), 
                                                        mole.findrate("H2*,e-=>H,H-"), 
                                                        mole.findrate("H,e-=>H-,PHOTON"),
                                                        hmi.HMinus_induc_rec_rate*findspecieslocal("H-")->den,
                                                        mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[0]].den,
                                                        mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[1]].den,
                                                        findspecieslocal("H-")->den,
                                                        hmi.H2_total,
                                                        hmi.rel_pop_LTE_Hmin,
                                                        hmi.rel_pop_LTE_H2g,
                                                        hmi.rel_pop_LTE_H2s
                                );
                }
        }
 
        if( (trace.lgTrace && trace.lgTraceMole) )
        {
                if( hmi.H2_rate_destroy > SMALLFLOAT )
                {
                        fprintf( ioQQQ, 
                          " H2 destroy rate=%.2e DIS;%.3f bat;%.3f h2dis;%.3f photoionize_rate;%.3f h2h2p;%.3f E-h;%.3f hmi.h2hph3p;%.3f sec;%.3f\n", 
                                                         hmi.H2_rate_destroy, 
                                                         hmi.H2_Solomon_dissoc_rate_used_H2g / hmi.H2_rate_destroy, 
                                                         mole.findrate("H2,e-=>H,H-") / (hmi.H2_total*hmi.H2_rate_destroy), 
                                                         mole.findrate("H,H2=>H,H,H") / (hmi.H2_total*hmi.H2_rate_destroy), 
                                                         h2.photoionize_rate / hmi.H2_rate_destroy, 
                                                         mole.findrate("H2,H+=>H2+,H") / (hmi.H2_total* hmi.H2_rate_destroy), 
                                                         (mole.findrate("H2,e-=>H,H,e-")+mole.findrate("H2*,e-=>H,H,e-")) / (hmi.H2_total* hmi.H2_rate_destroy), 
                                                         mole.findrate("H2,H+=>H3+") / (hmi.H2_total* hmi.H2_rate_destroy) ,
                                                         secondaries.csupra[ipHYDROGEN][0]*2.02 / hmi.H2_rate_destroy
                          );
                }
                else
                {
                        fprintf( ioQQQ, " Destroy H2: rate=0\n" );
                }
        }
 
        {
                /* following should be set true to print populations */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC )
                {
                        if( DEBUG_LOC && (nzone > 570) ) 
                        {
                                fprintf(ioQQQ,"Temperature %g\n",phycon.te);
                                fprintf(ioQQQ," Net mol ion rate [%g %g] %g\n",mole.source[ipHYDROGEN][1],mole.sink[ipHYDROGEN][1],
                                                                mole.source[ipHYDROGEN][1]-mole.sink[ipHYDROGEN][1]*mole.species[unresolved_atom_list[ipHYDROGEN]->ipMl[1]].den);
                        }
                }
        }
}
 
STATIC void mole_update_limiting_reactants()
{
        DEBUG_ENTRY( "mole_update_limiting_reactants()" );
 
        /* largest fraction of atoms in molecules */
        for( long i=0; i<mole_global.num_calc; ++i )
        {
                mole.species[i].xFracLim = 0.;
                mole.species[i].nAtomLim = -1;
                for (molecule::nAtomsMap::iterator atom = mole_global.list[i]->nAtom.begin();
                          atom != mole_global.list[i]->nAtom.end(); ++atom)
                {
                        long nelem = atom->first->el->Z-1;
                        if( dense.lgElmtOn[nelem])
                        {
                                realnum dens_elemsp = (realnum)( mole.species[i].den * atom->second * atom->first->frac );
                                if (mole.species[i].xFracLim*dense.gas_phase[nelem] < dens_elemsp ) 
                                {
                                        mole.species[i].xFracLim = dens_elemsp/dense.gas_phase[nelem];
                                        mole.species[i].nAtomLim = nelem;
                                }
                        }
                }
                //fprintf( ioQQQ, "DEBUGGG mole_update_limiting_reactants %li\t%s\t%i\t%.12e\n", i, mole_global.list[i]->label.c_str(), mole.species[i].nAtomLim, mole.species[i].xFracLim );
        }
        
        return;
}