NewMoleNickObservations

The following is a series of tests performed with the input script pdr_co_fully.in

r1207 seems to work fine, at least it stops for the intended reason. There are botches, but looking at the output the differences seem likely to be bugs fixed by Robin, and likely introduced by Nick.

r1208 does not work. The only change made, according to SVN, is that H- charge transfer was moved out of ion_recomb.cpp. Putting H- CT back into r1207 of ion_recomb.cpp has the expected outcome of getting pdr_co_fully.in to work. The comments state that this was originally removed because the code in ion_recomb was double-counting this process. Is it possible to leave this in ion_recomb and remove it from the molecular network?

r1209 does not work. When H- CT from ion_recomb.cpp r1207 is put into r1209, the script still does not work. Therefore, the changes made between r1208 and r1209 introduced new problems. The files that changed between r1208 and r1209 are hydrolevelpop.cpp and iso_ionize_recombine.cpp. The reasons for the changes are "Move molecular effects into level matrix solve, following pattern of ion_solver" The code were:

HYDROLEVELPOP.CPP

335                      Codereview(); /* Check signs */
336 	  	         if( conv.nTotalIoniz > 1 || iteration > 1 )
337 	  	         {
338 	  	 
339 	  	                 /* these are the external source and sink terms */
340 	  	                 /* source first */
341 	  	                 creation[ipH1s] += mole.source[nelem][i]/SDIV(dense.xIonDense[nelem][nelem+1]);
342 	  	 
343 	  	                 for( ipLo=ipH1s; ipLo<iso.numLevels_local[ipH_LIKE][nelem]; ++ipLo )
344 	  	                 {
345 	  	                         z[ipLo][ipLo] += mole.sink[nelem][ipH_LIKE];
346 	  	                 }
347 	  	         }

and, for ISO_IONIZE_RECOMBINE.CPP

{{{
/* The following terms in co.c destroy He+ and thereby need to be included * here. Also included in the sum of co.hep_destroy are the contributions * to the recombination due to hmole_step.*/

                 hep_destroy = (float) CO_sink_rate("He+"); 	 

                 iso.Ratecont2Level[ipISO][nelem][ipHe1s1S] += hep_destroy; 	 

                 /*fprintf(ioQQQ,"DEBUG co hep destroy %e\

",co.hep_destroy );*/ }}}

So the big changes that led to the model not working involve H- charge transfer and reorganizing how He interacts with the chemistry.