Changeset 1656 for branches/newmole/source

Timestamp:

12/10/07 10:00:15 (12 months ago)

Author:

rjrw

Message:

Remove some leftover code which over-wrote H- FB rate output.

Move mole_effects to mole_co_drive, a more logical home for it.

Location:

branches/newmole/source

Files:

• mole.h (modified) (1 diff)
• mole_co_drive.cpp (modified) (3 diffs)
• mole_newton_step.cpp (modified) (2 diffs)

branches/newmole/source/mole.h

@@ -247 +247 @@
                                                                                                  bool lgHeader, bool lgData, double depth);
 
-extern void mole_effects(void);
-
 extern void total_molecule_elems(realnum total[LIMELM]);
 

branches/newmole/source/mole_co_drive.cpp

@@ -7 +7 @@
 #include "hmi.h"
 #include "conv.h"
+#include "gammas.h"
+#include "mole_co_atom.h"
+#include "opacity.h"
 #include "phycon.h"
+#include "physconst.h"
+#include "radius.h"
+#include "secondaries.h"
+#include "timesc.h"
 #include "trace.h"
 #include "thermal.h"
@@ -33 +40 @@
 STATIC void MolecSetup (
         char *chReason );
+STATIC void mole_effects(void);
 
 static void CO_null_results(void);
@@ -401 +409 @@
         return;
 }
+STATIC void mole_effects()
+{
+        double 
+                h2phet,
+                plte;
+        realnum abundan;
+
+        double *b, **c;
+        long int i;
+
+        DEBUG_ENTRY( "mole_effects()" );
+
+        b = (double *)MALLOC((size_t)mole.num_total*sizeof(double));
+        c = (double **)MALLOC((size_t)mole.num_total*sizeof(double *));
+        c[0] = (double *)MALLOC((size_t)mole.num_total*
+                                                                                                        mole.num_total*sizeof(double));
+        for(i=1;i<mole.num_total;i++)
+        {
+                c[i] = c[i-1]+mole.num_total;
+        }
+
+        mole_eval_balance(mole.num_total,b,c);  
+
+        co.CODissHeat = (realnum)(CO_findrate("PHOTON,CO=>C,O")*1e-12);
+
+        thermal.heating[0][9] = co.CODissHeat;
+
+        /* the real multi-level model molecule */
+        abundan = (realnum) findspecies("CO")->hevmol;
+        /* IonStg and nelem were set to 2, 0 in makelevlines */
+        ASSERT( C12O16Rotate[0].Hi->IonStg < LIMELM+2 );
+        ASSERT( C12O16Rotate[0].Hi->nelem-1 < LIMELM+2 );
+        dense.xIonDense[ C12O16Rotate[0].Hi->nelem-1][C12O16Rotate[0].Hi->IonStg-1] = abundan;
+
+        CO_PopsEmisCool(&C12O16Rotate, nCORotate,abundan , "12CO",
+                &CoolHeavy.C12O16Rot,&CoolHeavy.dC12O16Rot );
+        /*if( nzone>400 )
+                fprintf(ioQQQ,"DEBUG co cool\t%.2f\t%.4e\t%li\t%.4e\t%.4e\n",
+                        fnzone,CoolHeavy.C12O16Rot,nCORotate,abundan,phycon.te );*/
+
+        abundan = (realnum) findspecies("CO")->hevmol/co.C12_C13_isotope_ratio;
+        /* IonStg and nelem were set to 3, 0 in makelevlines */
+        ASSERT( C13O16Rotate[0].Hi->IonStg < LIMELM+2 );
+        ASSERT( C13O16Rotate[0].Hi->nelem-1 < LIMELM+2 );
+        dense.xIonDense[ C13O16Rotate[0].Hi->nelem-1][C13O16Rotate[0].Hi->IonStg-1] = abundan;
+
+        CO_PopsEmisCool(&C13O16Rotate, nCORotate,abundan ,"13CO",
+                &CoolHeavy.C13O16Rot,&CoolHeavy.dC13O16Rot );
+
+        /* 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      */
+        {
+                int ipMH2, ipMH2s;
+                ipMH2 = findspecies("H2")->index;
+                ipMH2s = findspecies("H2*")->index;
+                hmi.H2_rate_destroy = 0;
+                if (ipMH2 != -1)
+                        hmi.H2_rate_destroy += -findspecies("H2")->hevmol * c[ipMH2][ipMH2];
+                if (ipMH2s != -1)
+                        hmi.H2_rate_destroy += -findspecies("H2*")->hevmol * c[ipMH2s][ipMH2s];
+                if (ipMH2 != -1 && ipMH2s != -1) 
+                        hmi.H2_rate_destroy += 
+                                - findspecies("H2")->hevmol * c[ipMH2][ipMH2s]
+                                - findspecies("H2*")->hevmol * c[ipMH2s][ipMH2];
+                hmi.H2_rate_destroy /= SDIV(hmi.H2_total);
+        }
+
+        /* establish local timescale for H2 molecule destruction */
+        if(findspecies("H2")->index != -1 && -c[findspecies("H2")->index][findspecies("H2")->index] > SMALLFLOAT )
+        {
+                /* units are now seconds */
+                timesc.time_H2_Dest_here = -1./c[findspecies("H2")->index][findspecies("H2")->index];
+        }
+        else
+        {
+                timesc.time_H2_Dest_here = 0.;
+        }
+        
+        /* local timescale for H2 formation 
+         * both grains and radiative attachment */
+        timesc.time_H2_Form_here = (CO_findrk("H,H,grnh=>H2*")+CO_findrk("H,H,grnh=>H2")) * 
+                /* this corrects for fact that we the timescale for H2 to form from an atomic gas.
+                 * The rate becomes very small when gas is fully molecular, and ratio of total hydrogen
+                 * to atomic hydrogen corrections for this. */
+                dense.gas_phase[ipHYDROGEN]/SDIV(dense.xIonDense[ipHYDROGEN][0]);
+        /* timescale is inverse of this rate */
+        if( timesc.time_H2_Form_here > SMALLFLOAT )
+        {
+                /* units are now seconds */
+                timesc.time_H2_Form_here = 1./timesc.time_H2_Form_here;
+        }
+        else
+        {
+                timesc.time_H2_Form_here = 0.;
+        }
+
+
+        /* total H2 - all forms */
+        hmi.H2_total = (realnum) findspecies("H2*")->hevmol + findspecies("H2")->hevmol;
+        /* first guess at ortho and para densities */
+        h2.ortho_density = 0.75*hmi.H2_total;
+        h2.para_density = 0.25*hmi.H2_total;
+
+        /* NB the first index must be kept parallel with nelem and ionstag in
+         * H2Lines EmLine struc,
+         * since that struc expects to find the abundances here */
+        /* >>chng 04 feb 19, had been ipMH2g, chng to total */
+        dense.xIonDense[LIMELM+2][0] = hmi.H2_total;
+
+
+        /* heating due to H2 dissociation */
+
+        /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
+        /*hmi.HeatH2Dish_TH85 = (float)(1.36e-23*findspecies("H2")->hevmol*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 * findspecies("H2")->hevmol +
+                 hmi.H2_Solomon_dissoc_rate_used_H2s * findspecies("H2*")->hevmol);
+        
+        /* TH85 deexcitation heating */
+        hmi.HeatH2Dexc_TH85 = ((CO_findrate("H2*,H2=>H2,H2") + CO_findrate("H2*,H=>H2,H"))
+                                                                                                 - (CO_findrate("H2,H2=>H2*,H2") + CO_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 * findspecies("H2")->hevmol +
+                         hmi.H2_Solomon_dissoc_rate_used_H2s * findspecies("H2*")->hevmol);
+                
+                /* Burton et al. 1990 heating */
+                hmi.HeatH2Dexc_BHT90 = ((CO_findrate("H2*,H2=>H2,H2") + CO_findrate("H2*,H=>H2,H")) 
+                                                                                                                - (CO_findrate("H2,H2=>H2*,H2") + CO_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 * findspecies("H2")->hevmol +
+                         hmi.H2_Solomon_dissoc_rate_used_H2s * findspecies("H2*")->hevmol);
+                /* Bertoldi & Draine heating, same as TH85 */
+                hmi.HeatH2Dexc_BD96 = ((CO_findrate("H2*,H2=>H2,H2") + CO_findrate("H2*,H=>H2,H"))
+                                                                                                         - (CO_findrate("H2,H2=>H2*,H2") + CO_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 * findspecies("H2")->hevmol +
+                         hmi.H2_Solomon_dissoc_rate_used_H2s * findspecies("H2*")->hevmol ) + 
+                        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 = pow(k_f4,2) * pow( 10. , 2.2211 * log_density - 29.8506);
+                }
+                else 
+                { 
+                        log_density = log10(dense.gas_phase[ipHYDROGEN]); 
+                        f4 = pow(k_f4,2) * 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 = ((CO_findrate("H2*,H2=>H2,H2") + CO_findrate("H2*,H=>H2,H"))
+                                                                                                                 - (CO_findrate("H2,H2=>H2*,H2") + CO_findrate("H2,H=>H2*,H"))) * 4.17e-12 * f3 + 
+                        f4 * (COmole[element_list[ipHYDROGEN]->ipMl[0]]->hevmol/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( h2.lgH2ON  && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
+                 {
+                 deexc_htwo = hmi.Average_collH2;
+                 deexc_hneut = hmi.Average_collH;
+                 }
+                 else
+                 {
+                 deexc_htwo = (1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6.;
+                 deexc_hneut =  (1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.;
+                 } */
+        
+        /* Leiden hacks try to turn off H2*, which is all unphysical.  do not include heating
+         * due to H2 deexcitation since H2* is bogus */
+        /* >>chng 05 aug 12, do not turn off vibrational heating when Leiden hack is in place
+         * other codes included heating but did not include H2s on chemistry 
+         hmi.HeatH2Dexc_TH85 *= hmi.lgLeiden_Keep_ipMH2s;*/
+        /*fprintf(ioQQQ,
+                "DEBUG hmole H2 deex heating:\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
+                fnzone,
+                hmi.HeatH2Dexc_TH85,
+                thermal.htot,
+                findspecies("H2*")->hevmol,
+                H2star_deexcit, 
+                findspecies("H2")->hevmol,
+                H2star_excit,
+                Hmolec_old[ipMH2g],
+                Hmolec_old[ipMH],
+                phycon.te);*/
+        
+        /* >>chng 05 jul 11, TE, each term must have unit cm-3 s-1,
+         * more terms added */
+        hmi.H2_rate_create = 
+                CO_findrate("H,H,grnh=>H2*")+
+                CO_findrate("H,H,grnh=>H2") + 
+                CO_findrate("H,H-=>H2,e-")+
+                CO_findrate("H,H-=>H2*,e-") +
+                CO_findrate("H,H,H=>H2,H") + 
+                CO_findrate("H,H2+=>H+,H2*") + 
+                CO_findrate("H,H=>H2") + 
+                CO_findrate("H,H3+=>H2,H2+") + 
+                CO_findrate("H2+,H-=>H2,H") +
+                CO_findrate("H,H,H2=>H2,H2") + 
+                CO_findrate("H3+,H-=>H2,H,H") +
+                CO_findrate("H3+,H-=>H2,H2") +
+                CO_findrate("H2,H3+=>H+,H2,H2") +
+                CO_findrate("e-,H3+=>H2,H") +
+                CO_findrate("H3+,PHOTON=>H2,H+") +
+                CO_findrate("H,CH=>C,H2") +
+                CO_findrate("H,CH+=>C+,H2") +
+                CO_findrate("H,CH2=>CH,H2") +
+                CO_findrate("H,CH3+=>CH2+,H2") +
+                CO_findrate("H,OH=>O,H2") +
+                CO_findrate("H-,HCO+=>CO,H2") +
+                CO_findrate("H-,H3O+=>H2O,H2") +
+                CO_findrate("H-,H3O+=>OH,H2,H") +
+                CO_findrate("H+,CH2=>CH+,H2") +
+                CO_findrate("H+,SiH=>Si+,H2") +
+                CO_findrate("H2+,CH=>CH+,H2") + 
+                CO_findrate("H2+,CH2=>CH2+,H2") + 
+                CO_findrate("H2+,CO=>CO+,H2") + 
+                CO_findrate("H2+,H2O=>H2O+,H2") + 
+                CO_findrate("H2+,O2=>O2+,H2") + 
+                CO_findrate("H2+,OH=>OH+,H2") + 
+                CO_findrate("H3+,C=>CH+,H2") + 
+                CO_findrate("H3+,CH=>CH2+,H2") + 
+                CO_findrate("H3+,CH2=>CH3+,H2") + 
+                CO_findrate("H3+,OH=>H2O+,H2") + 
+                CO_findrate("H3+,H2O=>H3O+,H2") + 
+                CO_findrate("H3+,CO=>HCO+,H2") + 
+                CO_findrate("H3+,O=>OH+,H2") + 
+                CO_findrate("H3+,SiH=>SiH2+,H2") + 
+                CO_findrate("H3+,SiO=>SiOH+,H2") +
+                CO_findrate("H,CH3=>CH2,H2") +
+                CO_findrate("H,CH4+=>CH3+,H2") +
+                CO_findrate("H,CH5+=>CH4+,H2") +
+                CO_findrate("H2+,CH4=>CH3+,H2,H") + 
+                CO_findrate("H2+,CH4=>CH4+,H2") + 
+                CO_findrate("H3+,CH3=>CH4+,H2") + 
+                CO_findrate("H3+,CH4=>CH5+,H2") + 
+                CO_findrate("H+,CH4=>CH3+,H2") +        
+                CO_findrate("H+,HNO=>NO+,H2") +
+                CO_findrate("H+,HS=>S+,H2") +
+                CO_findrate("H,HS+=>S+,H2") +
+                CO_findrate("H3+,NH=>NH2+,H2") + 
+                CO_findrate("H3+,NH2=>NH3+,H2") + 
+                CO_findrate("H3+,NH3=>NH4+,H2") + 
+                CO_findrate("H3+,CN=>HCN+,H2") + 
+                CO_findrate("H3+,NO=>HNO+,H2") + 
+                CO_findrate("H3+,S=>HS+,H2") + 
+                CO_findrate("H3+,CS=>HCS+,H2") + 
+                CO_findrate("H3+,NO2=>NO+,OH,H2") + 
+                CO_findrate("H2+,NH=>NH+,H2") + 
+                CO_findrate("H2+,NH2=>NH2+,H2") + 
+                CO_findrate("H2+,NH3=>NH3+,H2") + 
+                CO_findrate("H2+,CN=>CN+,H2") + 
+                CO_findrate("H2+,HCN=>HCN+,H2") + 
+                CO_findrate("H2+,NO=>NO+,H2") +
+                CO_findrate("H3+,Cl=>HCl+,H2")+
+                CO_findrate("H3+,HCl=>H2Cl+,H2")+
+                CO_findrate("H2+,C2=>C2+,H2")+  
+                CO_findrate("H-,NH4+=>NH3,H2")+
+                CO_findrate("H3+,HCN=>HCNH+,H2");
+
+        /*>>chng 05 jul 01, GS, h2s ionization by cosmic ray added*/
+        /* >>chng 05 jun 29, TE, used in new punch H2 destruction file*/
+        hmi.CR_reac_H2g = CO_findrk("H2,CRPHOT=>H2+,e-") + CO_findrk("H2,CRPHOT=>H,H+,e-") 
+                + CO_findrk("H2,CRPHOT=>H,H") + CO_findrk("H2,CRPHOT=>H+,H-") + CO_findrk("H2,CRPHOT=>H+,H,e-");
+        hmi.CR_reac_H2s = CO_findrk("H2*,CRPHOT=>H2+,e-") + CO_findrk("H2*,CRPHOT=>H,H+,e-") +
+                CO_findrk("H2*,CRPHOT=>H,H") + CO_findrk("H2*,CRPHOT=>H+,H-") + CO_findrk("H2*,CRPHOT=>H+,H,e-");
+
+        if( findspecies("H-")->hevmol > 0. && hmi.rel_pop_LTE_Hmin > 0. )
+        {
+                hmi.hmidep = (double)findspecies("H-")->hevmol/ 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*findspecies("H-")->hevmol - hmi.HMinus_induc_rec_cooling;
+
+        /* H2+  +  HNU  =>  H+  + H */
+        h2phet = GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.); /* Now used entirely for side-effect.  Yuk */
+        hmi.h2plus_heat = thermal.HeatNet*findspecies("H2+")->hevmol;
+
+        /* 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 = findspecies("H2")->hevmol/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 = findspecies("H2+")->hevmol/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 = findspecies("H3+")->hevmol/hmi.rel_pop_LTE_H3p;
+        }
+        else
+        {
+                hmi.h3pdep = 1.;
+        }
+
+        free(c[0]);
+        free(c);
+        free(b);
+        return;
+}

branches/newmole/source/mole_newton_step.cpp

@@ -14 +14 @@
 #include "ionbal.h"
 #include "dense.h"
+#include "gammas.h"
+#include "opacity.h"
 #include "secondaries.h"
 #include "mole.h"
 #include "mole_priv.h"
-#include "opacity.h"
 #include "rfield.h"
 #include "thermal.h"
-#include "timesc.h"
 #include "trace.h"
 #include "phycon.h"
 #include "doppvel.h"
 #include "thirdparty.h"
-#include "gammas.h"
 #include "h2.h"
 #include "dynamics.h"
 #include "conv.h"
-#include "radius.h"
 #include "hextra.h" 
 #include "hmi.h"
 #include "taulines.h"
 #include "coolheavy.h"
-#include "mole_co_atom.h"
 
 /* HP cc cannot compile following except in -O1 mode */
@@ -631 +628 @@
 
         CO_return_cached_species();
-
-        hmi.H2_total = 0.;
-        /* this is where the transition struc expects to find the H2 abundance */
-        dense.xIonDense[LIMELM+2][0] = 0.;
-        hmi.hmihet = 0.;
-        hmi.h2plus_heat = 0.;
-        hmi.H2Opacity = 0.;
-        hmi.hmicol = 0.;
-        hmi.hmidep = 1.;
-
-        if (0) 
-        {
-                fprintf(stdout,"SILICON %ld snk %g %g\n",nzone,mole.sink[ipSILICON][0],mole.sink[ipSILICON][1]);
-                fprintf(stdout,"SILICON src %g %g\n",mole.source[ipSILICON][0],mole.source[ipSILICON][1]);
-
-                fprintf(stdout,"HYDROGEN snk %g %g\n",mole.sink[ipHYDROGEN][0],mole.sink[ipHYDROGEN][1]);
-                fprintf(stdout,"HYDROGEN src %g %g\n",mole.source[ipHYDROGEN][0],mole.source[ipHYDROGEN][1]);
-                fprintf(stdout,"CARBON snk %g %g\n",mole.sink[ipCARBON][0],mole.sink[ipCARBON][1]);
-                fprintf(stdout,"CARBON src %g %g\n",mole.source[ipCARBON][0],mole.source[ipCARBON][1]);
-                fprintf(stdout,"NEON snk %g %g\n",mole.sink[ipNEON][0],mole.sink[ipNEON][1]);
-                fprintf(stdout,"NEON src %g %g\n",mole.source[ipNEON][0],mole.source[ipNEON][1]);
-                fprintf(stdout,"NEON dense %g %g\n",findspecies("Ne")->hevmol,findspecies("Ne+")->hevmol);
-        }
         
         return;
 }
 
-void mole_effects()
-{
-        double 
-                h2phet,
-                plte;
-        realnum abundan;
-
-        double *b, **c;
-        long int i;
-
-        DEBUG_ENTRY( "mole_effects()" );
-
-        b = (double *)MALLOC((size_t)mole.num_total*sizeof(double));
-        c = (double **)MALLOC((size_t)mole.num_total*sizeof(double *));
-        c[0] = (double *)MALLOC((size_t)mole.num_total*
-                                                                                                        mole.num_total*sizeof(double));
-        for(i=1;i<mole.num_total;i++)
-        {
-                c[i] = c[i-1]+mole.num_total;
-        }
-
-        mole_eval_balance(mole.num_total,b,c);  
-
-        co.CODissHeat = (realnum)(CO_findrate("PHOTON,CO=>C,O")*1e-12);
-
-        thermal.heating[0][9] = co.CODissHeat;
-
-        /* the real multi-level model molecule */
-        abundan = (realnum) findspecies("CO")->hevmol;
-        /* IonStg and nelem were set to 2, 0 in makelevlines */
-        ASSERT( C12O16Rotate[0].Hi->IonStg < LIMELM+2 );
-        ASSERT( C12O16Rotate[0].Hi->nelem-1 < LIMELM+2 );
-        dense.xIonDense[ C12O16Rotate[0].Hi->nelem-1][C12O16Rotate[0].Hi->IonStg-1] = abundan;
-
-        CO_PopsEmisCool(&C12O16Rotate, nCORotate,abundan , "12CO",
-                &CoolHeavy.C12O16Rot,&CoolHeavy.dC12O16Rot );
-        /*if( nzone>400 )
-                fprintf(ioQQQ,"DEBUG co cool\t%.2f\t%.4e\t%li\t%.4e\t%.4e\n",
-                        fnzone,CoolHeavy.C12O16Rot,nCORotate,abundan,phycon.te );*/
-
-        abundan = (realnum) findspecies("CO")->hevmol/co.C12_C13_isotope_ratio;
-        /* IonStg and nelem were set to 3, 0 in makelevlines */
-        ASSERT( C13O16Rotate[0].Hi->IonStg < LIMELM+2 );
-        ASSERT( C13O16Rotate[0].Hi->nelem-1 < LIMELM+2 );
-        dense.xIonDense[ C13O16Rotate[0].Hi->nelem-1][C13O16Rotate[0].Hi->IonStg-1] = abundan;
-
-        CO_PopsEmisCool(&C13O16Rotate, nCORotate,abundan ,"13CO",
-                &CoolHeavy.C13O16Rot,&CoolHeavy.dC13O16Rot );
-
-        /* 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      */
-        {
-                int ipMH2, ipMH2s;
-                ipMH2 = findspecies("H2")->index;
-                ipMH2s = findspecies("H2*")->index;
-                hmi.H2_rate_destroy = 0;
-                if (ipMH2 != -1)
-                        hmi.H2_rate_destroy += -findspecies("H2")->hevmol * c[ipMH2][ipMH2];
-                if (ipMH2s != -1)
-                        hmi.H2_rate_destroy += -findspecies("H2*")->hevmol * c[ipMH2s][ipMH2s];
-                if (ipMH2 != -1 && ipMH2s != -1) 
-                        hmi.H2_rate_destroy += 
-                                - findspecies("H2")->hevmol * c[ipMH2][ipMH2s]
-                                - findspecies("H2*")->hevmol * c[ipMH2s][ipMH2];
-                hmi.H2_rate_destroy /= SDIV(hmi.H2_total);
-        }
-
-        /* establish local timescale for H2 molecule destruction */
-        if(findspecies("H2")->index != -1 && -c[findspecies("H2")->index][findspecies("H2")->index] > SMALLFLOAT )
-        {
-                /* units are now seconds */
-                timesc.time_H2_Dest_here = -1./c[findspecies("H2")->index][findspecies("H2")->index];
-        }
-        else
-        {
-                timesc.time_H2_Dest_here = 0.;