Changes related to

a) Incorrect convergence behaviour. I added a check for both the
global averaged ionized and neutral fractions where the code before
was only checking for the global ionized fraction. Once the volume
becomes mostly ionized, checking only for the global ionized fraction
gives wrong convergence results.

b) Corrections of several mistakes introduced in the August 2013
 - importing of the isothermal and S_star_nominal variables from
   c2ray_parameters
Added an MPI_BARRIER to deal with a potential problem on stampede.
This did not seem to solve it but does not hurt.
Added two MPI_BARRIER calls to make sure that the different MPI
processes remain more or less in sync around the broadcasting
of the initial conditions.

This change was triggered by an MPI problem on stampede. However,
it seems that adding these calls did not solve the problem. Still,
they do not hurt so we can leave them in.
Added the definition for the allocatable array phiheat as well as
its allocation statement. This array is meant to hold the heating
rate in case of a non-isothermal run. All statements involving
phiheat should be wrapped with
if (.not.isothermal)
statements.

The non-isothermal version of C2Ray_3Dm is still under development.
These changes are a step in the direction of making the code capable
of running both isothermal and non-isothermal.