Changeset 1848 for branches/newmole/source/mole_drive.cpp

Timestamp:

03/15/08 13:18:03 (10 months ago)

Author:

rjrw

Message:

Tidy & correct H2 rate diagnostics

Files:

• branches/newmole/source/mole_drive.cpp (modified) (7 diffs)

branches/newmole/source/mole_drive.cpp

@@ -235 +235 @@
         realnum abundan;
 
-        double *b, **c;
+        valarray<double> b(mole.num_total); 
+        double **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*
@@ -249 +249 @@
         }
 
-        mole_eval_balance(mole.num_total,b,c);  
+        mole_eval_balance(mole.num_total,&b[0],c);      
 
         for (long int nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
@@ -299 +299 @@
 
         /* save rate H2 is destroyed units s-1 */
+
+        /* total H2 - all forms */
+        hmi.H2_total = (realnum) findspecies("H2*")->den + findspecies("H2")->den;
+        /* first guess at ortho and para densities */
+        h2.ortho_density = 0.75*hmi.H2_total;
+        h2.para_density = 0.25*hmi.H2_total;
+
         /* >>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 */
@@ -320 +327 @@
 
         /* establish local timescale for H2 molecule destruction */
-        if(findspecies("H2")->index != -1 && -c[findspecies("H2")->index][findspecies("H2")->index] > SMALLFLOAT )
+        if(hmi.H2_rate_destroy > SMALLFLOAT )
         {
                 /* units are now seconds */
-                timesc.time_H2_Dest_here = -1./c[findspecies("H2")->index][findspecies("H2")->index];
+                timesc.time_H2_Dest_here = 1./hmi.H2_rate_destroy;
         }
         else
@@ -331 +338 @@
         
         /* local timescale for H2 formation 
-         * both grains and radiative attachment */
-        timesc.time_H2_Form_here = (mole_findrk("H,H,grnh=>H2*")+mole_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]);
+         * 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,grnh=>H2*")+mole_findrate("H,H,grnh=>H2"));
         /* 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;
+                timesc.time_H2_Form_here =  hmi.H2_total/timesc.time_H2_Form_here;
         }
         else
@@ -347 +352 @@
                 timesc.time_H2_Form_here = 0.;
         }
-
-
-        /* total H2 - all forms */
-        hmi.H2_total = (realnum) findspecies("H2*")->den + findspecies("H2")->den;
-        /* 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
@@ -708 +706 @@
         free(c[0]);
         free(c);
-        free(b);
 
         mole_h_rate_diagnostics();