/* 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 */
/*coolr main routine to call others, to evaluate total cooling */
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "taulines.h"
#include "wind.h"
#include "coolheavy.h"
#include "radius.h"
#include "converge.h"
#include "h2.h"
#include "opacity.h"
#include "ionbal.h"
#include "dense.h"
#include "trace.h"
#include "dynamics.h"
#include "rfield.h"
#include "grainvar.h"
#include "atoms.h"
#include "called.h"
#include "hmi.h"
#include "numderiv.h"
#include "magnetic.h"
#include "phycon.h"
#include "bit32.h"
#include "lines_service.h"
#include "hyperfine.h"
#include "iso.h"
#include "thermal.h"
#include "cooling.h"
/*fndneg search cooling array to find negative values */
static void fndneg(void);
/*fndstr search cooling stack to find strongest values */
static void fndstr(double tot, 
  double dc);

/* set true to debug deriv of heating and cooling */
#define	PRT_DERIV	FALSE

void coolr(double *tot)
{
	long int ion,
	  i,
	  ipISO,
	  nelem;

	static long int nhit = 0, 
	  nzSave=0;

	static float TeEvalCS = 0.;
	static float TeUsedBrems=-1.f;
	static int nzoneUsedBrems=-1;

	double cs ,
	  deriv, 
	  factor, 
	  qn, 
	  rothi=-SMALLFLOAT, 
	  rotlow=-SMALLFLOAT, 
	  x;

	static double oltcool=0., 
	  oldtemp=0.;

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

	/* returns tot, the total cooling,
	 * and dc, the derivative of the cooling */

	/* routine atom_level2( t10 )
	 * routine atom_level3( abund , t10,t21,t20)
	 * tsq1 = 1. / (te**2)
	 * POPEXC( O12,g1,g2,A21,excit,abund) ; result already*a21
	 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2)
	 * AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41)
	 * FIVEL( G(1-5) , ex(wn,1-5), cs12,cs13,14,15,23,24,25,34,35,45,
	 *  A21,31,41,51,32,42,52,43,53,54, pop(1-5), abund) */

	if( trace.lgTrace )
		fprintf( ioQQQ, "   COOLR TE:%.4e zone %li %li Cool:%.4e Heat:%.4e eden:%.4e edenTrue:%.4e\n", 
		phycon.te, 
		nzone, conv.nPres2Ioniz ,
		thermal.ctot , thermal.htot,dense.eden,dense.EdenTrue );

	/* now zero out the cooling stack */
	CoolZero();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  0 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
	if( gv.lgGrainPhysicsOn )
	{
		/* grain heating and cooling */
		/* grain recombination cooling, evaluated elsewhere
		* can either heat or cool the gas, do cooling here */
		CoolAdd("dust",0,MAX2(0.,gv.GasCoolColl));

		/* grain cooling proportional to temperature ^3/2 */
		thermal.dCooldT += MAX2(0.,gv.GasCoolColl)*3./(2.*phycon.te);

		/* these are the various heat agents from grains */
		thermal.heating[0][13] = 0.;
		thermal.heating[0][14] = 0.;
		thermal.heating[0][25] = 0.;

		/* options to force gas heating or cooling by grains to zero - for tests only ! */
		if( gv.lgDustOn && gv.lgDHetOn )
		{
			long int nd;
			/* if grains hotter than gas then collisions with gas act
			* to heat the gas, add this in here
			* a symmetric statement appears in COOLR, where cooling is added on */
			thermal.heating[0][14] = MAX2(0.,-gv.GasCoolColl);

			/* now add up heating agents over all grain bins */
			for( nd=0; nd < gv.nBin; nd++ )
			{

				/* rate dust heats gas by photoelectric effect */
				thermal.heating[0][13] += gv.bin[nd]->GasHeatPhotoEl;

				/* thermionic is gas heating due to thermionic emissions */
				thermal.heating[0][25] += gv.bin[nd]->thermionic;

			}
		}
	}

	/* find net heating - cooling due to large H2 molecule */
	H2_Cooling();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  1 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* molecular molecules molecule cooling */
	if( hmi.lgNoH2Mole )
	{
		hmi.hmicol = 0.;
		CoolHeavy.brems_cool_hminus = 0.;
		/* line cooling within simple H2 molecule - zero when big used */
		CoolHeavy.h2line = 0.;
		/*  H + H+ => H2+ cooling */
		CoolHeavy.H2PlsCool = 0.;
		CoolHeavy.HD = 0.;

		/* thermal.heating[0][8] is heating due to collisions within X of H2 */
		thermal.heating[0][8] = 0.;
		/* thermal.heating[0][15] is H minus heating*/
		thermal.heating[0][15] = 0.;
		/* thermal.heating[0][16] is H2+ heating */
		thermal.heating[0][16] = 0.;
		hmi.HeatH2Dish_used = 0.;
		hmi.HeatH2Dexc_used = 0.;
	}

	else
	{

		/* save it */
		thermal.heating[0][15] = hmi.hmihet;
		thermal.heating[0][16] = hmi.h2plus_heat;
		/* now get heating from H2 molecule, either simple or from big one */
		if( h2.lgH2ON  && hmi.lgBigH2_evaluated && hmi.lgH2_Thermal_BigH2 )
		{
			/* these are explictly from big H2 molecule,
			 * first is heating due to radiative pump of excited states, followed by
			 * radiative decay into continuum of X, followed by dissociation of molecule
			 * with kinetic energy, typically 0.25 - 0.5 eV per event */
			hmi.HeatH2Dish_used = hmi.HeatH2Dish_BigH2;
			hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BigH2;
			/*fprintf(ioQQQ,"DEBUG big %.2f\t%.5e\t%.2e\t%.2e\t%.2e\n", 
				fnzone , phycon.te, hmi.HeatH2Dexc_used,
				hmi.H2_total , hmi.H2_star_BigH2);*/
			/* negative sign because right term is really deriv of heating,
			 * but will be used below as deriv of cooling */
			hmi.deriv_HeatH2Dexc_used = -hmi.deriv_HeatH2Dexc_BigH2;
		}
		else
		{
			/* these come from approximations in TH85, see comments above */
			hmi.HeatH2Dish_used = hmi.HeatH2Dish_TH85;
			hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_TH85;
			hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_TH85;
		}

		/* heating due to photodissociation heating */
		thermal.heating[0][17] = hmi.HeatH2Dish_used;

		/* heating (usually cooling in big H2) due to collisions within X */
		/* add to heating is net heating is positive */
		thermal.heating[0][8] = MAX2(0.,hmi.HeatH2Dexc_used);
		/* add to cooling if net heating is negative */
		CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used));
		/*fprintf(ioQQQ,"DEBUG coolh2\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
			fnzone, phycon.te, dense.eden, hmi.H2_total, thermal.ctot, -hmi.HeatH2Dexc_used );*/
		/* add to net derivative */
		/*thermal.dCooldT += MAX2(0.,-hmi.HeatH2Dexc_used)* ( 30172. * thermal.tsq1 - thermal.halfte );*/
		/* >>chng 04 jan 25, check sign to prevent cooling from entering here,
		 * also enter neg sign since going into cooling stack (bug), in headsum
		 * same term adds to deriv of heating */
		if( hmi.HeatH2Dexc_used < 0. )
			thermal.dCooldT -= hmi.deriv_HeatH2Dexc_used;

		/*  H + H+ => H2+ cooling */
		CoolHeavy.H2PlsCool = (float)(MAX2((2.325*phycon.te-1875.)*1e-20,0.)*
		  dense.xIonDense[ipHYDROGEN][0]*dense.xIonDense[ipHYDROGEN][1]*1.66e-11);

		if( h2.lgH2ON )
		{
			/* this is simplified approximation to H2 rotation cooling,
			 * big molecule goes this far better */
			CoolHeavy.h2line = 0.;
		}
		else
		{
			/*  rate for rotation lines from 
			*  >>refer	h2	cool	Lepp, S., & Shull, J.M. 1983, ApJ, 270, 578 */
			x = phycon.alogte - 4.;
			if( phycon.te > 1087. )
			{
				rothi = 3.90e-19*sexp(6118./phycon.te);
			}
			else
			{
				rothi = pow(10.,-19.24 + 0.474*x - 1.247*x*x);
			}

			/*  low density rotation cooling */
			/*&qn = pow(MAX2(hmi.Hmolec[ipMH2g],1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);*/
			qn = pow(MAX2(hmi.H2_total,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);
			if( phycon.te > 4031. )
			{
				rotlow = 1.38e-22*sexp(9243./phycon.te)*qn;
			}
			else
			{
				rotlow = pow(10.,-22.90 - 0.553*x - 1.148*x*x)*qn;
			}

			if( rotlow > 0. )
			{
				/*CoolHeavy.h2line = hmi.Hmolec[ipMH2g]*rothi/(1. + rothi/rotlow);*/
				CoolHeavy.h2line = hmi.H2_total*rothi/(1. + rothi/rotlow);
			}
			else
			{
				CoolHeavy.h2line = 0.;
			}
		}

		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC && nzone>187&& iteration > 1/**/)
			{
				fprintf(ioQQQ,"h2coolbug\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
					phycon.te, 
					CoolHeavy.h2line, 
					hmi.H2_total, 
					hmi.Hmolec[ipMHm], 
					hmi.HMinus_photo_rate,
					rothi,
					rotlow );
			}
		}
		
		/* >>chng 02 mar 07, add DH cooling using rates (eqn 6) from
		 * >>refer	HD	cooling	Puy, D., Grenacher, L, & Jetzer, P., 1999, A&A, 345, 723 */
		factor = sexp(128.6/phycon.te);
		/*CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2(hmi.Hmolec[ipMH2g]) * phycon.sqrte *
			factor/(1416.+phycon.sqrte*hmi.Hmolec[ipMH2g] * (1. + 3.*factor));*/
		CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2(hmi.H2_total) * phycon.sqrte *
			/*factor/(1416.+phycon.sqrte*hmi.Hmolec[ipMH2g] * (1. + 3.*factor));*/
			factor/(1416.+phycon.sqrte*hmi.H2_total * (1. + 3.*factor));
	}


	/*  H- FB; H + e -> H- + hnu */
	/*  H- FF is in with H ff */
	CoolAdd("H-fb",0,hmi.hmicol);
	/* >>chng 96 nov 15, fac of 2 in deriv to help convergence in very dense
	 * models where H- is important, this takes change in eden into
	 * partial account */
	thermal.dCooldT += 2.*hmi.hmicol*phycon.teinv;
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  2 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	CoolAdd("H2ln",0,CoolHeavy.h2line);
	/* >>chng 00 oct 21, added coef of 3.5, sign had been wrong */
	/*thermal.dCooldT += CoolHeavy.h2line*phycon.teinv;*/
	/* >>chng 03 mar 17, change 3.5 to 15 as per behavior in primal.in */
	/*thermal.dCooldT += 3.5*CoolHeavy.h2line*phycon.teinv;*/
	/* >>chng 03 may 18, from 15 to 30 as per behavior in primal.in - overshoots happen */
	/*thermal.dCooldT += 15.*CoolHeavy.h2line*phycon.teinv;*/
	/*>>chng 03 oct 03, from 30 to 3, no more overshoots in primalin */
	/*thermal.dCooldT += 30.*CoolHeavy.h2line*phycon.teinv;*/
	thermal.dCooldT += 3.0*CoolHeavy.h2line*phycon.teinv;

	{
		/* propblems with H2 cooling */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC /*&& nzone>300 && iteration > 1*/)
		{
			fprintf(ioQQQ,"coolr debuggg\t%.2f\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n",
			fnzone, 
			phycon.te, 
			hmi.H2_total , 
			CoolHeavy.h2line,
			hmi.Hmolec[ipMHm] , 
			dense.eden);
		}
	}

	CoolAdd("HDro",0,CoolHeavy.HD);
	thermal.dCooldT += CoolHeavy.HD*phycon.teinv;

	CoolAdd("H2+ ",0,CoolHeavy.H2PlsCool);
	thermal.dCooldT += CoolHeavy.H2PlsCool*phycon.teinv;

	/* cooling due to C12O16 */
	CoolAdd("CO C",12,CoolHeavy.C12O16Rot);
	thermal.dCooldT += CoolHeavy.dC12O16Rot; 
	/* >>chng 00 oct 25, add C13O16 cooling */
	/* cooling due to C13O16 */
	CoolAdd("CO C",13,CoolHeavy.C13O16Rot);
	thermal.dCooldT += CoolHeavy.dC13O16Rot; 

	/* heating due to three-body, will be incremented in iso_cool*/
	thermal.heating[0][3] = 0.;
	/* heating due to hydrogen lines */
	thermal.heating[0][23] = 0.;
	/* heating due to photoionization of all excited states of hydrogen species */
	thermal.heating[0][1] = 0.;

	/* isoelectronic species cooling, mainly lines, and level ionization */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			/* must always call iso_cool since we must zero those variables
			* that would have been set had the species been present */
			iso_cool( ipISO , nelem );
		}
	}

	/*lint -e777 float test equals */
	/* >>chng 02 jun 18, don't reevaluate needlessly */
	/* >>chng 03 nov 28, even faster - special logic for when ff is pretty
	 * negligible - eval of ff is pretty slow */
	/* >>chng 04 feb 19, must not test on temp not changing, since ionization
	 * can change at constant temperature 
	 * >>chng 04 sep 14, above introduced bug since brems never reevaluated
	 * now test is zone or temp has changed */
	if( (fabs(CoolHeavy.brems_cool_net)/SDIV(thermal.ctot)>0.001 ) || 
		conv.lgSearch || 
		fabs((phycon.te-TeUsedBrems)/phycon.te-1.)>0.01 ||
		nzone!=nzoneUsedBrems )
	/*lint +e777 float test equals */
	{
		double BremsThisIon;
		double OpacityThisIon;
		long int limit;
		/* free-free free free brems emission for all ions */

		TeUsedBrems = phycon.te;
		nzoneUsedBrems = nzone;
		/* highest frequency where we have non-zero botzmann factors */
		limit = MIN2( rfield.ipMaxBolt , rfield.nflux );

		CoolHeavy.brems_cool_hminus = 0.;
		CoolHeavy.brems_cool_h = 0.;
		CoolHeavy.brems_cool_metals = 0.;
		CoolHeavy.brems_cool_he = 0.;
		CoolHeavy.brems_heat_total = 0.;

		{
			static double *TotBremsMetalsAllIons;
			double bfac = (dense.eden/1e20)/phycon.sqrte/1e10, bhfac, ofac;
			static int lgFirstTime=TRUE;

			ASSERT(rfield.ipEnergyBremsThin < rfield.nupper);
			ASSERT(limit < rfield.nupper);
			if( lgFirstTime )
			{
				/* rfield.nupper will not change in one coreload, so just malloc this once */
				TotBremsMetalsAllIons = (double *)MALLOC((unsigned long)rfield.nupper*sizeof(double));
				lgFirstTime = FALSE;
			}
			for( i=0; i<limit; i++ )
			{
				TotBremsMetalsAllIons[i] = 0.;
			}
			
			/* the optically thick part */
			nelem = ipHYDROGEN;
			ion = 1;
			bhfac = bfac*dense.xIonDense[nelem][ion];
			for( i=0; i<rfield.ipEnergyBremsThin; i++ )
			{
				/* first do hydrogen opacity */
				OpacityThisIon =  bhfac*(1. - rfield.ContBoltz[i])*rfield.gff[ion][i] * opac.OpacStack[i-1+opac.ipBrems];
				opac.FreeFreeOpacity[i] = OpacityThisIon;

				/* now h minus */
				OpacityThisIon *= iso.Pop2Ion[ipH_LIKE][nelem][ipH1s] * opac.OpacStack[i-1+opac.iphmra];
				opac.FreeFreeOpacity[i] += OpacityThisIon;
			}	

			/* helium */
			/* now do helium, both He+ and He++ */
			nelem = ipHELIUM;

			for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
			{
				ofac = bfac*POW2( (double)ion ) *(dense.xIonDense[nelem][ion]);
				for( i=0; i<rfield.ipEnergyBremsThin; i++ )
				{
					OpacityThisIon = ofac*(1. - rfield.ContBoltz[i])*rfield.gff[ion][i] * opac.OpacStack[i-1+opac.ipBrems];
					opac.FreeFreeOpacity[i] += OpacityThisIon;
				}
			}

			/* hydrogen and helium were done above */
			for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
			{
				/* for H, IonLow is 0 and IonHigh is 1 */
				/* MAX2 occurs because we want to start at first ion (or above)
				 * and do not want atom */
				for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
				{
					ofac = bfac* POW2( (double)ion )*(dense.xIonDense[nelem][ion]);
					for( i=0; i<rfield.ipEnergyBremsThin; i++ )
					{
						OpacityThisIon = ofac*(1. - rfield.ContBoltz[i])*rfield.gff[ion][i] * opac.OpacStack[i-1+opac.ipBrems];
						opac.FreeFreeOpacity[i] += OpacityThisIon;
					}
				}
			}
			
			/* ipEnergyBremsThin is index to energy where gas goes optically thin to brems,
			 * so this loop is over optically thin frequencies 
			 * do not add optically thick part to outward beam since self absorbed */

			/* do hydrogen first, before main loop since want to break out as separate
			 * coolant, and what to add on H- brems */
			nelem = ipHYDROGEN;
			ion = 1;
			
			bhfac = POW2( (double)ion )*dense.xIonDense[nelem][ion]*CoolHeavy.lgFreeOn* dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
			ofac = bfac*(dense.xIonDense[nelem][ion]);
			for( i=rfield.ipEnergyBremsThin; i < limit; i++ )
			{
				
				/* in all following CoolHeavy.lgFreeOn is flag set with 'no free-free' to
				 * turn off brems heating and cooling */
				BremsThisIon = bhfac*rfield.gff[ion][i]*rfield.widflx[i]*rfield.ContBoltz[i];
				ASSERT( BremsThisIon >= 0. );
				CoolHeavy.brems_cool_h += BremsThisIon;
				
				OpacityThisIon =  ofac*rfield.gff[ion][i] *
					(1. - rfield.ContBoltz[i])* opac.OpacStack[i-1+opac.ipBrems];
				opac.FreeFreeOpacity[i] = OpacityThisIon;
				
				/* for case of hydrogen, do H- brems - OpacStack contains the ratio
				 * of the H- to H brems cross section - multiply by this and the fraction in ground */
				CoolHeavy.brems_cool_hminus += BremsThisIon * opac.OpacStack[i-1+opac.iphmra]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
				
				OpacityThisIon *= iso.Pop2Ion[ipH_LIKE][nelem][ipH1s]*opac.OpacStack[i-1+opac.iphmra];
				opac.FreeFreeOpacity[i] += OpacityThisIon;
			}

			/* now do helium, both He+ and He++ */
			nelem = ipHELIUM;
			for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
			{
				bhfac = POW2( (double)ion )*dense.xIonDense[nelem][ion]*CoolHeavy.lgFreeOn * dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
				ofac = POW2( (double)ion )*bfac*(dense.xIonDense[nelem][ion]);
				for( i=rfield.ipEnergyBremsThin; i < limit; i++ )
				{
					/* eff. charge is ion, so first rfield.gff argument must be "ion".	*/
					BremsThisIon = bhfac*rfield.gff[ion][i]*rfield.widflx[i]*rfield.ContBoltz[i];
					ASSERT( BremsThisIon >= 0. );				
					CoolHeavy.brems_cool_he += BremsThisIon;

					OpacityThisIon = ofac*(1. - rfield.ContBoltz[i])*rfield.gff[ion][i]*opac.OpacStack[i-1+opac.ipBrems];
					opac.FreeFreeOpacity[i] += OpacityThisIon;
				}
			}

			/* hydrogen and helium were done above */
			for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
			{
				/* for H, IonLow is 0 and IonHigh is 1 */
				/* MAX2 occurs because we want to start at first ion (or above)
				 * and do not want atom */
				for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
				{
					/* eff. charge is ion, so first rfield.gff argument must be "ion".	*/
					
					bhfac = POW2( (double)ion )*dense.xIonDense[nelem][ion];
					ofac = bfac*POW2( (double)ion )*(dense.xIonDense[nelem][ion]);
					for( i=rfield.ipEnergyBremsThin; i < limit; i++ )
					{
						BremsThisIon = bhfac*rfield.gff[ion][i];
						/*ASSERT( rfield.gff[ion][i] > 0. );*/
						ASSERT( BremsThisIon >= 0. );
						TotBremsMetalsAllIons[i] += BremsThisIon;

						OpacityThisIon = ofac*(1. - rfield.ContBoltz[i])*rfield.gff[ion][i] * 
							opac.OpacStack[i-1+opac.ipBrems];
						opac.FreeFreeOpacity[i] += OpacityThisIon;
					}
				}
			}
			
			for( i=rfield.ipEnergyBremsThin; i<limit;  i++ )
			{
				/* now convert to physical units */
				TotBremsMetalsAllIons[i] *= dense.eden*1.032e-11*rfield.widflx[i]*rfield.ContBoltz[i]/phycon.sqrte*EN1RYD*CoolHeavy.lgFreeOn;
				ASSERT( TotBremsMetalsAllIons[i] >= 0.);
				
				/* the total cooling due to brems */
				CoolHeavy.brems_cool_metals += TotBremsMetalsAllIons[i];
				
				/* the total heating due to brems */
				CoolHeavy.brems_heat_total += opac.FreeFreeOpacity[i]*rfield.flux[i]*rfield.anu[i]*EN1RYD*CoolHeavy.lgFreeOn;
			}
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC && nzone>60 /*&& iteration > 1*/)
				{
					double sumfield = 0., sumtot=0., sum1=0., sum2=0.;
					for( i=rfield.ipEnergyBremsThin; i<limit;  i++ )
					{
						sumtot += opac.FreeFreeOpacity[i]*rfield.flux[i]*rfield.anu[i];
						sumfield += rfield.flux[i]*rfield.anu[i];
						sum1 += opac.FreeFreeOpacity[i]*rfield.flux[i]*rfield.anu[i];
						sum2 += opac.FreeFreeOpacity[i]*rfield.flux[i];
					}
					fprintf(ioQQQ,"DEBUG brems heat\t%.2f\t%.3e\t%.3e\t%.3e\t%e\t%.3e\t%.3e\n",
						fnzone,
						CoolHeavy.brems_heat_total,
						sumtot/SDIV(sumfield) ,
						sum1/SDIV(sum2),
						phycon.te , 
						rfield.gff[1][1218],
						opac.FreeFreeOpacity[1218]);
				}
			}
		}
	}

	
	/* these two terms are both large, nearly canceling, near lte */
	CoolHeavy.brems_cool_net = 
		CoolHeavy.brems_cool_h + 
		CoolHeavy.brems_cool_he + 
		CoolHeavy.brems_cool_hminus + 
		CoolHeavy.brems_cool_metals - 
		CoolHeavy.brems_heat_total;
	/*fprintf(ioQQQ,"DEBUG brems\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
		fnzone,
		phycon.te,
		CoolHeavy.brems_cool_net,
		CoolHeavy.brems_cool_h ,
		CoolHeavy.brems_cool_he ,
		CoolHeavy.brems_cool_hminus,
		CoolHeavy.brems_cool_metals ,
		CoolHeavy.brems_heat_total);*/

	/* net free free brems cooling, count as cooling if positive */
	CoolAdd( "FF c" , 0, MAX2(0.,CoolHeavy.brems_cool_net) );

	/* now stuff into heating array if negative */
	thermal.heating[0][11] = MAX2(0.,-CoolHeavy.brems_cool_net);

	/* >>chng 96 oct 30, from HFFNet to just FreeFreeCool,
	 * since HeatSum picks up CoolHeavy.brems_heat_total */
	thermal.dCooldT += CoolHeavy.brems_cool_h*thermal.halfte;
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  3 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* >>chng 02 jun 21, net cooling already includes this */
	/* end of brems cooling */

	/* heavy element recombination cooling, do not count hydrogenic since
	 * already done above, also helium singlets have been done */
	/* >>chng 02 jul 21, put in charge dependent rec term */
	CoolHeavy.heavfb = 0.;
	for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			/*>>chng 02 jul 21, upper bound now uses NISO,
			 * note that detailed iso seq atoms are done in iso_cool */
			/* >>chng 03 sep 03, do not do zeroed stages of ionization */
			long limit_lo = MAX2( 1 , dense.IonLow[nelem] );
			/* there is a -1 on the NISO since loop will be up to <= not < */
			long limit_hi = MIN2( nelem-NISO+1-1 , dense.IonHigh[nelem] );
			for( ion=limit_lo; ion<=limit_hi; ++ion )
			/*for( ion=1; ion < nelem-NISO+1 ; ion++ )*/
			{
				/* factor of 0.9 is roughly correct for nebular conditions, see
				 * >>refer	H	rec cooling	LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */
				/* note that ionbal.RadRecomRateCoef is the rate coef, cm3 s-1, needs eden */
				/* >>chng 02 nov 07, move rec arrays around, this now has ONLY rad rec,
				 * previously had included grain rec and three body */
				/* recombination cooling for iso-seq atoms are done in iso_cool */
				double one = dense.xIonDense[nelem][ion] * ionbal.RadRecomRateCoef[nelem][ion-1]*
					dense.eden * phycon.te * BOLTZMANN;
				/*fprintf(ioQQQ,"debugggfb\t%li\t%li\t%.3e\t%.3e\t%.3e\n", nelem, ion, one, 
					dense.xIonDense[nelem][ion] , ionbal.RadRecomRateCoef[nelem][ion]);*/
				CoolHeavy.heavfb += one;
			}
		}
	}

	/*fprintf(ioQQQ,"debuggg hvFB\t%i\t%.2f\t%.2e\t%.2e\n",iteration, fnzone,CoolHeavy.heavfb, dense.eden);*/

	CoolAdd("hvFB",0,CoolHeavy.heavfb);
	thermal.dCooldT += CoolHeavy.heavfb*.113*phycon.teinv;
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  4 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* electron-electron brems, approx form from 
	 * >>refer	ee	brems	Stepney and Guilbert, MNRAS 204, 1269 (1983)
	 * ok for T<10**9 */
	CoolHeavy.eebrm = POW2(dense.eden*phycon.te*1.84e-21);

	/* >>chng 97 mar 12, added deriv */
	thermal.dCooldT += CoolHeavy.eebrm*thermal.halfte;
	CoolAdd("eeff",0,CoolHeavy.eebrm);

	/* add advective heating and cooling */
	/* this is cooling due to loss of matter from this region */
	CoolAdd("adve",0,dynamics.Cool );
	/* >>chng02 dec 04, rm factor of 8 in front of dCooldT */
	thermal.dCooldT += dynamics.dCooldT;
	/* local heating due to matter moving into this location */
	thermal.heating[1][5] = dynamics.Heat;
	thermal.dHeatdT += dynamics.dHeatdT;

	/* total Compton cooling */
	CoolHeavy.tccool = rfield.cmcool*phycon.te;
	CoolAdd("Comp",0,CoolHeavy.tccool);
	thermal.dCooldT += rfield.cmcool;

	/* option for "extra" cooling, expressed as power-law */
	if( thermal.lgCExtraOn )
	{
		CoolHeavy.cextxx = 
			(float)(thermal.CoolExtra*pow(phycon.te/1e4f,thermal.cextpw));
	}
	else
	{
		CoolHeavy.cextxx = 0.;
	}
	CoolAdd("Extr",0,CoolHeavy.cextxx);

	/* cooling due to wind expansion, only for winds expansion cooling */
	if( wind.windv > 0. )
	{
		dynamics.dDensityDT = (float)(wind.AccelTot/wind.windv + 2.*wind.windv/
		  radius.Radius);
		CoolHeavy.expans = dense.pden*phycon.te*BOLTZMANN*dynamics.dDensityDT;
	}
	else
	{
		dynamics.dDensityDT = 0.;
		CoolHeavy.expans = 0.;
	}
	CoolAdd("Expn",0,CoolHeavy.expans);
	thermal.dCooldT += CoolHeavy.expans/phycon.te;

	/* cyclotron cooling */
	/* coef is 4/3 /8pi /c * vtherm(elec) */
	CoolHeavy.cyntrn = 4.5433e-25f*magnetic.pressure*PI8*dense.eden*phycon.te;
	CoolAdd("Cycl",0,CoolHeavy.cyntrn);
	thermal.dCooldT += CoolHeavy.cyntrn/phycon.te;

	/* heavy element collisional ionization
	 * derivative should be zero since increased coll ion rate
	 * decreases neutral fraction by proportional amount */
	CoolAdd("Hvin",0,CoolHeavy.colmet);

	/* all the fine structure lines added by Will Goddard */
	if( fabs(phycon.te-TeEvalCS)/phycon.te > 0.05 )
	{

		/* H 21 cm emission/population,
		* cs will be sum of e cs and H cs converted from rate */
		cs = (H21cm_electron( phycon.te )*dense.eden + H21cm_H_atom( phycon.te )* dense.xIonDense[ipHYDROGEN][0]) * 3./	phycon.cdsqte;
		PutCS(  cs , &HFLines[0] );

		/* H 21 cm done above, now loop over remaining lines to get collision strengths */
		for( i=1; i < nHFLines; i++ )
		{
			cs = HyperfineCS( i );
			/* now generate the collision strength and put into the line array */
			PutCS(  cs , &HFLines[i] );
		}
		TeEvalCS = phycon.te;
	}

	/* do level pops for H 21 cm which is a special case since Lya pumping in included */
	H21_cm_pops(&HFLines[0] );
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  5 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* find total cooling due to hyperfine structure lines */
	hyperfine.cooling_total = HFLines[0].cool;

	/* now do level pops for all except 21 cm  */
	for( i=1; i < nHFLines; i++ )
	{
		/* remember current gas-phase abundances */
		float save = dense.xIonDense[HFLines[i].nelem-1][HFLines[i].IonStg-1];

		/* bail if no abundance */
		if( save<=0. ) continue;

		/* set gas-phase abundance to total times isotope ratio */
		dense.xIonDense[HFLines[i].nelem-1][HFLines[i].IonStg-1] *= hyperfine.HFLabundance[i] ;

		/* use the collision strength generated above and find pops and cooling */
		atom_level2( &HFLines[i] );

		/* put the correct gas-phase abundance back in the array */
		dense.xIonDense[HFLines[i].nelem-1][HFLines[i].IonStg-1] = save;

		/* find total cooling due to hyperfine structure lines */
		hyperfine.cooling_total += HFLines[i].cool;
	}
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  6 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Carbon cooling */
	CoolCarb();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  C %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Nitrogen cooling */
	CoolNitr();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  N %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Oxygen cooling */
	CoolOxyg();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  7 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Fluorine cooling */
	CoolFluo();

	/* Neon cooling */
	CoolNeon();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt Ne %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Magnesium cooling */
	CoolMagn();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  8 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Sodium cooling */
	CoolSodi();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt Na %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Aluminum cooling */
	CoolAlum();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt Al %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Silicon cooling */
	CoolSili();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt  9 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Phosphorus */
	CoolPhos();

	/* Sulphur cooling */
	CoolSulf();

	/* Chlorine cooling */
	CoolChlo();

	/* Argon cooling */
	CoolArgo();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt 10 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Potasium cooling */
	CoolPota();

	/* Calcium cooling */
	CoolCalc();

	/* Scandium cooling */
	CoolScan();

	/* Titanium cooling */
	CoolTita();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt 11 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Vanadium cooling */
	CoolVana();

	/* Chromium cooling */
	CoolChro();

	/* Iron cooling */
	CoolIron();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt 12 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* Manganese cooling */
	CoolMang();

	/* Cobalt cooling */
	CoolCoba();

	/* Nickel cooling */
	CoolNick();

	/* Zinc cooling */
	CoolZinc();
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt 13 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* do all the thousands of lines Dima added with g-bar approximation */
	CoolDima();

	/* now add up all the coolants */
	CoolSum(tot);
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dcdt 14 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);

	/* negative cooling */
	if( *tot <= 0. )
	{
		fprintf( ioQQQ, " COOLR; cooling is <=0, this is impossible.\n" );
		ShowMe();
		puts( "[Stop in coolr]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* bad derivative */
	if( thermal.dCooldT == 0. )
	{
		fprintf( ioQQQ, " COOLR; cooling slope <=0, this is impossible.\n" );
		if( (*tot > 0. && bit32.lgBit32) && dense.gas_phase[ipHYDROGEN] < 1e-4 )
		{
			fprintf( ioQQQ, " Probably due to very low density.\n" );
		}
		ShowMe();
		puts( "[Stop in coolr]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( trace.lgTrace )
	{
		fndstr(*tot,thermal.dCooldT);
	}

	/* lgTSetOn true for constant temperature model */
	if( (((((!thermal.lgTSetOn) && *tot < 0.) && called.lgTalk) && 
	  thermal.lgColNeg) && thermal.lgCNegChk) && nzone > 0 )
	{
		fprintf( ioQQQ, 
			" Negative cooling, zone%4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n", 
		  nzone, 
		  *tot, 
		  iso.cLya_cool[ipH_LIKE][ipHYDROGEN], 
		  iso.coll_ion[ipH_LIKE][ipHYDROGEN], 
		  phycon.te );
		fndneg();
	}

	/* possibility of getting emperical cooling derivative
	 * normally false, set true with 'set numerical derivatives' command */
	if( NumDeriv.lgNumDeriv )
	{
		/*lint -e777 */
		if( ((nzone > 2 && nzone == nzSave) && oldtemp != phycon.te) && nhit > 4 )
		/*lint +e777 */
		{
			/* hnit is number of trys on this zone - use to stop numerical problems
			 * do not evaluate numerical deriv until well into solution */
			deriv = (oltcool - *tot)/(oldtemp - phycon.te);
			thermal.dCooldT = deriv;
		}
		else
		{
			deriv = thermal.dCooldT;
		}
		if( nzone != nzSave )
			nhit = 0;

		nzSave = nzone;
		nhit += 1;
		oltcool = *tot;
		oldtemp = phycon.te;
	}

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

/*  */
#ifdef EPS
#	undef EPS
#endif
#define	EPS	0.01

/*fndneg search cooling array to find negative values */
static void fndneg(void)
{
	long int i;
	double trig;

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

	trig = fabs(thermal.htot*EPS);
	for( i=0; i < thermal.ncltot; i++ )
	{
		if( thermal.cooling[i] < 0. && fabs(thermal.cooling[i]) > trig )
		{
			fprintf( ioQQQ, " negative line=%s %.2f fraction of heating=%.3f\n", 
			  thermal.chClntLab[i], thermal.collam[i], thermal.cooling[i]/
			  thermal.htot );
		}

		if( thermal.heatnt[i] > trig )
		{
			fprintf( ioQQQ, " heating line=%s %.2f fraction of heating=%.3f\n", 
			  thermal.chClntLab[i], thermal.collam[i], thermal.heatnt[i]/
			  thermal.htot );
		}
	}


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

/*fndstr search cooling stack to find strongest values */
static void fndstr(double tot, 
  double dc)
{
	char chStrngLab[5];
	long int i;
	float wl;
	double str, 
	  strong;

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

	strong = 0.;
	wl = -FLT_MAX;
	for( i=0; i < thermal.ncltot; i++ )
	{
		if( fabs(thermal.cooling[i]) > strong )
		{
			/* this is the wavelength of the coolant, 0 for a continuum*/
			wl = thermal.collam[i];
			/* make sure labels are all valid*/
			ASSERT( strlen( thermal.chClntLab[i] ) == 4 );
			strcpy( chStrngLab, thermal.chClntLab[i] );
			strong = fabs(thermal.cooling[i]);
		}
	}

	str = strong;

	fprintf( ioQQQ, 
		"   fndstr cool: TE=%10.4e Ne=%10.4e C=%10.3e dC/dT=%10.2e ABS(%s %.1f)=%.2e nz=%ld\n", 
	  phycon.te, dense.eden, tot, dc, chStrngLab
	  , wl, str, nzone );

	/* option for extensive printout of lines */
	if( trace.lgCoolTr )
	{
		float ratio;

		/* flag all significant coolants, first zero out the array */
		coolpr(ioQQQ,(char*)thermal.chClntLab[0],1,0.,"ZERO");

		/* push all coolants onto the stack */
		for( i=0; i < thermal.ncltot; i++ )
		{
			/* usually positive, although can be neg for coolants that heats, 
			 * only do positive here */
			ratio = (float)(thermal.cooling[i]/thermal.ctot);
			if( ratio >= EPS )
			{
				/*>>chng 99 jan 27, only cal when ratio is significant */
				coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], ratio,"DOIT");
			}
		}

		/* complete the printout for positive coolants */
		coolpr(ioQQQ,"DONE",1,0.,"DONE");

		/* now do cooling that was counted as a heat source if significant */
		if( thermal.heating[0][22]/thermal.ctot > 0.05 )
		{
			fprintf( ioQQQ, 
				"     All coolant heat greater than%6.2f%% of the total will be printed.\n", 
			  EPS*100. );

			coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
			for( i=0; i < thermal.ncltot; i++ )
			{
				ratio = (float)(thermal.heatnt[i]/thermal.ctot);
				if( fabs(ratio) >=EPS )
				{
					coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
					  ratio,"DOIT");
				}
			}
			coolpr(ioQQQ,"DONE",1,0.,"DONE");
		}
	}

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

	return;
}
#undef	PRT_DERIV

