[FLASH-BUGS] Bug in iso7 reaction network

[Chris Malone]

I was recently trying to modify the iso7 reaction network when I discovered a 
bug in the calculation of the dense jacobian element pertaining to the triple 
alpha reaction.  When calculating the change in composition we have:

dydt(ihe4) = -3.0e0 * y(ihe4) * y(ihe4) * y(ihe4) * ratdum(ir3a) + ...

The -3.0e0 is from the fact that we are losing 3 helium nuclei per triple alpha 
reaction and the reaction rate has already been divided by 3! = 6 to prevent 
double counting.  Now to calculte the dense jacobian matrix element, dfdy(ihe4, 
ihe4), in diso7 (in FLASH version 2) we should simply take the above expression 
and differentiate with respect to helium composition giving:

dfdy(ihe4, ihe4) = -9.0e0 * y(ihe4) * y(ihe4) * ratdum(ir3a) + ...

However, in the FLASH implementation of iso7 the coefficient is -6.0e0 not 
-9.0e0.  I brought this to the attention of Frank Timmes and he said indeed it 
was a bug and that he was not suprised that the bug propagated through to FLASH 
because the paper he wrote that describes this reaction network also had the 
error.

Frank also commented that the bug does not affect the results of the 
calculation for the change in helium abundance but rather the efficiency of the 
algorithm.

Chris