/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*PrtComment analyze model, generating comments on its features */
/*badprt print out coolants if energy not conserved */
/*chkCaHeps check whether CaII K and H epsilon overlap */
/*outsum sum outward continuum beams */
/*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
#include "cddefines.h"
#include "physconst.h"
#include "doppvel.h"
#include "cddrive.h"
#include "lines_service.h"
#include "iso.h"
#include "continuum.h"
#include "hyperfine.h"
#include "dense.h"
#include "grainvar.h"
#include "version.h"
#include "rt.h"
#include "he.h"
#include "ionbal.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "lines.h"
#include "trace.h"
#include "map.h"
#include "hmi.h"
#include "punch.h"
#include "h2.h"
#include "converge.h"
#include "dynamics.h"
#include "opacity.h"
#include "geometry.h"
#include "elementnames.h"
#include "ca.h"
#include "broke.h"
#include "jeans_mass.h"
#include "pressure.h"
#include "co.h"
#include "atoms.h"
#include "abund.h"
#include "negdrg.h"
#include "rfield.h"
#include "colden.h"
#include "phycon.h"
#include "timesc.h"
#include "hextra.h"
#include "radius.h"
#include "iterations.h"
#include "fudgec.h"
#include "called.h"
#include "reason.h"
#include "magnetic.h"
#include "helike.h"
#include "wind.h"
#include "secondaries.h"
#include "struc.h"
#include "oxy.h"
#include "input.h"
#include "thermal.h"
#include "atmdat.h"
#include "warnings.h"
#include "agecheck.h"
#include "prt.h"

/*chkCaHeps check whether CaII K and H epsilon overlap */
static void chkCaHeps(double *totwid);

/*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
static void prt_smooth_predictions(void);

/*badprt print out coolants if energy not conserved */
static void badprt(double total);

/*outsum sum outward continuum beams */
static void outsum(double *outtot, 
  double *outin, 
  double *outout);

void PrtComment(void)
{
	char chLbl[11], 
	  chLine[INPUT_LINE_LENGTH];

	int lgAbort,
	  lgThick,
	  lgQHeat;

	long int _o, 
	  dum1,
	  dum2,
	  i, 
	  imas, 
	  ipLo ,
	  ipHi ,
	  ipISO,
	  nelem, 
	  isav, 
	  j, 
	  nc, 
	  nd, 
	  nline, 
	  nn, 
	  nneg, 
	  ns, 
	  nw;

	double big_ion_jump, 
	  absint ,
	  aj, 
	  alpha, 
	  big, 
	  bigm, 
	  sum_coolants, 
	  comfrc, 
	  flow_energy,
	  differ, 
	  error, 
	  flur, 
	  freqn, 
	  rate, 
	  ratio, 
	  sum_recom_lines, 
	  rel, 
	  small, 
	  tauneg, 
	  ts ,
	  HBeta, 
	  relfl ,
	  relhm,  
	  fedest,
	  GetHeat, 
	  SumNeg, 
	  c4363, 
	  t4363, 
	  outin, 
	  outout, 
	  totwid, 
	  outtot;

	double VolComputed , VolExpected , ConComputed , ConExpected;

#	ifdef DEBUG_FUN
	fputs( "<+>PrtComment()\n", debug_fp );
#	endif

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " PrtComment called.\n" );
	}

	/* 
	 * enter all comments cautions warnings and surprises into a large
	 * stack of statements
	 * at end of this routine is master call to printing routines
	 */
	iterations.lgIterAgain = FALSE;

	/* initialize the line saver */
	wcnint();

	if( version.nBetaVer > 0 )
	{
		sprintf( chLine, 
			"  !This is beta test version%4ld and is intended for testing only.", 
		  version.nBetaVer );
		bangin(chLine);
	}

	/* this flag set by call to fixit routine,
	 * to show that parts of the code need repair. */
	if( broke.lgFixit )
	{
		sprintf( chLine, "  !The code needs to be fixed - search for fixit()." );
		bangin(chLine);
	}

	/* this flag set by call to CodeReview routine,
	 * to show that parts of the code need to be reviewed. */
	if( broke.lgCheckit )
	{
		sprintf( chLine, "  !New code needs to be reviewed - search for CodeReview()." );
		bangin(chLine);
	}

	/* say why calculation stopped */
	if( conv.lgBadStop )
	{
		/* this case stop probably was not intended */
		sprintf( warnings.chRgcln[0], " W-Calculation stopped because %21.21s Iteration%3ld of %ld", 
		  reason.chReason, iteration, iterations.itermx + 1 );
		sprintf( chLine, " W-Calculation stopped because %21.21s", 
		  reason.chReason );
		warnin(chLine);
		sprintf( chLine, " W-This was not intended." );
		warnin(chLine);
	}
	else
	{
		/* for iterate to convergence, print reason why it was not converged on 3rd and higher iterations */
		if( (conv.lgAutoIt && iteration != iterations.itermx + 1) && 
		  iteration > 2 )
		{
			sprintf( warnings.chRgcln[0], 
				"   Calculation stopped because %21.21s Iteration %ld of %ld, not converged due to %s", 
			  reason.chReason, 
			  iteration, 
			  iterations.itermx + 1, 
			  conv.chNotConverged  );
		}
		else
		{
			sprintf( warnings.chRgcln[0], 
				"   Calculation stopped because %21.21s Iteration %ld of %ld", 
			  reason.chReason, iteration, iterations.itermx + 1 );
		}
	}

	/* check whether stopped because default number of zones hit,
	 * not intended?? */
	if( (!geometry.lgZoneSet) && geometry.lgZoneTrp )
	{
		conv.lgBadStop = TRUE;
		sprintf( chLine, 
			" W-Calculation stopped because default number of zones reached.  Was this intended???" );
		warnin(chLine);
		sprintf( chLine, 
			" W-Default limit can be increased while retaining this check with the SET NEND command." );
		warnin(chLine);
	}

	/* check whether stopped because zones too thin - only happens when glob set
	 * and depth + dr = depth
	 * not intended */
	if( radius.lgDrMinUsed || radius.lgdR2Small )
	{
		conv.lgBadStop = TRUE;
		sprintf( chLine, 
			" W-Calculation stopped zone thickness became too thin.   This was not intended." );
		warnin(chLine);
		sprintf( chLine, 
			" W-The most likely reason was an uncontrolled oscillation." );
		warnin(chLine);
		ShowMe();
	}
	if( radius.lgdR2Small )
	{
		sprintf( chLine, 
			" W-This happened because the globule scale became very small relative to the depth." );
		warnin(chLine);
		sprintf( chLine, 
			" W-This problem is described in Hazy." );
		warnin(chLine);
	}

	/* possible that last zone does not have stored temp - if so
	 * then make it up - this happens for some bad stops */
	ASSERT( nzone < struc.nzlim );

	if( struc.testr[nzone-1] == 0. )
		struc.testr[nzone-1] = struc.testr[nzone-2];

	if( struc.ednstr[nzone-1] == 0. )
		struc.ednstr[nzone-1] = struc.ednstr[nzone-2];

	/* give indication of geometry */
	rel = radius.depth/radius.rinner;
	if( rel < 0.1 )
	{
		sprintf( warnings.chRgcln[1], "   The geometry is plane-parallel." );
	}
	else if( rel >= 0.1 && rel < 3. )
	{
		sprintf( warnings.chRgcln[1], "   The geometry is a thick shell." );
	}
	else
	{
		sprintf( warnings.chRgcln[1], "   The geometry is spherical." );
	}

	/* levels of warnings: Warning   (possibly major problems)
	 *                     Caution   (not likely to invalidate the results)
	 *                     [      !] surprise, but not a problem
	 *                     [nothing] interesting note
	 */
	/* this flag set by call to routine broken ( ); 
	 * and show that the code is broken. */
	if( broke.lgBroke )
	{
		sprintf( chLine, " W-The code is broken - search for broken()." );
		warnin(chLine);
	}

	/* incorrect electron density detected in radinc during the calculation */
	if( dense.lgEdenBad )
	{
		if( dense.nzEdenBad == nzone )
		{
			sprintf( chLine, " C-The assumed electron density was incorrect for the last zone." );
			caunin(chLine);
			sprintf( chLine, " C-Did a temperature discontinuity occur??" );
			caunin(chLine);
		}
		else
		{
			sprintf( chLine, " W-The assumed electron density was incorrect during the calculation.  This is bad." );
			warnin(chLine);
			ShowMe();
			warnin(chLine);
		}
	}

	if( conv.lgAbort )
	{
		sprintf( chLine, " W-The calculation aborted.  Something REALLY went wrong!" );
		warnin(chLine);
	}

	/* thermal map was done but results were not ok */
	if( map.lgMapDone && !map.lgMapOK )
	{
		sprintf( chLine, "  !The thermal map had changes in slope - check map output." );
		bangin(chLine);
	}

	if( fudgec.nfudge > 0 )
	{
		sprintf( chLine, " !Fudge factors in place.  Why?" );
		bangin(chLine);
	}

	if( dense.gas_phase[ipHYDROGEN] > 1.1e13 )
	{
		if( dense.gas_phase[ipHYDROGEN] > 1e15 )
		{
			sprintf( chLine, " C-Density greater than 10**15, heavy elements are very uncertain." );
			caunin(chLine);
		}
		else
		{
			sprintf( chLine, " C-Density greater than 10**13" );
			caunin(chLine);
		}
	}

	/* HBeta is used later in the code to check on line intensites */
	if( cdLine("Pump",4861,&relfl,&absint)<=0 )
	{
		fprintf( ioQQQ, " Did not find flur H-beta\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	/* now find total Hbeta */
	if( cdLine( "totl",4861,&HBeta,&absint)<=0 )
	{
		fprintf( ioQQQ, " Did not find totl H-beta\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	/* check on continuum pumping of Balmer lines */
	if( HBeta > 0. )
	{
		flur = relfl/HBeta;
		if( flur > 0.1 )
		{
			sprintf( chLine, "  !Continuum fluorescent production of H-beta was very important." );
			bangin(chLine);
		}
		else if(flur > 0.01 )
		{
			sprintf( chLine, "   Continuum fluorescent production of H-beta was significant." );
			notein(chLine);
		}
	}

	/* check if there were problems with initial wind velocity */
	if( wind.windv0 > 0. && ((!wind.lgWindOK) || wind.windv < 1e6) )
	{
		sprintf( chLine, " C-Wind velocity below sonic point; solution is not valid." );
		caunin(chLine);
	}

	/* now confirm that mass flux here is correct */
	if( wind.windv0 != 0. )
	{
		rel = wind.emdot/(wind.windv*dense.gas_phase[ipHYDROGEN])/radius.r1r0sq;
		if( fabs(1.-rel)> 0.02 )
		{
			sprintf( chLine, " C-Wind mass flux error is %g%%",fabs(1.-rel)*100. );
			caunin(chLine);
			fprintf(ioQQQ,"DEBUG emdot\t%.3e\t%.3e\t%.3e\t%.3e\n",
				wind.emdot , wind.windv*dense.gas_phase[ipHYDROGEN],wind.windv,dense.gas_phase[ipHYDROGEN]);
		}
	}

	/* check that we didn't overrun zone scale */
	if( nzone >= struc.nzlim )
	{
		TotalInsanity();
	}

	/* check whether energy is conserved
	 * following is outward continuum */
	outsum(&outtot,&outin,&outout);
	/* sum_recom_lines is the sum of all recombination line energy */
	sum_recom_lines = totlin('r');
	if( sum_recom_lines == 0. )
	{
		sprintf( chLine, "  !Total recombination line energy is 0." );
		bangin(chLine);
	}

	/* sum_coolants is the sum of all coolants */
	sum_coolants = totlin('c');
	if( sum_coolants == 0. )
	{
		sprintf( chLine, "  !Total cooling is zero." );
		bangin(chLine);
	}

	if( wind.windv < 0. )
	{
		/* approximate energy density coming out in wind
		 * should ideally include flow into grid and internal energies */
		flow_energy = (2.5*struc.GasPressure[0]+0.5*struc.DenMass[0]*wind.windv0*wind.windv0)
			*(-wind.windv0);
	}
	else
	{
		flow_energy = 0.;
	}

	/* TotalLumin is total incident photon luminosity, evaluated in setup
	 * sum_coolants is evaluated here, is sum of all coolants */
	if( (sum_coolants + sum_recom_lines + flow_energy ) > continuum.TotalLumin && (!thermal.lgTSetOn) )
	{
		if( (hextra.cryden == 0.) && ((hextra.TurbHeat+DoppVel.DispScale) == 0.) && !secondaries.lgCSetOn )
		{

			sprintf( chLine, 
				" W-Radiated luminosity (cool+rec+wind=%.2e+%.2e+%.2e) is greater than that in incident radiation (TotalLumin=%8.2e).  Power radiated was %8.2e", 
			  sum_coolants, sum_recom_lines, flow_energy , continuum.TotalLumin, thermal.power );
			warnin(chLine);
			/* write same thing directly to output (above will be sorted later) */
			fprintf( ioQQQ, "\n\n This calculation DID NOT CONSERVE ENERGY!\n" );

			/* the case b command can cause this problem - say so if case b was set */
			if( opac.lgCaseB )
				fprintf( ioQQQ, "\n The CASE B command was entered - this can have unphysical effects, including non-conservation of energy.\n Why was it needed?\n\n" );

			/* print out significant heating and cooling */
			badprt(continuum.TotalLumin);

			sprintf( chLine, " W-Something is really wrong: the ratio of radiated to incident luminosity  is %10.2e.", 
			  (sum_coolants + sum_recom_lines)/continuum.TotalLumin );
			warnin(chLine);

			/* this can be caused by the force te command */
			if( thermal.ConstTemp > 0. )
			{
				fprintf( ioQQQ, " This may have been caused by the FORCE TE command.\n" );
				fprintf( ioQQQ, " Remove it and run again.\n" );
			}
			else
			{
				ShowMe();
			}
			warnin(chLine);
		}
	}

	/* warn if cosmic rays and bmag field both present */
	if( hextra.cryden*magnetic.lgB > 0. )
	{
		sprintf( chLine, 
			" W-Magnetic field and cosmic rays are both present.  This is not correctly treated." );
		warnin(chLine);
	}

	/* PrtComment if test code is in place */
	if( lgTestCodeCalled )
	{
		sprintf( chLine, "  !Test code is in place." );
		bangin(chLine);
	}

	/* lgComUndr set to .true. if compton cooling rate underflows to 0 */
	if( rfield.lgComUndr )
	{
		sprintf( chLine, 
			"  !Compton cooling rate underflows to zero.  Is this important?" );
		bangin(chLine);
	}

	/* make note if input stream contained an underscore, which was converted into a space */
	if( input.lgUnderscoreFound )
	{
		sprintf( chLine, 
			"  !Some input lines contained underscores, these were changed to spaces." );
		bangin(chLine);
	}

	/* make note if input stream contained a left or right bracket, which was converted into a space */
	if( input.lgBracketFound )
	{
		sprintf( chLine, 
			"  !Some input lines contained [ or ], these were changed to spaces." );
		bangin(chLine);
	}

	/* lgHionRad set to .true. if no hydrogen ionizing radiation */
	if( rfield.lgHionRad )
	{
		sprintf( chLine, 
			"  !There is no hydrogen-ionizing radiation.  Was this intended?" );
		bangin(chLine);
	}

	/* check whether certain zones were thermally unstable */
	if( thermal.nUnstable > 0 )
	{
		sprintf( chLine, 
			"   Derivative of net cooling negative and so possibly thermally unstable in%4ld zones.", 
		  thermal.nUnstable );
		notein(chLine);
	}

	/* generate a bang if a large fraction of the zones were unstable */
	if( nzone > 1 && 
		(float)(thermal.nUnstable)/(float)(nzone) > 0.25 )
	{
		sprintf( chLine, 
			"  !A large fraction of the zones were possibly thermally unstable,%4ld out of%4ld", 
		  thermal.nUnstable, nzone );
		bangin(chLine);
	}

	/* comment if negative coolants were ever significant */
	if( thermal.CoolHeatMax > 0.2 )
	{
		sprintf( chLine, 
			"  !Negative cooling reached %6.1f%% of the local heating, due to %4.4s %.1f", 
		  thermal.CoolHeatMax*100., thermal.chCoolHeatMax, thermal.wlCoolHeatMax );
		bangin(chLine);
	}
	else if( thermal.CoolHeatMax > 0.05 )
	{
		sprintf( chLine, 
			"   Negative cooling reached %6.1f%% of the local heating, due to %4.4s %.2f", 
		  thermal.CoolHeatMax*100., thermal.chCoolHeatMax, thermal.wlCoolHeatMax );
		notein(chLine);
	}

	/* check if advectionn heating was important */
	if( dynamics.HeatMax > 0.05 )
	{
		sprintf( chLine, 
			"  !Advection heating reached %.2f%% of the local heating.", 
		  dynamics.HeatMax*100. );
		bangin(chLine);
	}
	else if( dynamics.HeatMax > 0.005 )
	{
		sprintf( chLine, 
			"   Advection heating reached %.2f%% of the local heating.", 
		  dynamics.HeatMax*100. );
		notein(chLine);
	}

	/* check if advectionn cooling was important */
	if( dynamics.CoolMax > 0.05 )
	{
		sprintf( chLine, 
			"  !Advection cooling reached %.2f%% of the local cooling.", 
		  dynamics.CoolMax*100. );
		bangin(chLine);
	}
	else if( dynamics.CoolMax > 0.005 )
	{
		sprintf( chLine, 
			"   Advection cooling reached %.2f%% of the local heating.", 
		  dynamics.CoolMax*100. );
		notein(chLine);
	}

	/* check if thermal ionization of ground state of hydrogen was important */
	if( hydro.HCollIonMax > 0.10 )
	{
		sprintf( chLine, 
			"  !Thermal collisional ionization of H reached %.2f%% of the local ionization rate.", 
		  hydro.HCollIonMax*100. );
		bangin(chLine);
	}
	else if( hydro.HCollIonMax > 0.005 )
	{
		sprintf( chLine, 
			"   Thermal collisional ionization of H reached %.2f%% of the local ionization rate.", 
		  hydro.HCollIonMax*100. );
		notein(chLine);
	}

	/* check if lookup table for Hummer & Storey case B was exceeded */
	if( !atmdat.lgHCaseBOK[1][ipHYDROGEN]  )
	{
		sprintf( chLine, 
			"   Te-ne bounds of Case B lookup table exceeded, H I Case B line intensities set to zero." );
		notein(chLine);
	}
	if( !atmdat.lgHCaseBOK[1][ipHELIUM]  )
	{
		sprintf( chLine, 
			"   Te-ne bounds of Case B lookup table exceeded, He II Case B line intensities set to zero." );
		notein(chLine);
	}

	/* check if secondary ionization of hydrogen was important */
	if( secondaries.SecHIonMax > 0.10 )
	{
		sprintf( chLine, 
			"  !Suprathermal collisional ionization of H reached %.2f%% of the local H ionization rate.", 
		  secondaries.SecHIonMax*100. );
		bangin(chLine);
	}
	else if( secondaries.SecHIonMax > 0.005 )
	{
		sprintf( chLine, 
			"   Suprathermal collisional ionization of H reached %.2f%% of the local H ionization rate.", 
		  secondaries.SecHIonMax*100. );
		notein(chLine);
	}

	/* check if H2 vib-deexcitation heating was important */
	if( hmi.HeatH2DexcMax > 0.05 )
	{
		sprintf( chLine, 
			"  !H2 vib deexec heating reached %.2f%% of the local heating.", 
		  hmi.HeatH2DexcMax*100. );
		bangin(chLine);
	}
	else if( hmi.HeatH2DexcMax > 0.005 )
	{
		sprintf( chLine, 
			"   H2 vib deexec heating reached %.2f%% of the local heating.", 
		  hmi.HeatH2DexcMax*100. );
		notein(chLine);
	}

	/* check if H2 vib-deexcitation heating was important */
	if( hmi.CoolH2DexcMax > 0.05 )
	{
		sprintf( chLine, 
			"  !H2 deexec cooling reached %.2f%% of the local heating.", 
		  hmi.CoolH2DexcMax*100. );
		bangin(chLine);
	}
	else if( hmi.CoolH2DexcMax > 0.005 )
	{
		sprintf( chLine, 
			"   H2 deexec cooling reached %.2f%% of the local heating.", 
		  hmi.CoolH2DexcMax*100. );
		notein(chLine);
	}

	/* check if charge transfer ionization of hydrogen was important */
	if( atmdat.HIonFracMax > 0.10 )
	{
		sprintf( chLine, 
			"  !Charge transfer ionization of H reached %.1f%% of the local H ionization rate.", 
		  atmdat.HIonFracMax*100. );
		bangin(chLine);
	}
	else if( atmdat.HIonFracMax > 0.005 )
	{
		sprintf( chLine, 
			"   Charge transfer ionization of H reached %.2f%% of the local H ionization rate.", 
		  atmdat.HIonFracMax*100. );
		notein(chLine);
	}

	/* check if charge transfer heating cooling was important */
	if( atmdat.HCharHeatMax > 0.05 )
	{
		sprintf( chLine, 
			"  !Charge transfer heating reached %.2f%% of the local heating.", 
		  atmdat.HCharHeatMax*100. );
		bangin(chLine);
	}
	else if( atmdat.HCharHeatMax > 0.005 )
	{
		sprintf( chLine, 
			"   Charge transfer heating reached %.2f%% of the local heating.", 
		  atmdat.HCharHeatMax*100. );
		notein(chLine);
	}

	if( atmdat.HCharCoolMax > 0.05 )
	{
		sprintf( chLine, 
			"  !Charge transfer cooling reached %.2f%% of the local heating.", 
		  atmdat.HCharCoolMax*100. );
		bangin(chLine);
	}
	else if( atmdat.HCharCoolMax > 0.005 )
	{
		sprintf( chLine, 
			"   Charge transfer cooling reached %.2f%% of the local heating.", 
		  atmdat.HCharCoolMax*100. );
		notein(chLine);
	}

	/* did Lya mase? */
	if( hydro.xLaMase < 0. )
	{
		if( hydro.xLaMase < -0.01 )
		{
			sprintf( chLine, "  !H Lya mased!!  Smallest optical depth was %.2e", 
			  hydro.xLaMase );
			bangin(chLine);
		}
		else
		{
			sprintf( chLine, "   H Lya mased.  Smallest optical depth was %.2e", 
			  hydro.xLaMase );
			notein(chLine);
		}
	}

	/* check whether photo from up level of Mg2 2798 ever important */
	if( atoms.xMg2Max > 0.1 )
	{
		sprintf( chLine, 
			"  !Photoionization of upper level of Mg II 2798 reached %.1f%% of the total Mg+ photo rate.", 
		  atoms.xMg2Max*100. );
		bangin(chLine);
	}
	else if( atoms.xMg2Max > 0.01 )
	{
		sprintf( chLine, 
			"   Photoionization of upper level of Mg II 2798 reached %.1f%% of the total Mg+ photo rate.", 
		  atoms.xMg2Max*100. );
		notein(chLine);
	}

	/* check whether photo from up level of [O I] 6300 ever important */
	if( oxy.poimax > 0.1 )
	{
		sprintf( chLine, 
			"  !Photoionization of upper levels of [O I] reached %.1f%% of the total O destruction rate.", 
		  oxy.poimax*100. );
		bangin(chLine);
	}
	else if( oxy.poimax > 0.01 )
	{
		sprintf( chLine, 
			"   Photoionization of upper levels of [O I] reached %.1f%% of the total O destruction rate.", 
		  oxy.poimax*100. );
		notein(chLine);
	}

	/* check whether photo from up level of [O III] 5007 ever important */
	if( (oxy.poiii2Max + oxy.poiii3Max) > 0.1 )
	{
		sprintf( chLine, 
			"  !Photoionization of upper levels of [O III] reached %.1f%% of the total O++ photo rate.", 
		  (oxy.poiii2Max + oxy.poiii3Max)*100. );
		bangin(chLine);
	}
	else if( (oxy.poiii2Max + oxy.poiii3Max) > 0.01 )
	{
		sprintf( chLine, 
			"   Photoionization of upper levels of [O III] reached %.1f%% of the total O++ photo rate.", 
		  (oxy.poiii2Max + oxy.poiii3Max)*100. );
		notein(chLine);
	}

	/* check whether photoionization of He 2trip S was important */
	if( he.frac_he0dest_23S > 0.1 )
	{
		sprintf( chLine, 
			"  !Destruction of He 2TriS reached %.1f%% of the total He0 dest rate"
			" at zone %li, %.1f%% of that was photoionization.", 
		  he.frac_he0dest_23S*100, 
		  he.nzone, 
		  he.frac_he0dest_23S_photo*100.  );
		bangin(chLine);
	}
	else if( he.frac_he0dest_23S > 0.01 )
	{
		sprintf( chLine, 
			"   Destruction of He 2TriS reached %.1f%% of the total He0 dest rate"
			" at zone %li, %.1f%% of that was photoionization.", 
		  he.frac_he0dest_23S*100, 
		  he.nzone, 
		  he.frac_he0dest_23S_photo*100.  );
		notein(chLine);
	}

	if( !helike.lgCritDensLMix && dense.lgElmtOn[ipHELIUM] )
	{
		sprintf( chLine,
			"   Critical density for l-mixing of He I not reached. "
			" More resolved levels are needed for accurate He I line ratios.");
		notein(chLine);
	}


	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
	{
		if( helike.lgFullSize[nelem] && dense.lgElmtOn[nelem] )
		{
			sprintf( chLine, "   Helike: Z=%li is a complete model due to continuum lowering.  Max n is %li\n",
				nelem+1,
				helike.n_HighestResolved[nelem]+helike.nCollapsed[nelem]);
			notein(chLine);
		}
	}

	/* frequency array may not have been defined for all energies */
	if( !rfield.lgMMok )
	{
		sprintf( chLine, 
			" C-Continuum not defined in extreme infrared - Compton scat, grain heating, not treated properly?" );
		caunin(chLine);
	}

	if( !rfield.lgHPhtOK )
	{
		sprintf( chLine, 
			" C-Continuum not defined at photon energies which ionize excited states of H, important for H- and ff heating." );
		caunin(chLine);
	}

	if( !rfield.lgXRayOK )
	{
		sprintf( chLine, 
			" C-Continuum not defined at X-Ray energies - Compton scattering and Auger ionization wrong?" );
		caunin(chLine);
	}

	if( !rfield.lgGamrOK )
	{
		sprintf( chLine, 
			" C-Continuum not defined at gamma-ray energies - pair production and Compton scattering OK?" );
		caunin(chLine);
	}

	if( continuum.lgCon0 )
	{
		sprintf( chLine, " C-Continuum zero at some energies." );
		caunin(chLine);
	}

	if( continuum.lgCoStarInterpolationCaution )
	{
		sprintf( chLine , " C-CoStarInterpolate interpolated between non-adjoining tracks, this may not be valid." );
		caunin(chLine);
	}

	if( rfield.lgOcc1Hi )
	{
		sprintf( chLine, 
			"  !The continuum occupation number at 1 Ryd is greater than unity." );
		bangin(chLine);
	}

	/* this flag set true it set dr forced first zone to be too big */
	if( radius.lgDR2Big )
	{
		sprintf( chLine, 
			" C-The thickness of the first zone was set larger than optimal by a SET DR command." );
		caunin(chLine);
		sprintf( chLine, 
			" C-Consider using the STOP THICKNESS command instead." );
		caunin(chLine);
	}

	/* check whether non-col excit of 4363 was important */
	if( cdLine("TOTL",4363,&t4363,&absint)<=0 )
	{
		fprintf( ioQQQ, " PrtComment could not find total O III 4363 with cdLine.\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	if( cdLine("Coll",4363,&c4363,&absint)<=0 )
	{
		fprintf( ioQQQ, " PrtComment could not find collis O III 4363 with cdLine.\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	/* only print this comment if 4363 is significant and density low */
	if( HBeta > 1e-35 )
	{
		/* >>chng 02 feb 27, lower ratio from -4 to -5, and raise density from 7 to 8 */
		if( t4363/HBeta > 1e-5 && dense.gas_phase[ipHYDROGEN] < 1e8 )
		{
			ratio = (t4363 - c4363)/t4363;
			if( ratio > 0.01 )
			{
				sprintf( chLine, 
					"  !Non-collisional excitation of [O III] 4363 reached %6.2f%% of the total.", 
				  ratio*100. );
				bangin(chLine);
			}
			else if( ratio > 0.001 )
			{
				sprintf( chLine, 
					"   Non-collisional excitation of [O III] 4363 reached %6.2f%% of the total.", 
				  ratio*100. );
				notein(chLine);
			}
		}
	}

	/* check for plasma shielding */
	if( rfield.lgPlasNu )
	{
		sprintf( chLine, 
			"  !The largest plasma frequency was %.2e Ryd = %.2e micron  The continuum is set to 0 below this.", 
		  rfield.plsfrqmax,
		  /* wavelength in microns */
		  RYDLAM/rfield.plsfrqmax/1e4);
		bangin(chLine);
	}

	if( rfield.occmax > 0.1 )
	{
		if( rfield.occmnu > 1e-4 )
		{
			sprintf( chLine, 
				"  !The largest continuum occupation number was %10.3e at %10.3e Ryd.", 
			  rfield.occmax, rfield.occmnu );
			bangin(chLine);
		}
		else
		{
			/* not surprising if occ num bigger than 1 around 1e-5 Ryd,
			 * since this is the case for 3K background */
			sprintf( chLine, 
				"   The largest continuum occupation number was %10.3e at%10.3e Ryd.", 
			  rfield.occmax, rfield.occmnu );
			notein(chLine);
		}
	}

	if( rfield.occmax > 1e4 && rfield.occ1nu > 0. )
	{
		/* occ1nu is energy (ryd) where continuum occ num falls below 1 */
		if( rfield.occ1nu < 0.0912 )
		{
			sprintf( chLine, 
				"   The continuum occupation number fell below 1 at %10.3e microns.", 
			  0.0912/rfield.occ1nu );
			notein(chLine);
		}
		else if( rfield.occ1nu < 1. )
		{
			sprintf( chLine, 
				"   The continuum occupation number fell  below 1 at %10.3e Angstroms.", 
			  912./rfield.occ1nu );
			notein(chLine);
		}
		else
		{
			sprintf( chLine, 
				"   The continuum occupation number fell  below 1 at %10.3e Ryd.", 
			  rfield.occ1nu );
			notein(chLine);
		}
	}

	if( rfield.tbrmax > 1e3 )
	{
		sprintf( chLine, 
			"  !The largest continuum brightness temperature was %10.3eK at %10.3e Ryd.", 
		  rfield.tbrmax, rfield.tbrmnu );
		bangin(chLine);
	}

	if( rfield.tbrmax > 1e4 )
	{
		/* tbr4nu is energy (ryd) where continuum bright temp falls < 1e4 */
		if( rfield.tbr4nu < 0.0912 )
		{
			sprintf( chLine, 
				"   The continuum brightness temperature fell below 10,000K at %10.3e microns.", 
			  0.0912/rfield.tbr4nu );
			notein(chLine);
		}
		else if( rfield.tbr4nu < 1. )
		{
			sprintf( chLine, 
				"   The continuum brightness temperature fell below 10,000K at %10.3e Angstroms.", 
			  912./rfield.tbr4nu );
			notein(chLine);
		}
		else
		{
			sprintf( chLine, 
				"   The continuum brightness temperature fell below 10,000K at %10.3e Ryd.", 
			  rfield.tbr4nu );
			notein(chLine);
		}
	}

	/* turbulence AND constant pressure do not make sense */
	if( DoppVel.TurbVel > 0. && strcmp(dense.chDenseLaw,"CPRE") == 0 )
	{
		sprintf( chLine, 
			"  !Both constant pressure and turbulence makes no physical sense???" );
		bangin(chLine);
	}

	/* filling factor AND constant pressure do not make sense */
	if( geometry.FillFac < 1. && strcmp(dense.chDenseLaw,"CPRE") == 0 )
	{
		sprintf( chLine, 
			"  !Both constant pressure and a filling factor makes no physical sense???" );
		bangin(chLine);
	}

	/* grains and solar abundances do not make sense */
	if( gv.lgDustOn && abund.lgAbnSolar )
	{
		sprintf( chLine, 
			"  !Grains are present, but the gas phase abundances were left at the solar default.  This is not physical." );
		bangin(chLine);
	}

	/* check if depletion command set but no grains, another silly thing to do */
	if( abund.lgDepln && !gv.lgDustOn )
	{
		sprintf( chLine, 
			"  !Grains are not present, but the gas phase abundances were depleted.  This is not physical." );
		bangin(chLine);
	}

	if( gv.lgDustOn )
	{
		long nBin=0L,nFail=0L;
		for( nd=0; nd < gv.nBin; nd++ )
		{
			if( gv.bin[nd]->QHeatFailures > 0L )
			{
				++nBin;
				nFail += gv.bin[nd]->QHeatFailures;
			}
		}
		if( nFail > 0 )
		{
			sprintf( chLine,
				 "  !The grain quantum heating treatment failed to converge %ld time(s) in %ld bin(s).", nFail, nBin );
			bangin(chLine);
		}
	}

	/* check if PAHs were present in the ionized region, with constant abundances */
	if( gv.lgDustOn )
	{
		int lgPAHsPresent_and_constant = FALSE;
		for( nd=0; nd < gv.nBin; nd++ )
		{
			lgPAHsPresent_and_constant = lgPAHsPresent_and_constant || 
				/* it is ok to have PAHs in the ionized region if the abundances vary */
				(gv.bin[nd]->lgPAHsInIonizedRegion && !gv.bin[nd]->lgDustVary);
		}
		if( lgPAHsPresent_and_constant )
		{
			sprintf( chLine,
				 " C-PAH's were present in the ionized region, this has never been observed in H II regions." );
			caunin(chLine);
		}
	}

	/* constant temperature greater than continuum energy density temperature */
	if( thermal.lgTSetOn && thermal.ConstTemp*1.0001 < phycon.TEnerDen )
	{
		sprintf( chLine, 
			" C-The continuum energy density temperature (%g K)"
			" is greater than the electron temperature (%g K).",
			phycon.TEnerDen , thermal.ConstTemp);
		caunin(chLine);
		sprintf( chLine, " C-This is unphysical." );
		caunin(chLine);
	}

	/* remark that grains not present but energy density was low */
	if( !gv.lgDustOn && phycon.TEnerDen < 800. )
	{
		sprintf( chLine, 
			"   Grains were not present but might survive in this environment (energy density temperature was %.2eK)", 
		  phycon.TEnerDen );
		notein(chLine);
	}

	/* call routine that will check age of cloud */
	AgeCheck();

	/* check on CaH and H-epsilon overlapping
	 * need to do this since need to refer to lines arrays */
	chkCaHeps(&totwid);
	if( totwid > 121. )
	{
		sprintf( chLine, "   H-eps and CaH overlap." );
		notein(chLine);
	}

	/* warning that something was turned off */
	if( !phycon.lgPhysOK )
	{
		sprintf( chLine, "  !A physical process has been disabled." );
		bangin(chLine);
	}

	/* check on lifetimes of [O III] against photoionization, only for low den */
	if( dense.gas_phase[ipHYDROGEN] < 1e8 )
	{
		if( oxy.r5007Max > 0.0263f )
		{
			sprintf( chLine, 
				"  !Photoionization of upper level of [O III] 5007 reached %10.2e%% of the radiative lifetime.", 
			  oxy.r5007Max*100. );
			bangin(chLine);
		}
		else if( oxy.r5007Max > 0.0263f/10.f )
		{
			sprintf( chLine, 
				"   Photoionization of upper level of [O III] 5007 reached %10.2e%% of the radiative lifetime.", 
			  oxy.r5007Max*100. );
			notein(chLine);
		}
		if( oxy.r4363Max > 1.78f )
		{
			sprintf( chLine, 
				"  !Photoionization of upper level of [O III] 4363 reached %10.2e%% of the radiative lifetime.", 
			  oxy.r4363Max*100. );
			bangin(chLine);
		}
		else if( oxy.r4363Max > 1.78f/10.f )
		{
			sprintf( chLine, 
				"   Photoionization of upper level of [O III] 4363 reached %10.2e%% of the radiative lifetime.", 
			  oxy.r4363Max*100. );
			notein(chLine);
		}
	}

	/* check whether total heating and cooling matched
	 * >>chng 97 mar 28, added GrossHeat, heat in terms normally heat-cool */
	error = fabs(thermal.power-thermal.totcol)/((thermal.power + thermal.totcol)/2.);
	if( thermal.lgTSetOn )
	{
		if( error > 0.05 )
		{
			sprintf( chLine, 
				"  !Heating - cooling mismatch =%5.1f%%. Caused by constant temperature assumption. ", 
			  error*100. );
			bangin(chLine);
		}
	}

	else
	{
		if( error > 0.05 && error < 0.2 )
		{
			sprintf( chLine, " C-Heating - cooling mismatch =%.1f%%. What\'s wrong?", 
			  error*100. );
			caunin(chLine);
		}
		else if( error >= 0.2 )
		{
			sprintf( chLine, " W-Heating - cooling mismatch =%.2e%%. What\'s wrong????", 
			  error*100. );
			warnin(chLine);
		}
	}

	/* say if Ly-alpha photo of Ca+ excited levels was important */
	if( ca.Ca2RmLya > 0.01 )
	{
		sprintf( chLine, 
			"   Photoionization of Ca+ 2D level by Ly-alpha reached %6.1f%% of the total rate out.", 
		  ca.Ca2RmLya*100. );
		notein(chLine);
	}

	/* check if Lya alpha ever hotter than gas */
	if( hydro.nLyaHot > 0 )
	{
		if( hydro.TLyaMax/hydro.TeLyaMax > 1.05 )
		{
			sprintf( chLine, 
				"  !The excitation temp of Lya exceeded the electron temp, largest value was%10.2eK (gas temp there was%10.2eK, zone%4ld)", 
			  hydro.TLyaMax, hydro.TeLyaMax, hydro.nZTLaMax );
			bangin(chLine);
		}
	}

	/* check if line absorption heating was important */

	/* get all negative lines, check if line absorption significant heat source
	 * this is used in "final" for energy budget print out */
	if( cdLine("Line",0,&SumNeg,&absint)<=0 )
	{
		fprintf( ioQQQ, " did not get sumneg cdLine\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	/* this is total heating */
	if( cdLine("TotH",0,&GetHeat,&absint)<=0 )
	{
		fprintf( ioQQQ, " did not get GetHeat cdLine\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	if( GetHeat > 0. )
	{
		SumNeg /= GetHeat;
		if( SumNeg > 0.1 )
		{
			sprintf( chLine, 
				"  !Line absorption heating reached %.2f%% of the global heating.", 
			  SumNeg*100. );
			bangin(chLine);
		}
		else if( SumNeg > 0.01 )
		{
			sprintf( chLine, 
				"   Line absorption heating reached %.2f%% of the global heating.", 
			  SumNeg*100. );
			notein(chLine);
		}
	}

	/* this is check of extra lines added with g-bar */
	if( input.lgSetNoBuffering )
	{
		sprintf( chLine, 
			"  !NO BUFFERING command was entered - this increases exec time by LARGE amounts.");
		bangin(chLine);
	}

	/* this is check of extra lines added with g-bar */
	if( thermal.GBarMax > 0.1 )
	{
		ASSERT( thermal.ipMaxExtra > 0 );
		strcpy( chLbl, chLineLbl(&TauLine2[thermal.ipMaxExtra-1]) );

		sprintf( chLine, 
			"  !G-bar cooling lines reached %.2f%% of the local cooling.  Line=%10.10s", 
		  thermal.GBarMax*100., chLbl );
		bangin(chLine);
	}

	else if( thermal.GBarMax > 0.01 )
	{
		strcpy( chLbl, chLineLbl(&TauLine2[thermal.ipMaxExtra-1]) );

		sprintf( chLine, 
			"   G-bar cooling lines reached %.2f%% of the local cooling.  Line=%10.10s", 
		  thermal.GBarMax*100., chLbl );
		notein(chLine);
	}

	/* this is check of hyperfine structure lines*/
	if( hyperfine.cooling_max > 0.1 )
	{
		sprintf( chLine, 
			"  !Hyperfine structure line cooling reached %.2f%% of the local cooling.", 
		  hyperfine.cooling_max*100.);
		bangin(chLine);
	}

	else if( hyperfine.cooling_max > 0.01 )
	{
		sprintf( chLine, 
			"   Hyperfine structure line cooling reached %.2f%% of the local cooling.", 
		  hyperfine.cooling_max*100. );
		notein(chLine);
	}

	/* line absorption heating reached more than 10% of local heating?
	 * HeatLineMax is largest heating(1,23)/htot */
	if( thermal.HeatLineMax > 0.1 )
	{
		if( thermal.levlmax == 1 )
		{
			/* main block of lines */
			/* >>chng 01 may 05, removed chGetLbl routine, which was here,
			 * replaced with chLineLbl routine and address of TauLines 
			 * should be no change in functionality */
			strcpy( chLbl, chLineLbl(&TauLines[thermal.ipHeatlmax] ) );
		}
		else if( thermal.levlmax == 2 )
		{
			/* level 2 lines */
			strcpy( chLbl, chLineLbl(&TauLine2[thermal.ipHeatlmax]) );
		}
		else if( thermal.levlmax == 3 )
		{
			/* C12O16 lines */
			strcpy( chLbl, chLineLbl(&C12O16Rotate[thermal.ipHeatlmax]) );
		}
		else if( thermal.levlmax == 4 )
		{
			/* C13O16 lines */
			strcpy( chLbl, chLineLbl(&C13O16Rotate[thermal.ipHeatlmax]) );
		}
		else
		{
			fprintf( ioQQQ, " PrtComment has insane levlmax,=%5ld\n", 
			  thermal.levlmax );
		}
		sprintf( chLine, 
			"  !Line absorption heating reached %.2f%% of the local heating - largest by level%2ld line %10.10s", 
		  thermal.HeatLineMax*100., thermal.levlmax, chLbl );
		bangin(chLine);
	}

	else if( thermal.HeatLineMax > 0.01 )
	{
		sprintf( chLine, 
			"   Line absorption heating reached %.2f%% of the local heating.", 
		  thermal.HeatLineMax*100. );
		notein(chLine);
	}

	if( ionbal.CompHeating_Max > 0.05 )
	{
		sprintf( chLine, 
			"  !Bound Compton heating reached %.2f%% of the local heating.", 
		  ionbal.CompHeating_Max*100. );
		bangin(chLine);
	}
	else if( ionbal.CompHeating_Max > 0.01 )
	{
		sprintf( chLine, 
			"   Bound Compton heating reached %.2f%% of the local heating.", 
		  ionbal.CompHeating_Max*100. );
		notein(chLine);
	}

	/* did Ly-alpha mase */
	if( hydro.lgHLyaMased )
	{
		sprintf( chLine, "   Ly-alpha mased." );
		notein(chLine);
	}

	if( dense.gas_phase[ipHYDROGEN] < 1e7 )
	{
		/* check on IR fine structure lines */
		lgThick = FALSE;
		tauneg = 0.;
		alpha = 0.;
		/* loop from 3, since 0 is dummy, 1 and 2 are spin-flip transitions of H and He */
		for( i=3; i <= nLevel1; i++ )
		{
			/* define IR as anything longward of 1 micron */
			if( TauLines[i].EnergyWN < 10000. )
			{
				if( TauLines[i].TauIn > 1. )
				{
					lgThick = TRUE;
					alpha = MAX2(alpha,(double)TauLines[i].TauIn);
				}
				else if( TauLines[i].TauIn < (float)tauneg )
				{
					tauneg = TauLines[i].TauIn;
					strcpy( chLbl, chLineLbl(&TauLines[i]) );
				}
			}
		}
		/* now print results, were any fine structure lines optically thick? */
		if( lgThick )
		{
			sprintf( chLine, 
				"  !Some infrared fine structure lines are optically thick:  largest tau was%10.2e", 
			  alpha );
			bangin(chLine);
		}
		/* did any fine structure lines mase? */
		if( tauneg < -0.01 )
		{
			sprintf( chLine, 
				"  !Some fine structure lines mased: line %s had optical depth %.2e", 
			  chLbl, tauneg );
			bangin(chLine);
		}
	}

	/* were any other lines masing? */
	tauneg = 0.;
	alpha = 0.;
	for( i=1; i <= nLevel1; i++ )
	{
		/* define UV as anything shortward of 1 micron */
		if( TauLines[i].EnergyWN >= 10000. )
		{
			if( TauLines[i].TauIn < (float)tauneg )
			{
				tauneg = TauLines[i].TauIn;
				strcpy( chLbl, chLineLbl(&TauLines[i]) );
			}
		}
	}

	/* did any lines mase? */
	if( tauneg < -0.01 )
	{
		sprintf( chLine, 
			"  !Some lines mased: most negative ion and tau were:%10.10s%10.2e", 
		  chLbl, tauneg );
		bangin(chLine);
	}

	/* this is check that at least a second iteration was done with sphere static,
	 * the test is orverridden with the (OK) option on the sphere static command,
	 * which sets geometry.lgStaticNoIt TRUE */
	if( geometry.lgStatic && iterations.lgLastIt && (iteration == 1) && 
		!geometry.lgStaticNoIt)
	{
		sprintf( chLine, " C-I must iterate when SPHERE STATIC is set." );
		caunin(chLine);
		iterations.lgIterAgain = TRUE;
	}

	/* caution if continuum is punched but only one iteration performed */
	if( punch.lgPunContinuum && iteration == 1 && iterations.lgLastIt)
	{
		sprintf( chLine, " C-I must iterate when punch continuum output is done." );
		caunin(chLine);
		iterations.lgIterAgain = TRUE;
	}

	/* say that we may be wasting time is punch continuum false but quantum 
	 * heating is turned on, first var set true when punch continuum is
	 * set, and second is set with "NO QHEAT" option on grains */
	lgQHeat = FALSE;
	for( nd=0; nd < gv.nBin; nd++ )
	{
		if( gv.bin[nd]->lgQHeat )
			lgQHeat = TRUE;
	}
	if( !punch.lgPunContinuum &&  lgQHeat )
	{
		sprintf( chLine, "  !Grain quantum heating is enabled but continuum or spectrum is not punched." );
		caunin(chLine);
		sprintf( chLine, "  !Consider turning off quantum heating with NO QHEAT option on grains command to save time." );
		caunin(chLine);
	}

	/*TODO	2	extend to all iso and elem */
	/* how important was induced two photon?? */
	if( iso.TwoNu_induc_dn_max[ipH_LIKE][ipHYDROGEN] > 1. )
	{
		sprintf( chLine, "  !Rate of induced H 2-photon emission reached %.2e s^-1", 
		  iso.TwoNu_induc_dn_max[ipH_LIKE][ipHYDROGEN] );
		bangin(chLine);
	}

	else if( iso.TwoNu_induc_dn_max[ipH_LIKE][ipHYDROGEN] > 0.01 )
	{
		sprintf( chLine, "   Rate of induced H 2-photon emission reached %.2e s^-1", 
		  iso.TwoNu_induc_dn_max[ipH_LIKE][ipHYDROGEN] );
		notein(chLine);
	}


	/* how important was induced recombination? */
	if( hydro.FracInd > 0.01 )
	{
		sprintf( chLine, 
			"   Induced recombination was %5.1f%% of the total for H level%3ld", 
		  hydro.FracInd*100., hydro.ndclev );
		notein(chLine);
	}

	if( hydro.fbul > 0.01 )
	{
		sprintf( chLine, 
			"   Stimulated emission was%6.1f%% of the total for H transition%3ld -%3ld", 
		  hydro.fbul*100., hydro.nbul + 1, hydro.nbul );
		notein(chLine);
	}

	/* check whether Fe II destruction of La was important */
	if( cdLine("Fe 2",1216,&fedest,&absint)<=0 )
	{
		fprintf( ioQQQ, " Did not find Fe II Lya\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	/* find total Lya for comparison */
	if( cdLine("TOTL",1216,&relhm,&absint)<=0 )
	{
		fprintf( ioQQQ, " Did not find Lya\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	if( relhm > 0. )
	{
		ratio = fedest/(fedest + relhm);
		if( ratio > 0.1 )
		{
			sprintf( chLine, "  !Fe II destruction of Ly-a removed %.1f%% of the line.", 
			  ratio *100.);
			bangin(chLine);
		}
		else if( ratio > 0.01 )
		{
			sprintf( chLine, "   Fe II destruction of Ly-a removed %.1f%% of the line.", 
			  ratio );
			notein(chLine);
		}
	}

	if( cdLine("H-CT",6563,&relhm,&absint)<=0 )
	{
		fprintf( ioQQQ, " Comment did not find H-CT H-alpha\n" );
		puts( "[Stop in PrtComment]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	if( HBeta > 0. )
	{
		if( relhm/HBeta > 0.01 )
		{
			sprintf( chLine, 
				"  !Mutual neutralization production of H-alpha was significant." );
			bangin(chLine);
		}
	}

	/* note about very high population in H n=2 rel to ground, set in hydrogenic */
	if( hydro.lgHiPop2 )
	{
		sprintf( chLine, 
			"   The population of H n=2 reached%10.2e relative to the ground state.", 
		  hydro.pop2mx );
		notein(chLine);
	}

	/* check where diffuse emission error */
	for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
	{
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			if( iso.CaseBCheck[ipISO][nelem] > 1.5 )
			{
				sprintf( chLine, 
					"   Ratio of computed diffuse emission to case B reached %g for iso %li element %li",
					iso.CaseBCheck[ipISO][nelem] , ipISO , nelem+1 );
				notein(chLine);
			}
		}
	}

	if( thermal.thist > 1e9 )
	{
		if( thermal.thist > 5e9 )
		{
			sprintf( chLine, " W-Electrons were relativistic; High TE=%10.2e", 
			  thermal.thist );
			warnin(chLine);
		}
		else
		{
			sprintf( chLine, " C-Electrons were mildly relativistic; High TE=%10.2e", 
			  thermal.thist );
			caunin(chLine);
		}
	}

	rate = timesc.TimeErode*2e-26;
	if( rate > 1e-35 )
	{
		/*  2E-26 is roughly cross section for photoerosion
		 *  see 
		 * >>refer	all	photoerode	Boyd and Ferland ApJ Let 318, L21. */
		ts = (1./rate)/3e7;
		if( ts < 1e3 )
		{
			sprintf( chLine, "  !Timescale-photoerosion of Fe=%.2e yr", 
			  ts );
			bangin(chLine);
		}
		else if( ts < 1e9 )
		{
			sprintf( chLine, "   Timescale-photoerosion of Fe=%.2e yr", 
			  ts );
			notein(chLine);
		}
	}

	/* check whether Compton heating was significant */
	comfrc = rfield.comtot/thermal.power;
	if( comfrc > 0.01 )
	{
		sprintf( chLine, "   Compton heating was %5.1f%% of the total.", 
		  comfrc*100. );
		notein(chLine);
	}

	/* check on relative importance of induced compton heating */
	if( comfrc > 0.01 && rfield.cinrat > 0.05 )
	{
		sprintf( chLine, 
			"  !Induced Compton heating was %10.2e of the total Compton heating.", 
		  rfield.cinrat );
		bangin(chLine);
	}

	/* check whether equilibrium timescales are short rel to Hubble time */
	if( timesc.tcmptn > 5e17 )
	{
		if( comfrc > 0.05 )
		{
			sprintf( chLine, 
				" C-Compton cooling is significant and the equilibrium timescale (%10.2es) is longer than the Hubble time.", 
			  timesc.tcmptn );
			caunin(chLine);
		}
		else
		{
			sprintf( chLine, 
				"   Compton cooling equilibrium timescale (%10.2es) is longer than Hubble time.", 
			  timesc.tcmptn );
			notein(chLine);
		}
	}

	if( timesc.ttherm > 5e17 )
	{
		sprintf( chLine, 
			" C-Thermal equilibrium timescale, %10.2es, longer than Hubble time; this cloud is not time-steady.", 
		  timesc.ttherm );
		caunin(chLine);
	}

	/* check whether model large relative to Jeans length
	 * DMEAN is mean density (gm per cc)
	 * mean temp is weighted by mass density */
	if( log10(radius.depth) > jeans_mass.rjnmin )
	{
		/* AJMIN is minimum Jeans mass, log in grams */
		aj = pow(10.,jeans_mass.ajmmin - log10(SOLAR_MASS));
		if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
		{
			sprintf( chLine, 
				" C-Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)", 
			  pow(10.f,jeans_mass.rjnmin), aj );
			caunin(chLine);
		}
		else
		{
			sprintf( chLine, 
				"   Cloud thicker than smallest Jeans length=%8.2ecm; stability problems? (smallest Jeans mass=%8.2eMo)", 
			  pow(10.f,jeans_mass.rjnmin), aj );
			notein(chLine);
		}
	}

	/* check whether grains too hot to survive */
	for( nd=0; nd < gv.nBin; nd++ )
	{
		if( gv.bin[nd]->TeGrainMax > gv.bin[nd]->Tsublimat )
		{
			sprintf( chLine, 
				" W-Maximum temperature of grain%-12.12s was%10.2eK, above its sublimation temperature,%10.2eK.", 
			  gv.bin[nd]->chDstLab, gv.bin[nd]->TeGrainMax, 
			  gv.bin[nd]->Tsublimat );
			warnin(chLine);
		}
		else if( gv.bin[nd]->TeGrainMax > gv.bin[nd]->Tsublimat* 0.9 )
		{
			sprintf( chLine, 
				" C-Maximum temperature of grain%-12.12s was%10.2eK, near its sublimation temperature,%10.2eK.", 
			  gv.bin[nd]->chDstLab, gv.bin[nd]->TeGrainMax, 
			  gv.bin[nd]->Tsublimat );
			caunin(chLine);
		}
	}

	if( negdrg.lgNegGrnDrg )
	{
		sprintf( chLine, "  !Grain drag force <0." );
		bangin(chLine);
	}

	/* largest relative number of electrons donated by grains */
	if( gv.GrnElecDonateMax > 0.05 )
	{
		sprintf( chLine, 
			"  !Grains donated %5.1f%% of the total electrons in some regions.", 
		  gv.GrnElecDonateMax*100. );
		bangin(chLine);
	}
	else if( gv.GrnElecDonateMax > 0.005 )
	{
		sprintf( chLine, 
			"   Grains donated %5.1f%% of the total electrons in some regions.", 
		  gv.GrnElecDonateMax*100. );
		notein(chLine);
	}

	/* largest relative number of electrons on grain surface */
	if( gv.GrnElecHoldMax > 0.05 )
	{
		sprintf( chLine, 
			"  !Grains contained %5.1f%% of the total electrons in some regions.", 
		  gv.GrnElecHoldMax*100. );
		bangin(chLine);
	}
	else if( gv.GrnElecHoldMax > 0.005 )
	{
		sprintf( chLine, 
			"   Grains contained %5.1f%% of the total electrons in some regions.", 
		  gv.GrnElecHoldMax*100. );
		notein(chLine);
	}

	/* is photoelectric heating of gas by photioniz of grains important */
	if( gv.dphmax > 0.5 )
	{
		sprintf( chLine, 
			"  !Local grain-gas photoelectric heating rate reached %5.1f%% of the total.", 
		  gv.dphmax*100. );
		bangin(chLine);
	}
	else if( gv.dphmax > 0.05 )
	{
		sprintf( chLine, 
			"   Local grain-gas photoelectric heating rate reached %5.1f%% of the total.", 
		  gv.dphmax*100. );
		notein(chLine);
	}

	if( gv.TotalDustHeat/thermal.power > 0.01 )
	{
		sprintf( chLine, 
			"   Global grain photoelectric heating of gas was%5.1f%% of the total.", 
		  gv.TotalDustHeat/thermal.power*100. );
		notein(chLine);
		if( gv.TotalDustHeat/thermal.power > 0.25 )
		{
			sprintf( chLine, 
				"  !Grain photoelectric heating is VERY important." );
			bangin(chLine);
		}
	}

	/* grain-gas collis cooling of gas */
	if( gv.dclmax > 0.05 )
	{
		sprintf( chLine, 
			"   Local grain-gas cooling of gas rate reached %5.1f%% of the total.", 
		  gv.dclmax*100. );
		notein(chLine);
	}

	/* check how H2 chemistry network performed */
	if( h2.renorm_max > 1.05 )
	{
		if( h2.renorm_max > 1.2 )
		{
			sprintf( chLine, 
				"  !The large H2 molecule - main chemistry network renormalization factor reached %.2f.", 
				h2.renorm_max);
			bangin(chLine);
		}
		else
		{
			sprintf( chLine, 
				"   The large H2 molecule - main chemistry network renormalization factor reached %.2f.", 
				h2.renorm_max);
			notein(chLine);
		}
	}
	if( h2.renorm_min < 0.95 )
	{
		if( h2.renorm_min < 0.8 )
		{
			sprintf( chLine, 
				"  !The large H2 molecule - main chemistry network renormalization factor reached %.2f.", 
				h2.renorm_min);
			bangin(chLine);
		}
		else
		{
			sprintf( chLine, 
				"   The large H2 molecule - main chemistry network renormalization factor reached %.2f.", 
				h2.renorm_min);
			notein(chLine);
		}
	}

	/* check whether photodissociation of H_2^+ molecular ion was important */
	if( hmi.h2pmax > 0.10 )
	{
		sprintf( chLine, 
			"  !The local H2+ photodissociation heating rate reached %5.1f%% of the total heating.", 
		  hmi.h2pmax*100. );
		bangin(chLine);
	}

	else if( hmi.h2pmax > 0.01 )
	{
		sprintf( chLine, 
			"   The local H2+ photodissociation heating rate reached %.1f%% of the total heating.", 
		  hmi.h2pmax*100. );
		notein(chLine);
	}

	/* check whether photodissociation of molecular hydrogen (H2)was important */
	if( hmi.h2dfrc > 0.1 )
	{
		sprintf( chLine, 
			"  !The local H2 photodissociation heating rate reached %.1f%% of the total heating.", 
		  hmi.h2dfrc*100. );
		bangin(chLine);
	}
	else if( hmi.h2dfrc > 0.01 )
	{
		sprintf( chLine, 
			"   The local H2 photodissociation heating rate reached %.1f%% of the total heating.", 
		  hmi.h2dfrc*100. );
		notein(chLine);
	}

	/* check whether cooling by molecular hydrogen (H2) was important */
	if( hmi.h2line_cool_frac > 0.1 )
	{
		sprintf( chLine, 
			"  !The local H2 cooling rate reached %.1f%% of the local cooling.", 
		  hmi.h2line_cool_frac*100. );
		bangin(chLine);
	}
	else if( hmi.h2line_cool_frac > 0.01 )
	{
		sprintf( chLine, 
			"   The local H2 cooling rate reached %.1f%% of the local cooling.", 
		  hmi.h2line_cool_frac*100. );
		notein(chLine);
	}

	if( hmi.h2dtot/thermal.power > 0.01 )
	{
		sprintf( chLine, 
			"   Global H2 photodissociation heating of gas was %.1f%% of the total heating.", 
		  hmi.h2dtot/thermal.power*100. );
		notein(chLine);
		if( hmi.h2dtot/thermal.power > 0.25 )
		{
			sprintf( chLine, "   H2 photodissociation heating is VERY important." );
			notein(chLine);
		}
	}

	/* check whether photodissociation of carbon monoxide (co) was important */
	if( co.codfrc > 0.25 )
	{
		sprintf( chLine, 
			"  !Local CO photodissociation heating rate reached %.1f%% of the total.", 
		  co.codfrc*100. );
		bangin(chLine);
	}
	else if( co.codfrc > 0.05 )
	{
		sprintf( chLine, 
			"   Local CO photodissociation heating rate reached %.1f%% of the total.", 
		  co.codfrc*100. );
		notein(chLine);
	}

	/* say whether CO rotation cooling was important */
	if( co.COCoolBigFrac >0.1 )
	{
		if( co.lgCOCoolCaped ) 
		{
			/* this is bad, CO cooling important and atom capped */
			sprintf( chLine, 
				" C-Local CO cooling reached %.1f%% of the local cooling, but the CO molecule was too small.", 
			  co.COCoolBigFrac*100. );
			caunin(chLine);
			sprintf( chLine, 
				" C-Increase size with ATOM CO LEVELS xxx command.  There were %li levels.", 
			  nCORotate );
			caunin(chLine);
		}
		else
		{
			/* this is good, CO cooling important and atom not capped */
			sprintf( chLine, 
				"   Local CO rotation cooling reached %.1f%% of the local cooling.", 
			  co.COCoolBigFrac*100. );
			notein(chLine);
		}
	}

	if( co.codtot/thermal.power > 0.01 )
	{
		sprintf( chLine, 
			"   Global CO photodissociation heating of gas was %.1f%% of the total.", 
		  co.codtot/thermal.power*100. );
		notein(chLine);
		if( co.codtot/thermal.power > 0.25 )
		{
			sprintf( chLine, "   CO photodissociation heating is VERY important." );
			notein(chLine);
		}
	}

	if( thermal.lgEdnGTcm )
	{
		sprintf( chLine, 
			"   Energy density of radiation field was greater than the Compton temperature. Is this physical?" );
		notein(chLine);
	}

	/* was cooling due to induced recombination important? */
	if( hydro.cintot/thermal.power > 0.01 )
	{
		sprintf( chLine, "   Induced recombination cooling was %.1f%% of the total.", 
		  hydro.cintot/thermal.power*100. );
		notein(chLine);
	}

	/* check whether free-free heating was significant */
	if( thermal.FreeFreeTotHeat/thermal.power > 0.1 )
	{
		sprintf( chLine, "  !Free-free heating was %.1f%% of the total.", 
		  thermal.FreeFreeTotHeat/thermal.power*100. );
		bangin(chLine);
	}
	else if( thermal.FreeFreeTotHeat/thermal.power > 0.01 )
	{
		sprintf( chLine, "   Free-free heating was %.1f%% of the total.", 
		  thermal.FreeFreeTotHeat/thermal.power*100. );
		notein(chLine);
	}

	/* was heating due to H- absorption important? */
	if( hmi.hmitot/thermal.power > 0.01 )
	{
		sprintf( chLine, "   H- absorption heating was %.1f%% of the total.", 
		  hmi.hmitot/thermal.power*100. );
		notein(chLine);
	}

	/* water desruction rate was zero */
	if( co.lgH2Ozer )
	{
		sprintf( chLine, "   Water destruction rate zero." );
		notein(chLine);
	}

	/* numerical instability in matrix inverson routine */
	if( atoms.nNegOI > 0 )
	{
		sprintf( chLine, " C-O I negative level populations %ld times.", 
		  atoms.nNegOI );
		caunin(chLine);
	}

	/* check for negative optical depths,
	 * optical depth in excited state helium lines */
	small = 0.;
	imas = 0;
	isav = 0;
	j = 0;
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ipLo=ipH2p; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
				{
					if( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn < (float)small )
					{
						small = EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn;
						imas = ipHi;
						j = ipLo;
						isav = nelem;
					}
				}
			}
		}
	}

	if( small < -0.05 )
	{
		sprintf( chLine, 
			"  !Some hydrogenic lines mased, species was %2s%2ld, smallest tau was %.2e, transition %li-%li", 
			elementnames.chElementSym[isav], 
			isav+1,small, imas , j );
		bangin(chLine);
	}

	/* check for negative opacities */
	if( opac.lgOpacNeg )
	{
		sprintf( chLine, "  !Some opacities were negative - the SET NEGOPC command will punch which ones." );
		bangin(chLine);
	}

	/* now check continua */
	small = 0.;
	imas = 0;
	isav = 0;
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( i=0; i < iso.numLevels[ipH_LIKE][nelem]; i++ )
			{
				if( opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][i]-1] < -0.001 )
				{
					small = MIN2(small,(double)opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][i]-1]);
					imas = i;
					isav = nelem;
				}
			}
		}
	}

	if( small < -0.05 )
	{
		sprintf( chLine, "  !Some hydrogenic (%2s%2ld) continua optical depths were negative; smallest=%10.2e level=%3ld", 
			elementnames.chElementSym[isav], 
			isav+1,
		  small, imas );
		bangin(chLine);
	}

	/* check whether any continuum optical depths are negative */
	nneg = 0;
	tauneg = 0.;
	freqn = 0.;
	for( i=0; i < rfield.nflux; i++ )
	{
		if( opac.TauAbsGeo[0][i] < -0.001 )
		{
			nneg += 1;
			/* only remember the smallest freq, and most neg optical depth */
			if( nneg == 1 )
				freqn = rfield.anu[i];
			tauneg = MIN2(tauneg,(double)opac.TauAbsGeo[0][i]);
		}
	}

	if( nneg > 0 )
	{
		sprintf( chLine, "  !Some continuous optical depths <0.  The lowest freq was %.3e Ryd, and a total of%4ld", 
		  freqn, nneg );
		bangin(chLine);
		sprintf( chLine, "  !The smallest optical depth was %.2e", 
		  tauneg );
		bangin(chLine);
	}

	/* say if Balmer continuum optically thick */
	if( opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1] > 0.05 )
	{
		sprintf( chLine, "   The Balmer continuum optical depth was%10.2e.", 
		  opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1] );
		notein(chLine);
	}

	/* was correction for stimulated emission significant? */
	if( opac.stimax[0] > 0.02 && opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] > 0.2 )
	{
		sprintf( chLine, "   The Lyman continuum stimulated emission correction to optical depths reached%10.2e.", 
		  opac.stimax[0] );
		notein(chLine);
	}
	else if( opac.stimax[1] > 0.02 && opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1] > 0.1 )
	{
		sprintf( chLine, "   The Balmer continuum stimulated emission correction to optical depths reached%10.2e.", 
		  opac.stimax[1] );
		notein(chLine);
	}

	/* say if Paschen continuum optically thick */
	if( opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][3]-1] > 0.2 )
	{
		sprintf( chLine, 
			"   The Paschen continuum optical depth was %.2e.", 
		  opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][3]-1] );
		notein(chLine);
	}

	/* some comments about near IR total optical depth */
	if( opac.TauAbsGeo[0][0] > 1. )
	{
		sprintf( chLine, 
			"   The continuum optical depth at the lowest energy considered (%.3e Ryd) was %10.3e.", 
		  rfield.anu[0], opac.TauAbsGeo[0][0] );
		notein(chLine);
	}

	/* comment if optical depth to Rayleigh scattering is big
	 * cs from VAL 76 */
	if( colden.colden[ipCHI]*7e-24 > 0.01 )
	{
		sprintf( chLine, 
			"   The optical depth to Rayleigh scattering at 1300A is %.2e", 
		  colden.colden[ipCHI]*6.71e-24 );
		notein(chLine);
	}

	if( colden.colden[ipCH2PLS]*7e-18 > 0.1 )
	{
		sprintf( chLine, 
			"  !The optical depth to the H2+ molecular ion is %.2e", 
		  colden.colden[ipCH2PLS]*7e-18 );
		bangin(chLine);
	}
	else if( colden.colden[ipCH2PLS]*7e-18 > 0.01 )
	{
		sprintf( chLine, 
			"   The optical depth to the H2+ molecular ion is %.2e", 
		  colden.colden[ipCH2PLS]*7e-18 );
		notein(chLine);
	}

	/* warn if optically thick to H- absorption */
	if( opac.thmin > 0.1 )
	{
		sprintf( chLine, 
			"  !Optical depth to negative hydrogen ion is %.2e", 
		  opac.thmin );
		bangin(chLine);
	}
	else if( opac.thmin > 0.01 )
	{
		sprintf( chLine, 
			"   Optical depth to negative hydrogen ion is %.2e", 
		  opac.thmin );
		notein(chLine);
	}

	/* say if 3-dody rec coef function outside range of validity
	 * tripped if te/z**2 < 100 or approx 10**13: => effect >50% of radiative
	 * other intergers defined in source for da */
	if( ionbal.ifail > 0 && ionbal.ifail <= 10 )
	{
		sprintf( chLine, 
			"   3 body rec coef outside range %ld", ionbal.ifail );
		notein(chLine);
	}
	else if( ionbal.ifail > 10 )
	{
		sprintf( chLine, 
			" C-3 body rec coef outside range %ld", ionbal.ifail );
		caunin(chLine);
	}

	/* check whether energy density less than background */
	if( phycon.TEnerDen < 2.6 )
	{
		sprintf( chLine, 
			"  !Incident radiation field energy density is less than 2.7K.  Add background with FIREBALL command." );
		bangin(chLine);
	}

	/* incident radiation field is less than background Habing ISM field */
	if( rfield.lgHabing )
	{
		sprintf( chLine, 
			"  !The intensity of the incident radiation field is less than 10 times the Habing diffuse ISM field.  Is this OK?" );
		bangin(chLine);
		sprintf( chLine, 
			"  !   Consider adding diffuse ISM emission with TABLE ISM command." );
		bangin(chLine);
	}

	/* check whether cosmic rays on, but model thick to them */
	if( hextra.cryden > 0. && (colden.colden[ipCHI]/10. + colden.colden[ipCHII]) > 1e23 )
	{
		sprintf( chLine, 
			" C-Model is thick to cosmic rays, which are on." );
		caunin(chLine);
	}

	/* was ionization rate less than cosmic ray ionization rate in ISM? */
	if( hextra.cryden == 0. && iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s] < 1e-17 )
	{
		sprintf( chLine, 
			"  !Ionization rate fell below background cosmic ray ionization rate.  Should this be added too?" );
		bangin(chLine);
		sprintf( chLine, 
			"  !   Use the COSMIC RAY BACKGROUND command." );
		bangin(chLine);
	}

	/* more than 10% of H is in the H2 molecule */
	if( hmi.BiggestH2 > 0.1 )
	{
		sprintf( chLine, "  !The fraction of %s in %s reached %.1f%% at some point in the cloud.", 
			"H ",
			"H2   ",
			hmi.BiggestH2*100. );
		bangin(chLine);
	}
	else if( hmi.BiggestH2>0.01 )
	{
		sprintf( chLine, "   The fraction of %s in %s reached %.2f%% at some point in the cloud.", 
			"H ",
			"H2   ",
			hmi.BiggestH2*100. );
		notein(chLine);
	}
	else if( hmi.BiggestH2 > 1e-3 )
	{
		sprintf( chLine, "   The fraction of %s in %s reached %.3f%% at some point in the cloud.", 
			"H ",
			"H2   ",
			hmi.BiggestH2*100. );
		notein(chLine);
	}

	/* largest fraction in any heavy element molecule */
	for( i=0; i<NUM_HEAVY_MOLEC; ++i )
	{
		if( co.xMoleFracMax[i] > 0.1 )
		{
			sprintf( chLine, "  !The fraction of %s in %s reached %.1f%% at some point in the cloud.", 
				elementnames.chElementSym[co.nelem_hevmol[i]],
				co.chLab[i],
				co.xMoleFracMax[i]*100. );
			bangin(chLine);
		}
		else if( co.xMoleFracMax[i]>0.01 )
		{
			sprintf( chLine, "   The fraction of %s in %s reached %.2f%% at some point in the cloud.", 
				elementnames.chElementSym[co.nelem_hevmol[i]],
				co.chLab[i],
				co.xMoleFracMax[i]*100. );
			notein(chLine);
		}
		else if( co.xMoleFracMax[i] > 1e-3 )
		{
			sprintf( chLine, "   The fraction of %s in %s reached %.3f%% at some point in the cloud.", 
				elementnames.chElementSym[co.nelem_hevmol[i]],
				co.chLab[i],
				co.xMoleFracMax[i]*100. );
			notein(chLine);
		}
	}

	/* brems optical depth */
	if( rfield.EnergyBremsThin > 0.09 )
	{
		sprintf( chLine, "  !The cloud is optically thick at optical wavelengths, extending to %.3e Ryd =%10.3eA", 
		  rfield.EnergyBremsThin, RYDLAM/rfield.EnergyBremsThin );
		bangin(chLine);
	}
	else if( rfield.EnergyBremsThin > 0.009 )
	{
		sprintf( chLine, "   The continuum of the computed structure may be optically thick in the near infrared." );
		notein(chLine);
	}

	/* did model run away to very large radius? */
	if( radius.Radius > 1e23 && radius.Radius/radius.rinner > 10. )
	{
		sprintf( chLine, "   Is an outer radius of %.2e reasonable?", 
		  radius.Radius );
		notein(chLine);
	}

	/* following set true in RT_line_one_tauinc if maser capped at tau = -1 */
	if( rt.lgMaserCapHit )
	{
		sprintf( chLine, "   Laser maser optical depths capped in RT_line_one_tauinc." );
		notein(chLine);
	}

	/* following set true in adius_next if maser cap set dr */
	if( rt.lgMaserSetDR )
	{
		sprintf( chLine, "  !Line maser set zone thickness in some zones." );
		bangin(chLine);
	}

	/* lgPradCap is true if radiation pressure was capped on first iteration
	 * also check that this is a constant total pressure model */
	if( (pressure.lgPradCap && (strcmp(dense.chDenseLaw,"CPRE") == 0)) && 
	  pressure.lgRadPresON )
	{
		sprintf( chLine, "   Radiation pressure kept below gas pressure on this iteration." );
		notein(chLine);
	}

	if( pressure.RadBetaMax > 0.25 )
	{
		if( pressure.ipPradMax_line == 0 )
		{
			sprintf( chLine, 
				"  !The ratio of radiation to gas pressure reached %.2e at zone %li.  Caused by Lyman alpha.", 
			  pressure.RadBetaMax,
			  pressure.ipPradMax_nzone);
			bangin(chLine);
		}
		else
		{
			sprintf( chLine, 
				"  !The ratio of radiation to gas pressure reached %.2e at zone %li.  "
				"Caused by line number %ld, label %s", 
			  pressure.RadBetaMax, 
			  pressure.ipPradMax_nzone,
			  pressure.ipPradMax_line,
			  pressure.chLineRadPres );
			bangin(chLine);
		}
	}

	else if( pressure.RadBetaMax > 0.025 )
	{
		if( pressure.ipPradMax_line == 0 )
		{
			sprintf( chLine, 
				"   The ratio of radiation to gas pressure reached %.2e at zone %li.  Caused by Lyman alpha.", 
			  pressure.RadBetaMax,
			  pressure.ipPradMax_nzone);
			notein(chLine);
		}
		else
		{
			sprintf( chLine, 
				"   The ratio of radiation to gas pressure reached %.2e at zone %li.  "
				"Caused by line number %ld, label %s", 
			  pressure.RadBetaMax, 
			  pressure.ipPradMax_nzone,
			  pressure.ipPradMax_line,
			  pressure.chLineRadPres );
			notein(chLine);
		}
	}

	if( opac.telec >= 5. )
	{
		sprintf( chLine, " W-The model is optically thick to electron scattering; tau=%.2e", 
		  opac.telec );
		warnin(chLine);
	}
	else if( opac.telec > 2.0  )
	{
		sprintf( chLine, " C-The model is moderately optically thick to electron scattering; tau=%.1f", 
		  opac.telec );
		caunin(chLine);
	}
	else if( opac.telec > 0.1  )
	{
		sprintf( chLine, "  !The model has modest optical depth to electron scattering; tau=%.2f", 
		  opac.telec );
		bangin(chLine);
	}
	else if( opac.telec > 0.01 )
	{
		sprintf( chLine, "   The optical depth to electron scattering is %.3f", 
		  opac.telec );
		notein(chLine);
	}

	/* optical depth to 21 cm */
	if( HFLines[0].TauIn > 0.5 )
	{
		sprintf( chLine, "  !The optical depth in the H I 21 cm line is %.2e",HFLines[0].TauIn );
		bangin(chLine);
	}

	/* optical depth in the CO 1-0 transition */
	if( C12O16Rotate[0].TauIn > 0.5 )
	{
		sprintf( chLine, "  !The optical depth in the 12CO J=1-0 line is %.2e",C12O16Rotate[0].TauIn );
		bangin(chLine);
	}
	if( C13O16Rotate[0].TauIn > 0.5 )
	{
		sprintf( chLine, "  !The optical depth in the 13CO J=1-0 line is %.2e",C13O16Rotate[0].TauIn );
		bangin(chLine);
	}

	/* check on optical depth convergence of all hydrogenic lines */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			if( EmisLines[ipH_LIKE][nelem][3][ipH2p].TauIn > 0.2 )
			{
				differ = fabs(1.-EmisLines[ipH_LIKE][nelem][3][ipH2p].TauIn*
				  rt.DoubleTau/EmisLines[ipH_LIKE][nelem][3][ipH2p].TauTot)*100.;

				/* check whether H-alpha optical depth changed by much on last iteration
				 * no tolerance can be finer than autocv, the tolerance on the
				 * iterate to convergence comamnd.  It is 15% */
				if( ((iterations.lgLastIt && EmisLines[ipH_LIKE][nelem][3][ipH2p].TauIn > 0.8) && 
					differ > 20.) && wind.windv == 0. )
				{
					sprintf( chLine, 
						" C-This is the last iteration and %2s%2ld Bal(a) optical depth"
						" changed by%6.1f%% (was %10.2e). Try another iteration.", 
					  elementnames.chElementSym[nelem], 
					  nelem+1, differ, 
					  EmisLines[ipH_LIKE][nelem][3][ipH2p].TauTot );
					caunin(chLine);
					iterations.lgIterAgain = TRUE;
				}

				/* only check on Lya convergence if Balmer lines are thick */
				if( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn > 0. )
				{
					differ = fabs(1.-EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn*
					  rt.DoubleTau/EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot)*100.;

					/* check whether Lya optical depth changed on last iteration
					 * no tolerance can be finer than autocv, the tolerance on the
					 * iterate to convergence comamnd.  It is 15% */
					if( ((iterations.lgLastIt && EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn > 0.8) && 
						differ > 25.) && wind.windv == 0. )
					{
						sprintf( chLine, 
							" C-This is the last iteration and %2s%2ld Ly(a) optical depth"
							" changed by%7.0f%% (was %10.2e). Try another iteration.", 
						elementnames.chElementSym[nelem], 
						  nelem+1,differ, EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot );
						caunin(chLine);
						iterations.lgIterAgain = TRUE;
					}
				}
			}
		}
	}

	/* check whether sphere was set if dr/r large */
	if( radius.Radius/radius.rinner > 2. && !geometry.lgSphere )
	{
		sprintf( chLine, " C-R(out)/R(in)=%10.2e and SPHERE was not set.", 
		  radius.Radius/radius.rinner );
		caunin(chLine);
	}

	/* check if thin in hydrogen or helium continua, but assumed to be thick */
	if( iterations.lgLastIt && !opac.lgCaseB )
	{

		/* check if thin in Lyman continuum, and assumed thick */
		if( rfield.nflux > iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] )
		{
			if( opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] < 2. && 
				opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1] > 2. )
			{
				sprintf( chLine, " C-The Lyman continuum is thin, and I assumed"
					" that it was thick.  Try another iteration." );
				caunin(chLine);
				iterations.lgIterAgain = TRUE;
			}
		}

		/* check on the He+ ionizing continuum */
		if( rfield.nflux > iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s] )
		{
			if( (opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1] < 2. && 
				 opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1] > 2.) )
			{
				sprintf( chLine, 
					" C-The He II continuum is thin and I assumed that it was thick."
					"  Try another iteration." );
				caunin(chLine);
				iterations.lgIterAgain = TRUE;
			}
		}

		if( rfield.nflux > iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] )
		{
			if( (opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1] < 2. && 
				 opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1] > 2.) )
			{
				sprintf( chLine, 
					" C-The He I continuum is thin and I assumed that it was thick."
					"  Try another iteration." );
				caunin(chLine);
				iterations.lgIterAgain = TRUE;
			}
		}
	}

	/* check whether column density changed by much on this iteration */
	if( iteration > 1 )
	{
		if( colden.colden_old[ipCOL_HTOT] <= 0. )
		{
			fprintf( ioQQQ, " colden_old is insane in PrtComment.\n" );
			ShowMe();
			puts( "[Stop in PrtComment]" );
			ShowMe();
			cdEXIT(EXIT_FAILURE);
		}

		differ = fabs(1.-colden.colden[ipCOL_HTOT]/
			colden.colden_old[ipCOL_HTOT]);

		if( differ > 0.1 && differ <= 0.3 )
		{
			sprintf( chLine, 
				"   The H column density changed by %.2e%% between this and previous iteration.", 
			  differ*100. );
			notein(chLine);
		}

		else if( differ > 0.3 )
		{
			if( iterations.lgLastIt )
			{
				sprintf( chLine, 
					" C-The H column density changed by %.2e%% and this is the last iteration.  What happened?", 
				  differ*100. );
				caunin(chLine);
			}
			else
			{
				sprintf( chLine, 
					"  !The H column density changed by %.2e%%  What happened?", 
				  differ*100. );
				bangin(chLine);
			}
		}

		/* check on H2 column density, but only if significant fraction of H is molecular */
		if( colden.colden[ipCOLH2]/colden.colden[ipCOL_HTOT] > 1e-5 )
		{
			differ = fabs(1.-colden.colden[ipCOLH2]/
				SDIV(colden.colden_old[ipCOLH2]));

			if( differ > 0.1 && differ <= 0.3 )
			{
				sprintf( chLine, 
					"   The H2 column density changed by %.2e%% between this and previous iteration.", 
				differ*100. );
				notein(chLine);
			}

			else if( differ > 0.3 )
			{
				if( iterations.lgLastIt )
				{
					sprintf( chLine, 
						" C-The H2 column density changed by %.2e%% and this is the last iteration.  What happened?", 
					differ*100. );
					caunin(chLine);
				}
				else
				{
					sprintf( chLine, 
						"  !The H2 column density changed by %.2e%%  What happened?", 
					differ*100. );
					bangin(chLine);
				}
			}
		}
	}

	/* say if rad pressure caused by la and la optical depth changed too much */
	differ = fabs(1.-EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn/
	  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot)*100.;

	if( iterations.lgLastIt && (pressure.RadBetaMax > 0.1) && 
		(differ > 50.) && (pressure.ipPradMax_line == 1) && (pressure.lgRadPresON) && 
		wind.windv == 0. )
	{
		sprintf( chLine, " C-This is the last iteration, radiation pressure was significant, and the L-a optical depth changed by %7.2f%% (was%10.2e)", 
		  differ, EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot );
		caunin(chLine);
	}

	/* say if la rad pressure capped by thermalization length */
	if( pressure.lgPradDen )
	{
		sprintf( chLine, "   Line radiation pressure capped by thermalization length." );
		notein(chLine);
	}

	/* print te failures */
	nline = MIN2(conv.nTeFail,10);
	if( conv.nTeFail != 0 )
	{
		if( conv.failmx < 0.1 )
		{
			_o = sprintf( chLine, "   There were %ld minor temperature failures.  zones:", 
				conv.nTeFail );
			for( i=0; i < nline; i++ )
			{
				_o += sprintf( chLine+_o, " %ld", conv.ifailz[i] );
			}
			notein(chLine);
		}
		else
		{
			_o = sprintf( chLine, 
				"  !There were %ld temperature failures, and some were large. The largest was %.1f%%.  What happened?", 
			  conv.nTeFail, conv.failmx*100. );
			bangin(chLine);

			_o = sprintf( chLine , "  !The zones were" );
			for( i=0; i < nline; i++ )
			{
				_o += sprintf( chLine+_o, " %ld", conv.ifailz[i] );
			}
			bangin(chLine);

			if( struc.testr[0] > 8e4 && phycon.te < 6e5 )
			{
				sprintf( chLine, "  !I think they may have been caused by the change from hot to nebular gas phase.  The physics of this is unclear." );
				bangin(chLine);
			}
		}
	}

	/* check for temperature jumps */
	big_ion_jump = 0.;
	j = 0;
	for( i=1; i < nzone; i++ )
	{
		big = fabs(1.-struc.testr[i-1]/struc.testr[i]);
		if( big > big_ion_jump )
		{
			j = i;
			big_ion_jump = big;
		}
	}

	if( big_ion_jump > 0.2 )
	{
		/* this is a sanity check, but only do it if jump detected */
		if( j < 1 )
		{
			fprintf( ioQQQ, " j too small bigjump check\n" );
			puts( "[Stop in PrtComment]" );
			ShowMe();
			cdEXIT(EXIT_FAILURE);
		}

		if( big_ion_jump > 0.4 )
		{
			sprintf( chLine, " C-A temperature discontinuity occurred at zone %ld from %.2eK to %.2eK.", 
			  j, struc.testr[j-1], struc.testr[j] );
			caunin(chLine);
		}
		else if( big_ion_jump > 0.2 )
		{
			sprintf( chLine, "  !A temperature discontinuity occurred at zone %ld from %.2eK to %.2eK.", 
			  j, struc.testr[j-1], struc.testr[j] );
			bangin(chLine);
		}
	}

	/* check for largest error in local electron density */
	if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed )
	{
		/* this only produces a warning if not the very last zone */
		if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed*20. && dense.nzEdenBad != 
		  nzone )
		{
			sprintf( chLine, " W-The local error in the electron density reached %.1f%% at zone %ld", 
			  conv.BigEdenError*100, dense.nzEdenBad );
			warnin(chLine);
		}
		else if( fabs(conv.BigEdenError) > conv.EdenErrorAllowed*5. )
		{
			sprintf( chLine, " C-The local error in the electron density reached %.1f%% at zone %ld", 
			  conv.BigEdenError*100, dense.nzEdenBad );
			caunin(chLine);
		}
		else
		{
			sprintf( chLine, "   The local error in the electron density reached %.1f%% at zone %ld", 
			  conv.BigEdenError*100, dense.nzEdenBad );
			notein(chLine);
		}
	}

	/* check for temperature oscillations or fluctuations*/
	big_ion_jump = 0.;
	j = 0;
	for( i=1; i < (nzone - 1); i++ )
	{
		big = fabs( (struc.testr[i-1] - struc.testr[i])/struc.testr[i] );
		bigm = fabs( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] );

		/* this is sign of change in temperatre, we are looking for change in sign */
		rel = ( (struc.testr[i-1] - struc.testr[i])/struc.testr[i])*
			( (struc.testr[i] - struc.testr[i+1])/struc.testr[i] );

		if( rel < 0. && MIN2( bigm , big ) > big_ion_jump )
		{
			j = i;
			big_ion_jump = MIN2( bigm , big );
		}
	}

	if( big_ion_jump > 0.1 )
	{
		/* only do sanity check if jump detected */
		if( j < 1 )
		{
			fprintf( ioQQQ, " j too small bigjump2 check\n" );
			puts( "[Stop in PrtComment]" );
			ShowMe();
			cdEXIT(EXIT_FAILURE);
		}

		if( big_ion_jump > 0.3 )
		{
			sprintf( chLine, 
				" C-A temperature oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", 
			  j, big_ion_jump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] );
			caunin(chLine);
		}
		else if( big_ion_jump > 0.1 )
		{
			sprintf( chLine, 
				"  !A temperature oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", 
			  j, big_ion_jump*100., struc.testr[j-1], struc.testr[j], struc.testr[j+1] );
			bangin(chLine);
		}
	}

	/* check for eden oscillations */
	if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
	{
		j = 0;
		big_ion_jump = 0.;
		for( i=1; i < (nzone - 1); i++ )
		{
			big = (struc.ednstr[i-1] - struc.ednstr[i])/struc.ednstr[i];
			if( fabs(big) < conv.EdenErrorAllowed )
				big = 0.;
			bigm = (struc.ednstr[i] - struc.ednstr[i+1])/struc.ednstr[i];
			if( fabs(bigm) < conv.EdenErrorAllowed )
				bigm = 0.;
			if( big*bigm < 0. && 
				fabs(struc.ednstr[i-1]-struc.ednstr[i])/struc.ednstr[i] > big_ion_jump )
			{
				j = i;
				big_ion_jump = fabs(struc.ednstr[i-1]-struc.ednstr[i])/
				  struc.ednstr[i];
			}
		}
		if( j < 1 )

		/* only check on j if there was a big jump detected, number must be
		 * smallest jump */
		if( big_ion_jump > conv.EdenErrorAllowed*3. )
		{
			if( j < 1 )
			{
				fprintf( ioQQQ, " j too small bigjump3 check\n" );
				puts( "[Stop in PrtComment]" );
				ShowMe();
				cdEXIT(EXIT_FAILURE);
			}

			if( big_ion_jump > conv.EdenErrorAllowed*10. )
			{
				sprintf( chLine, " C-An electron density oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", 
				  j, big_ion_jump*100., struc.ednstr[j-1], struc.ednstr[j], 
				  struc.ednstr[j+1] );
				caunin(chLine);
			}
			else if( big_ion_jump > conv.EdenErrorAllowed*3. )
			{
				sprintf( chLine, "  !An electron density oscillation occurred at zone%4ld by%5.0f%% from%10.2e to %10.2e to%10.2e", 
				  j, big_ion_jump*100., struc.ednstr[j-1], struc.ednstr[j], 
				  struc.ednstr[j+1] );
				bangin(chLine);
			}
		}
	}

	/*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
	/* >>chng 03 dec 05, add this test */
	prt_smooth_predictions();

	/**********************************************************
	 * check that the volume integrates out ok                *
	 **********************************************************/

	/* this was the number 1 fed through the line integrators,
	 * the number 1e-10 is sent to linadd in lineset1 as follows:*/
	/*linadd( 1.e-10 , 1 , "Unit" , 'i' );*/
	i = cdLine( "Unit" , 1 , &rate , &absint );
	ASSERT( i> 0 );

	/* this is now the linear vol, rel to inner radius */
	VolComputed = LineSv[i].sumlin /  1e-10 ;

	/* >>chng 02 apr 22, do not zero this element out, so can be used to get vol */
	/* set stored number to zero so it does not appear on the emission line list
	LineSv[i].sumlin = 0.; */

	/* spherical or plane parallel case? */
	if( radius.Radius/radius.rinner > 1.0001 )
	{
		/* spherical case, 
		 * geometry.iEmissPower is usually 2,
		 * and can be reset to 1 (long slit) or 0 (beam) with 
		 * slit and beam options on aperture */
		VolExpected = geometry.covgeo*geometry.FillFac*radius.rinner/(geometry.iEmissPower+1)*
			( powi( radius.Radius/radius.rinner,geometry.iEmissPower+1 ) - 1. );
	}
	else
	{
		/* plane parallel case */
		/* next can be zero for very thin model, depth is always positive */
		VolExpected = geometry.covgeo*geometry.FillFac*(radius.depth-DEPTH_OFFSET);
	}

	/* now get the relative difference between computed and expected volumns */
	error = fabs(VolComputed - VolExpected)/VolExpected;

	/* we need to ignore this test if filling factor changes with radius, or
	 * cylinder geometry in place */
	if( radius.lgCylnOn || geometry.filpow!=0. )
	{
		error = 0.;
	}

	/* how large is relative error? */
	if( error > 0.001 )
	{
		sprintf( chLine, 
			" W-PrtComment insanity - Line unit integration did not verify \n");
		warnin(chLine);
		fprintf( ioQQQ,
			" PROBLEM PrtComment insanity - Line unit integration did not verify \n");
		fprintf( ioQQQ,
			" expected, derived vols were %g %g \n",
			VolExpected , VolComputed );
		fprintf( ioQQQ,
			" relative difference is %g, ratio is %g.\n",error,VolComputed/VolExpected);
		TotalInsanity();
	}

	/* next do same thing for fake continuum point propagated in highest energy cell, plus 1 
	 *  = 
	 * the variable rfield.ConEmitLocal[rfield.nflux]
	 * are set to 
	 * the number 1.e-10f * Dilution in RT_diffuse.  this is the outward
	 * local emissivity, per unit vol.  It is then added to the outward beams
	 * by the rest of the code, and then checked here.
	 *
	 * insanity will be detected if diffuse emission is thrown into the outward beam
	 * in MadeDiffuse.  this happens if the range of ionization encompasses the full
	 * continuum array, up to nflux.  */
	ConComputed = rfield.ConInterOut[rfield.nflux]/ 1e-10;
	/* correct for fraction that went out, as set in ZoneStart,
	 * this should now be the volume of the emitting region */
	ConComputed /= ( (1. + geometry.covrt)/2. ) ;

	/* we expect this to add up to the integral of unity over r^-2 */
	if( radius.Radius/radius.rinner < 1.001 )
	{
		/* plane parallel case, no dilution, use thickness */
		ConExpected = (radius.depth-DEPTH_OFFSET)*geometry.FillFac;
	}
	else
	{
		/* spherical case */
		ConExpected = radius.rinner*geometry.FillFac * (1. - radius.rinner/radius.Radius );
	}
	/* this is impossible */
	ASSERT( ConExpected > 0. );

	/* now get the relative difference between computed and expected volumns */
	error = fabs(ConComputed - ConExpected)/ConExpected;

	/* we need to ignore this test if filling factor changes with radius, or
	 * cylinder geometry in place */
	if( radius.lgCylnOn || geometry.filpow!=0. )
	{
		error = 0.;
	}

	/* how large is relative error? */
	if( error > 0.001 )
	{
		sprintf( chLine, 
			" W-PrtComment insanity - Continuum unit integration did not verify \n");
		warnin(chLine);
		fprintf( ioQQQ," PROBLEM PrtComment insanity - Continuum unit integration did not verify \n");
		fprintf( ioQQQ," exact vol= %g, derived vol= %g relative difference is %g \n",
			ConExpected , ConComputed ,error);
		fprintf( ioQQQ," ConInterOut= %g,  \n",
			rfield.ConInterOut[rfield.nflux]);
		TotalInsanity();
	}

	/* final prinout of warnings, cautions, notes, */
	cdNwcns(&lgAbort,&nw,&nc,&nn,&ns,&i,&j,&dum1,&dum2);

	warnings.lgWarngs = ( nw > 0 );
	warnings.lgCautns = ( nc > 0 );

	if( called.lgTalk )
	{
		/* print the title of the calculation */
		fprintf( ioQQQ, "   %s\n", input.chTitle  );
		/* say why the calculation stopped, and indicate the geometry*/
		cdReasonGeo(ioQQQ);
		/* print all warnings */
		cdWarnings(ioQQQ);
		/* all cautions */
		cdCautions(ioQQQ);
		/* surprises, begining with a ! */
		cdSurprises(ioQQQ);
		/* notes about the calculations */
		cdNotes(ioQQQ);
	}

	/* option to print warnings on special io */
	if( lgPrnErr )
	{
		cdWarnings(ioPrnErr);
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " PrtComment returns.\n" );
	}

#	ifdef DEBUG_FUN
	fputs( " <->PrtComment()\n", debug_fp );
#	endif
	return;
}

/*badprt print out coolants if energy not conserved */
static void badprt(
		/* total is total luminosity available in radiation */
		double total)
{
	/* following is smallest ratio to print */
#	define	RATIO	0.02
	char chInfo;
	long int i;
	float sum_coolants, 
	  sum_recom_lines;

#	ifdef DEBUG_FUN
	fputs( "<+>badprt()\n", debug_fp );
#	endif

	fprintf( ioQQQ, " badprt: all entries with greater than%6.2f%% of incident continuum follow.\n", 
	  RATIO*100. );

	fprintf( ioQQQ, " badprt: Intensities are relative to total energy radiated.\n" );

	/* now find sum of recombination lines */
	chInfo = 'r';
	sum_recom_lines = (float)totlin('r');
	fprintf( ioQQQ, 
		" Sum of energy in recombination lines is %.2e relative to total incident radiation is %.2e\n", 
		sum_recom_lines, 
		sum_recom_lines/MAX2(1e-30,total) );

	fprintf(ioQQQ," all strong information lines \n line  wl  ener/total\n");
	/* now print all strong lines */
	for( i=0; i < LineSave.nsum; i++ )
	{
		if( LineSv[i].chSumTyp == chInfo && LineSv[i].sumlin/total > RATIO )
		{
			fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
			prt_wl( ioQQQ, LineSv[i].wavelength );
			fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].sumlin/total, chInfo );
		}
	}

	fprintf(ioQQQ," all strong cooling lines \n line  wl  ener/total\n");
	chInfo = 'c';
	sum_coolants = (float)totlin('c');
	fprintf( ioQQQ, " Sum of coolants (abs) =%10.2e (rel)=%10.2e\n", 
	  sum_coolants, sum_coolants/MAX2(1e-30,total) );
	for( i=0; i < LineSave.nsum; i++ )
	{
		if( LineSv[i].chSumTyp == chInfo && LineSv[i].sumlin/total > RATIO )
		{
			fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
			prt_wl( ioQQQ, LineSv[i].wavelength );
			fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].sumlin/total, chInfo );
		}
	}

	fprintf(ioQQQ," all strong heating lines \n line  wl  ener/total\n");
	chInfo = 'h';
	fprintf( ioQQQ, " Sum of heat (abs) =%10.2e (rel)=%10.2e\n", 
	  thermal.power, thermal.power/MAX2(1e-30,total) );
	for( i=0; i < LineSave.nsum; i++ )
	{
		if( LineSv[i].chSumTyp == chInfo && LineSv[i].sumlin/total > RATIO )
		{
			fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
			prt_wl( ioQQQ, LineSv[i].wavelength );
			fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].sumlin/total, chInfo );
		}
	}


#	ifdef DEBUG_FUN
	fputs( " <->badprt()\n", debug_fp );
#	endif
	return;
}

/*chkCaHeps check whether CaII K and H epsilon overlap */
static void chkCaHeps(double *totwid)
{
	double conca, 
		conalog ,
	  conhe;

#	ifdef DEBUG_FUN
	fputs( "<+>chkCaHeps()\n", debug_fp );
#	endif

	*totwid = 0.;
	/* pumping of CaH ovelapping with Hy eps, 6-2 of H */
	if( iso.numLevels[ipH_LIKE][ipHYDROGEN] > 6 )
	{
		if( TauLines[ipT3969].TauIn > 0. && 
			EmisLines[ipH_LIKE][ipHYDROGEN][6][2].TauIn >   0. )
		{
			conca = pow(6.1e-5* TauLines[ipT3969].TauIn,0.5);
			conalog = log(TauLines[ipT3969].TauIn);
			conalog = sqrt(MAX2(1.,	conalog));
			conca = MAX2(conalog,conca);

			conalog = log(EmisLines[ipH_LIKE][ipHYDROGEN][6][2].TauIn);
			conalog = sqrt(MAX2(1.,conalog));
			conhe = pow(1.7e-6*EmisLines[ipH_LIKE][ipHYDROGEN][6][2].TauIn,0.5);
			conhe = MAX2(conalog, conhe);

			*totwid = 10.*conhe + 1.6*conca;
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->chkCaHeps()\n", debug_fp );
#	endif
	return;
}

/*outsum sum outward continuum beams */
static void outsum(double *outtot, 
  double *outin, 
  double *outout)
{
	long int i;

#	ifdef DEBUG_FUN
	fputs( "<+>outsum()\n", debug_fp );
#	endif

	*outin = 0.;
	*outout = 0.;
	for( i=0; i < rfield.nflux; i++ )
	{
		/* N.B. in following en1ryd prevents overflow */
		*outin += rfield.anu[i]*(rfield.flux[i]*EN1RYD);
		*outout += rfield.anu[i]*(rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])*
		  EN1RYD;
	}

	*outtot = *outin + *outout;

#	ifdef DEBUG_FUN
	fputs( " <->outsum()\n", debug_fp );
#	endif
	return;
}

/*prt_smooth_predictions check whether fluctuations in any predicted quantities occurred */
static void prt_smooth_predictions(void)
{
	long int i,
		nzone_oscillation,
		nzone_ion_jump,
		nzone_den_jump,
		nelem,
		ion;
	double BigOscillation ,
		big_ion_jump,
		big_jump,
		rel,
		big,
		bigm;

	char chLine[INPUT_LINE_LENGTH];

#	ifdef DEBUG_FUN
	fputs( "<+>prt_smooth_predictions()\n", debug_fp );
#	endif

	/* check for ionization oscillations or fluctuations and or jumps */
	nzone_oscillation = 0;
	nzone_ion_jump = 0;

	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ion=0; ion<=nelem+1; ++ion) 
			{
				BigOscillation = 0.;
				big_ion_jump = -15.;
				for( i=1; i < (nzone - 1); i++ )
				{

					/* only do check if all ions are positive */
					if( struc.xIonDense[nelem][ion][i-1]/struc.gas_phase[nelem][i-1]>struc.dr_ionfrac_limit &&
						struc.xIonDense[nelem][ion][i  ]/struc.gas_phase[nelem][i  ]>struc.dr_ionfrac_limit &&
						struc.xIonDense[nelem][ion][i+1]/struc.gas_phase[nelem][i+1]>struc.dr_ionfrac_limit )
					{

						/* this is check for oscillations */
						big = fabs( (struc.xIonDense[nelem][ion][i-1] - struc.xIonDense[nelem][ion][i])/struc.xIonDense[nelem][ion][i] );
						bigm = fabs( (struc.xIonDense[nelem][ion][i]  - struc.xIonDense[nelem][ion][i+1])/struc.xIonDense[nelem][ion][i] );

						/* this is sign of change in ionization, we are looking for change in sign */
						rel = ( (struc.xIonDense[nelem][ion][i-1] - struc.xIonDense[nelem][ion][i]  )/struc.xIonDense[nelem][ion][i])*
							  ( (struc.xIonDense[nelem][ion][i]   - struc.xIonDense[nelem][ion][i+1])/struc.xIonDense[nelem][ion][i] );

						if( rel < 0. && MIN2( bigm , big ) > BigOscillation )
						{
							nzone_oscillation = i;
							BigOscillation = MIN2( bigm , big );
						}

						/* check whether we tripped over an ionizaiton front - a major source
						 * of instability in a complete linearization code like this one */
						/* neg sign picks up only increases in ionization */
						rel = -log10( (struc.xIonDense[nelem][ion][i]/struc.gas_phase[nelem][i]) / 
							(struc.xIonDense[nelem][ion][i+1]/struc.gas_phase[nelem][i+1] ) );
						/* only do sig stages of ionization */
						if( rel > big_ion_jump )
						{
							big_ion_jump = rel;
							nzone_ion_jump = i;
						}
					}
				}
				/* end loop over zones, 
				 * check whether this ion and element underwent fluctuations or jump */

				if( BigOscillation > 0.2 )
				{
					/* only do sanity check if jump detected */
					if( nzone_oscillation < 1 )
					{
						fprintf( ioQQQ, " nzone_oscillation too small bigjump2 check\n" );
						puts( "[Stop in prt_smooth_predictions]" );
						ShowMe();
						cdEXIT(EXIT_FAILURE);
					}
					if( BigOscillation > 3. )
					{
						sprintf( chLine, 
							" W-An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e", 
							nzone_oscillation, 
							elementnames.chElementSym[nelem], ion+1,
							BigOscillation*100., 
							struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1], 
							struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation], 
							struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
							warnin(chLine);
					}

					else if( BigOscillation > 0.7 )
					{
						sprintf( chLine, 
							" C-An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e", 
							nzone_oscillation, 
							elementnames.chElementSym[nelem], ion+1,
							BigOscillation*100., 
							struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1], 
							struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation], 
							struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
							caunin(chLine);
					}
					else if( BigOscillation > 0.2 )
					{
						sprintf( chLine, 
							"  !An ionization oscillation occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e", 
							nzone_oscillation, 
							elementnames.chElementSym[nelem], ion+1,
							BigOscillation*100., 
							struc.xIonDense[nelem][ion][nzone_oscillation-1]/struc.gas_phase[nelem][nzone_oscillation-1], 
							struc.xIonDense[nelem][ion][nzone_oscillation]/struc.gas_phase[nelem][nzone_oscillation], 
							struc.xIonDense[nelem][ion][nzone_oscillation+1]/struc.gas_phase[nelem][nzone_oscillation+1] );
							bangin(chLine);
					}
				}

				/* big_ion_jump was a log above, convert to linear quantity */
				/* if no jump occurred then big_ion_jump is small and nzone_ion_jump is 0 */
				big_ion_jump = pow(10., big_ion_jump );
				if( big_ion_jump > 1.5 && nzone_ion_jump > 0 )
				{
					if( big_ion_jump > 10. )
					{
						sprintf( chLine, 
							" C-An ionization jump occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e", 
							nzone_ion_jump, 
							elementnames.chElementSym[nelem], ion+1,
							big_ion_jump*100., 
							struc.xIonDense[nelem][ion][nzone_ion_jump-1]/struc.gas_phase[nelem][nzone_ion_jump-1], 
							struc.xIonDense[nelem][ion][nzone_ion_jump]/struc.gas_phase[nelem][nzone_ion_jump], 
							struc.xIonDense[nelem][ion][nzone_ion_jump+1]/struc.gas_phase[nelem][nzone_ion_jump+1] );
							caunin(chLine);
					}
					else
					{
						sprintf( chLine, 
							"  !An ionization jump occurred at zone %ld, elem %.2s%2li, by %.0f%% from %.2e to %.2e to %.2e", 
							nzone_ion_jump, 
							elementnames.chElementSym[nelem], ion+1,
							big_ion_jump*100., 
							struc.xIonDense[nelem][ion][nzone_ion_jump-1]/struc.gas_phase[nelem][nzone_ion_jump-1], 
							struc.xIonDense[nelem][ion][nzone_ion_jump]/struc.gas_phase[nelem][nzone_ion_jump], 
							struc.xIonDense[nelem][ion][nzone_ion_jump+1]/struc.gas_phase[nelem][nzone_ion_jump+1] );
							bangin(chLine);
					}
				}
			}
		}
	}

	big_jump = -15;
	nzone_den_jump = 0;

	for( i=1; i < (nzone - 1); i++ )
	{
		/* this first check is on how the total hydrogen density has changed */
		rel = fabs(log10( struc.gas_phase[ipHYDROGEN][i] / 
			struc.gas_phase[ipHYDROGEN][i+1] ) );
		/* only do sig stages of ionization */
		if( rel > big_jump )
		{
			big_jump = rel;
			nzone_den_jump = i;
		}
	}

	/* check how stable density was */
	big_jump = pow( 10., big_jump );
	if( big_jump > 1.2 )
	{
		if( big_jump > 3. )
		{
			sprintf( chLine, 
				" C-The H density jumped at by %.0f%% at zone %ld, from %.2e to %.2e to %.2e", 
				big_jump*100., 
				nzone_den_jump, 
				struc.gas_phase[ipHYDROGEN][nzone_den_jump-1], 
				struc.gas_phase[ipHYDROGEN][nzone_den_jump], 
				struc.gas_phase[ipHYDROGEN][nzone_den_jump+1] );
				caunin(chLine);
		}
		else
		{
			sprintf( chLine, 
				"  !An H density jump occurred at zone %ld, by %.0f%% from %.2e to %.2e to %.2e", 
				nzone_den_jump, 
				big_jump*100., 
				struc.gas_phase[ipHYDROGEN][nzone_den_jump-1], 
				struc.gas_phase[ipHYDROGEN][nzone_den_jump], 
				struc.gas_phase[ipHYDROGEN][nzone_den_jump+1] );
				bangin(chLine);
		}
	}

	/* now do check on smoothness of radiation pressure */
	big_jump = -15;
	nzone_den_jump = 0;

	/* loop starts on zone 3 since dramatic fall in radiation pressure across first
	 * few zones is normal behavior */
	for( i=3; i < (nzone - 2); i++ )
	{
		/* this first check is on how the total hydrogen density has changed */
		rel = fabs(log10( SDIV(struc.PresRadCurr[i]) / 
			SDIV(0.5*(struc.PresRadCurr[i-1]+struc.PresRadCurr[i+1])) ) );
		/* only do sig stages of ionization */
		if( rel > big_jump )
		{
			big_jump = rel;
			nzone_den_jump = i;
		}
	}

	/* check how stable radiation pressure was, but only if signicant */
	if( pressure.RadBetaMax > 0.01 )
	{
		big_jump = pow( 10., big_jump );
		if( big_jump > 1.2 )
		{
			/* only make it a caution is pressure jumped, and we were trying
			* to do a constant pressure model */
			if( big_jump > 3. && strcmp(dense.chDenseLaw,"CPRE") == 0)
			{
				sprintf( chLine, 
					" C-The radiation pressure jumped at by %.0f%% at zone %ld, from %.2e to %.2e to %.2e", 
					big_jump*100., 
					nzone_den_jump, 
					struc.PresRadCurr[nzone_den_jump-1], 
					struc.PresRadCurr[nzone_den_jump], 
					struc.PresRadCurr[nzone_den_jump+1] );
					caunin(chLine);
			}
			else
			{
				sprintf( chLine, 
					"  !The radiation pressure jumped at by %.0f%% at zone %ld, from %.2e to %.2e to %.2e", 
					big_jump*100., 
					nzone_den_jump, 
					struc.PresRadCurr[nzone_den_jump-1], 
					struc.PresRadCurr[nzone_den_jump], 
					struc.PresRadCurr[nzone_den_jump+1] );
					bangin(chLine);
			}
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->prt_smooth_predictions()\n", debug_fp );
#	endif
	return;
}
