/* 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 */
/* DynaEndIter called at end of iteration when advection is turned on */
/* DynaStartZone called at start of zone calculation when advection is turned on */
/* DynaEndZone called at end of zone calculation when advection is turned on */
/* DynaIonize, called from ionize to evaluate advective terms for current conditions */
/* DynaPresChngFactor, called from PressureChange to evaluate new density needed for
 * current conditions and wind solution, returns ratio of new to old density */
/* DynaZero zero some dynamics variables, called from zero.c */
/* DynaCreateArrays allocate some space needed to save the dynamics structure variables, 
 * called from DynaCreateArrays */
/* DynaPrtZone - called to print zone results */
/* DynaPunch punch info related to advection */
/* DynaPunch, punch output for dynamics solutions */
#include "cddefines.h"
#include "struc.h"
#include "radius.h"
#include "trace.h"
#include "converge.h"
#include "timesc.h"
#include "dense.h"
#include "co.h"
#include "thermal.h"
#include "pressure.h"
#include "phycon.h"
#include "rt.h"
#include "wind.h"
#include "hmi.h"
#include "dynamics.h"
static int ipUpstream=0,iphUpstream=0,ipyUpstream=0;

/* 
 * >>chng 01 mar 16, incorporate advection within dynamical solutions
 * this file contains the routines that incorporeate effects of dynamics and advection upon the
 * thermal and ionization solutions.  
 *
 * This code was originally developed in March 2001 during
 * a visit to UNAM Morelia, in collaboration with Robin Williams, Jane Arthur, and Will Henney.
 * Development was continued by email, and in a meeting in July/Aug 2001 at U Cardiff
 * Further development June 2002, UNAM Morelia (Cloudy and the Duendes Verdes)
 */

/* turn this on for dynamics printout */
#define DIAG_PRINT FALSE
#define MAINPRINT FALSE

/* the interpolated upstream densities of all ionization stages,
 * [element][ion] */
static double UpstreamIon[LIMELM][LIMELM+2];
/* total abundance of each element per hydrogen */
static double UpstreamElem[LIMELM];

/* hydrogen molecules */
static double UpstreamMolec[N_H_MOLEC];

/* routine called at end of iteration to determine new step sizes */
static void DynaNewStep(void);

/* routine called at end of iteration to save values in previous iteration */
static void DynaSaveLast(void);

/* routine called to determine mass flux at given distance */
/* static float DynaFlux(double depth); */

/* lookahead distance, as derived in DynaEndIter */
static double Dyn_dr;

/* advected work */
static double AdvecSpecificEnthalpy;

/* the upstream hydrogen atoms per unit vol*/
static double Upstream_hden;

/* HI ionization structure from previous iteration */
static float *Old_histr/*[NZLIM]*/ ,
	/* Lyman continuum optical depth from previous iteration */
	*Old_xLyman_depth/*[NZLIM]*/,
	/* depth of position in previous iteration */
	*Old_depth/*[NZLIM]*/,
	/* old n_p density from previous iteration */
	*Old_hiistr/*[NZLIM]*/,
	/* old pressure from previous iteration */
	*Old_pressure/*[NZLIM]*/,
	/* H density - particles per unit vol */
	*Old_hden/*[NZLIM]*/ ,
	/* density - total grams per unit vol */
	*Old_DenMass/*[NZLIM]*/ ,
	/* sum of enthalpy and kinetic energy per gram */
	*EnthalpyDensity/*[NZLIM]*/,
	/* old electron density from previous iteration */
	*Old_ednstr/*[NZLIM]*/,
	/* sum of enthalpy and kinetic energy per gram */
	*Old_EnthalpyDensity/*[NZLIM]*/;

static float **Old_Molec;

/* the ionization fractions from the previous iteration */
static float ***Old_xIonDense;

/* the gas phase abundances from the previous iteration */
static float **Old_gas_phase;

/* the number of zones that were saved in the previous iteration */
static long int nOld_zone;

/* the number of zones that were saved in the previous iteration */
static float DivergePresInteg;

/* ============================================================================== */
/* DynaPresChngFactor, called from PressureChange to evaluate factor needed
 * to find new density needed for 
 * current conditions and wind solution, returns ratio of new to old density,
 * called when wind velocity is negative */

/* object is to get the local ram pressure
 * RamNeeded = pressure.PresTotlInit + pressure.PresGasCurr + pressure.PresInteg;
 * to equal the sum of the inital pressur at the illuminated face, the local gas pressure,
 * and the integrate radiative acceleration from the incident continuum 
 *
 * the local gas pressure is linear in density if the temperature is constant,
 *
 * the local ram pressure is inversely linear in density because of the relationship
 * between velocity and density introduced by the mass flux conservation 
 */
#define SUBSONIC   1
#define SUPERSONIC 2
/*#define FREE       3*/
#define STRONGD    4
#define ORIGINAL   5
#define SHOCK      6
#define ANTISHOCK  7
#define ANTISHOCK2 8
double DynaPresChngFactor(void)
{
	double factor,
		er,
		width;
	static double rsave = -1.,
		dp = -1.,
		dpp = -1.,
		erp = -1.,
		erpp = -1.;
		/* TODO	2	this should be setable at run time as per Robin Williams email
		rshock = 4e16; */
	static int lastzone = -1,
		zonepresmode,
		globalpresmode;
	int ipPRE;

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

	/* following makes sure old wind v only updated once per zone,
	 * difference between old and new velocity is used to control zoning in nextdr */
	/*lint -e777 float equal */
	if( radius.depth != rsave )
	/*lint +e777 float equal */
	{
		rsave = radius.depth;
		/* evaluate radiative acceleration */
		RT_radiative_acceleration();
	}

	/* update the current pressure */
	PresTotCurrent();

	/* update the current desired pressure */
	pressure.PresTotlCorrect = pressure.PresTotlInit + pressure.PresInteg*pressure.lgContRadPresOn
		+ DivergePresInteg;

	if( DIAG_PRINT)
		fprintf(stdout,"Pressure: init %g +rad %g +diverge %g = %g cf %g\n",
			pressure.PresTotlInit, pressure.PresInteg*pressure.lgContRadPresOn,
			DivergePresInteg, pressure.PresTotlCorrect, pressure.PresTotlCurr);
	/*fprintf(ioQQQ,"DEBUG\t%.2f\thden\t%.4e\tPtot\t%.4e\tPgas\t%.4e\n",
		fnzone, dense.gas_phase[ipHYDROGEN],pressure.PresTotlCurr,pressure.PresGasCurr );*/

	/* dynamics.lgSetPresMode is flag to indicate sane value of dynamics.chPresMode.  
	 * If set true with SET DYNAMICS PRESSURE MODE
	 * then do not override that choice */
	if( !dynamics.lgSetPresMode )
	{
		/* above set true if pressure mode was set with
		 * SET DYNAMICS PRESSURE MODE - if we got here
		 * it was not set, and must make a guess */
		if (pressure.PresGasCurr < pressure.PresRamCurr)
			strcpy( dynamics.chPresMode , "supersonic" );
		else
			strcpy( dynamics.chPresMode , "subsonic" );
		/* clear the flag - pressure mode has been set */
		dynamics.lgSetPresMode = TRUE;
	}

	/* if globally looking for transonic solution, then locally sometimes
	 * supersonic, sometimes subsonic - this branch sets global flag,
	 * which can also be set with SET DYNAMICS PRESSURE MODE.
	 * Under default conditions, ChPresMode was set in previous branch
	 * to sub or super sonic depending on current velocity on first time*/
	if( strcmp( dynamics.chPresMode , "original" ) == 0 )
	{
		globalpresmode = ORIGINAL;
		pressure.lgSonicPointAbortOK = TRUE;
	}
	else if( strcmp( dynamics.chPresMode , "subsonic" ) == 0 )
	{
		globalpresmode = SUBSONIC;
		pressure.lgSonicPointAbortOK = TRUE;
	}
	/* supersonic */
	else if( strcmp( dynamics.chPresMode , "supersonic" ) == 0 )
	{
		globalpresmode = SUPERSONIC;
		pressure.lgSonicPointAbortOK = TRUE;
	}
	/* strong d */
	else if( strcmp( dynamics.chPresMode , "strongd" ) == 0 )
	{
		globalpresmode = STRONGD;
		pressure.lgSonicPointAbortOK = FALSE;
	}
	else if( strcmp( dynamics.chPresMode , "shock" ) == 0 )
	{
		globalpresmode = SHOCK;
		pressure.lgSonicPointAbortOK = FALSE;
	}
	else if( strcmp( dynamics.chPresMode , "antishock-by-mach" ) == 0 )
	{
		globalpresmode = ANTISHOCK2;
		pressure.lgSonicPointAbortOK = FALSE;
	}
	else if( strcmp( dynamics.chPresMode , "antishock" ) == 0 )
	{
		globalpresmode = ANTISHOCK;
		pressure.lgSonicPointAbortOK = FALSE;
	}

	/* this sets pressure mode for the current zone based on global mode
	 * and local conditions */
	if (globalpresmode == ORIGINAL)
	{
		/* in this mode use comparison between ram and gas pressure to determine
		 * whether sub or super sonic */
		if (pressure.PresGasCurr > pressure.PresRamCurr)
		{
			zonepresmode = SUBSONIC;
		}
		else
		{
			zonepresmode = SUPERSONIC;
		}
	}
	else if (globalpresmode == STRONGD)
	{
		if (nzone <= 1)
			zonepresmode = SUPERSONIC;
	}
	else if (globalpresmode == SUBSONIC)
	{
		zonepresmode = SUBSONIC;
	}
	else if (globalpresmode == SUPERSONIC)
	{
		zonepresmode = SUPERSONIC;
	}
	else if (globalpresmode == SHOCK)
	{
		if (radius.depth < dynamics.ShockDepth)
		{
			zonepresmode = SUBSONIC;
		}
		else
		{
			zonepresmode = SUPERSONIC;
		}
	}
	else if (globalpresmode == ANTISHOCK)
	{
		if (radius.depth < dynamics.ShockDepth)
		{
			zonepresmode = SUPERSONIC;
		}
		else
		{
			zonepresmode = SUBSONIC;
		}
	}
	else if (globalpresmode == ANTISHOCK2)
	{
		/* WJH 19 May 2004: This version of the antishock mode will
		  insert the antishock at the point where the isothermal Mach
		  number falls below a certain value, dynamics.ShockMach */
		if ( fabs(wind.windv) > dynamics.ShockMach*sqrt(pressure.PresGasCurr/dense.xMassDensity))
		{
			zonepresmode = SUPERSONIC;
		}
		else
		{
			zonepresmode = SUBSONIC;
		}
	}
	else
	{
		printf("Need to set global pressure mode\n");
		exit(-1);
	}
		
	er = pressure.PresTotlCurr-pressure.PresTotlCorrect;
	/* fprintf(ioQQQ,"Ds %ld: %ld; %g %g %g %g %g %g %g %g\n",iteration,nzone,dense.gas_phase[ipHYDROGEN],er,pressure.PresTotlCorrect,pressure.PresTotlCurr,phycon.te,thermal.ctot,thermal.htot,phycon.EnthalpyDensity); */
	/* fprintf(ioQQQ,"Ds %ld: %ld; %g %g %g %g\n",iteration,nzone,dense.gas_phase[ipHYDROGEN],er,pressure.PresTotlCurr,phycon.te); */
	if (globalpresmode == ORIGINAL || lastzone != nzone || fabs(er-erp) < SMALLFLOAT) 
	{
		/* First time in this zone, or when last step made no change, 
		 * take step hopefully in the right direction...
		 * ...or at least somewhere */
		if ( zonepresmode == SUBSONIC ) 
		{
			/* when gas pressure dominated need to increase density to increase pressure */
			factor = pressure.PresTotlCorrect / pressure.PresTotlCurr;
			ipPRE = 0;
		}
		else
		{
			/* when ram pressure dominated need to decrease density to increase pressure */
			factor = pressure.PresTotlCurr / pressure.PresTotlCorrect;
			ipPRE = 1;
		}
		if (fabs(factor-1.) > 0.01)
		{
			factor = 1.+sign(0.01,factor-1.);
		}
		erp = er;
		dp = dense.gas_phase[ipHYDROGEN];
		erpp = -1.;
		dpp = -1.;
	}
	else
	{
#if 0
		printf("Ds: %d; %g %g %g; %g %g %g tot %g\n",nzone,dense.gas_phase[ipHYDROGEN],dp,dpp,er,erp,erpp,
					 pressure.PresTotlCorrect);
#endif
		if(1 || dpp < 0. || fabs((dense.gas_phase[ipHYDROGEN]-dp)*(dp-dpp)*(dpp-dense.gas_phase[ipHYDROGEN])) < SMALLFLOAT) 
		{
			/* second iteration on this zone, do linerar fit to previous Pres - rho curve
			 * Linear approximation to guess root with two independent samples */
			factor = (dense.gas_phase[ipHYDROGEN]*erp-dp*er)/(erp-er)/dense.gas_phase[ipHYDROGEN];
			/* Limit step length to `reasonable' extrapolation */
			width = fabs(1.-dp/dense.gas_phase[ipHYDROGEN]);
			if (width > 1e-2)
				width = 1e-2;

			/* Subsonic case: pressure ^ with density ^ => increase density further */
			/* Super "  case: pressure ^ with density v => decrease density further */

			/* printf("Presmode %d flag %g factor %g\n",zonepresmode,(er-erp)*(dense.gas_phase[ipHYDROGEN]-dp),factor); */
			/* WJH 21 May 2004: I think that this part is to force the solution
			 * towards the required branch when it appears to have the "wrong" value
			 * of dP/drho */
			if (zonepresmode == SUBSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) < 0)
			{				
				factor = 1+3*width;
			}
			else if (zonepresmode == SUPERSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) > 0)
			{
				factor = 1-3*width;
			}

			if (fabs(factor-1.) > 3*width)
			{
				factor = 1.+sign(3*width,factor-1.);
			}
			ipPRE = 2;
			if (fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT) 
			{
				dpp = dp;
				erpp = erp;
				dp = dense.gas_phase[ipHYDROGEN];
				erp = er;
			}
		}
		else
		{
			/* 3rd or more solutions in this zone so have extensive data
			 * on pressure - density relation 
			 * Use quadratic fit to last three errors to estimate optimum */
			double a, b, c, q, dmin, emin, dsol, f1 , f2;
			int iroot;
			a = er/(dense.gas_phase[ipHYDROGEN]-dp)/(dense.gas_phase[ipHYDROGEN]-dpp) +
				erp/(dp-dense.gas_phase[ipHYDROGEN])/(dp-dpp)+
				erpp/(dpp-dense.gas_phase[ipHYDROGEN])/(dpp-dp);
			b = (erp-erpp)/(dp-dpp) - a * (dp+dpp);
			c = erp - dp*(a*dp+b);
			dmin = (-0.5*b/a);
			emin = (a*dmin+b)*dmin+c;
			if (a < 0) 
			{
				a = -a;
				b = -b;
				c = -c;
			}
#if 0
			fprintf(ioQQQ,"Check 1: %g %g\n",er,(a*dense.gas_phase[ipHYDROGEN]+b)*dense.gas_phase[ipHYDROGEN]+c);
			fprintf(ioQQQ,"Check 2: %g %g\n",erp,(a*dp+b)*dp+c);
			fprintf(ioQQQ,"Check 3: %g %g\n",erpp,(a*dpp+b)*dpp+c);
#endif
			q = b*b-4*a*c;
			if (q < 0) 
			{
				/* Imaginary root, search for local minimum */
				/* printf("No root at %d (%g cf %g) => %g\n",nzone,q,b*b,dmin); */
				factor = dmin/dense.gas_phase[ipHYDROGEN];

				/*TODO	2	Switch to supersonic when bad enough failure for STRONGD --
				 * need to improve this logic when the p(rho) data is cleaner
				 * */
				if (globalpresmode == STRONGD && -q > 1e-3*b*b)
				{
					zonepresmode = SUBSONIC;
				}
			} 
			else
			{
				/* Look for nearest root */
				/* if (zonepresmode == SUPERSONIC || (zonepresmode != SUBSONIC && (dense.gas_phase[ipHYDROGEN]-dmin) < 0)) */

				if (zonepresmode == SUPERSONIC)
					iroot = 1;
				else
					iroot = 0;

				if (emin > 0)
					iroot = ! iroot;

				if (iroot)
				{
					/* Look for smaller root */
					if (b > 0) 
					{
						dsol = -(b+sqrt(q))/(2*a);
					} 
					else 
					{
						dsol = 2*c/(-b+sqrt(q));
					}
				}
				else
				{
					if (b < 0)
					{
						dsol = (-b+sqrt(q))/(2*a);
					}
					else
					{
						dsol = -2*c/(b+sqrt(q));
					}
				}
				factor = dsol/dense.gas_phase[ipHYDROGEN];
#if 0
				fprintf(ioQQQ,"Check 4: %g %g %d %g %g\n",dsol,(a*dsol+b)*dsol+c,iroot,dmin,emin);
#endif
			}
			/* Limit step length */
			f1 = fabs(1.-dpp/dense.gas_phase[ipHYDROGEN]);
			f2 = fabs(1.- dp/dense.gas_phase[ipHYDROGEN]);
			/*width = MAX2(fabs(1.-dpp/dense.gas_phase[ipHYDROGEN]),fabs(1.-dp/dense.gas_phase[ipHYDROGEN]));*/
			width = MAX2(f1,f2);
			/* width = MAX2(width,1e-2); */
			if (fabs(factor-1.) > 3*width)
			{
				factor = 1.+sign(3*width,factor-1.);
			}
			ipPRE = 3;
			if (fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT) 
			{
				dpp = dp;
				erpp = erp;
				dp = dense.gas_phase[ipHYDROGEN];
				erp = er;
			}
		}
	}		

#if 0
	printf("Out: %d; %g; %g %g; %g %g\n",nzone,factor*dense.gas_phase[ipHYDROGEN],dp,dpp,erp,erpp);
#endif
	lastzone = nzone;

	if( DIAG_PRINT )
		fprintf(ioQQQ,"windv %li r %g v %g f %g\n",
			nzone,radius.depth,wind.windv,DynaFlux(radius.depth));

	/* convergence trace at this level */
	if( trace.lgTrConvg>=1  )
	{
		fprintf( ioQQQ, 
			" DynaPresChngFactor mode %s scale factor %.3f vel %.3e\n",
			dynamics.chPresMode , factor , wind.windv );
	}

	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			char chPRE[][4] = {"gas" , "ram", "sec", "par" };
			fprintf(ioQQQ,
				"pre %s\tfac\t%.5f\n",
				chPRE[ipPRE],
				factor
				);
		}
	}

	/* identically zero velocities cannot occur */
	ASSERT( wind.windv != 0. );

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

/* ============================================================================== */
/* DynaIonize, called from ConvBase to evaluate advective terms for current conditions,
 * calculates terms to add to ionization balance equation */
void DynaIonize(void)
{
	long int nelem, ion, mol;
	double  timestep;

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

	/* the time (s) needed for gas to move Dyn_dr  */
	/* >>chng 02 dec 11 rjrw increase timestep when beyond end of previous zone, to allow -> eqm */
	timestep = -Dyn_dr/wind.windv;
	/* printf("%d %g %g %g %g\n",nzone,radius.depth,Dyn_dr,radius.depth-Old_depth[nOld_zone-1],timestep); */

	ASSERT(nzone<struc.nzlim );
	if( nzone > 0 )
		EnthalpyDensity[nzone-1] = (float)phycon.EnthalpyDensity;

	/* do nothing on first iteration or when looking beyond previous iteration */
	if( iteration < 2 || radius.depth > dynamics.oldFullDepth)
	{
		/* first iteration, return zero */
		dynamics.Cool = 0.;
		dynamics.Heat = 0.;
		dynamics.dCooldT = 0.;
		dynamics.dHeatdT = 0.;
		
		dynamics.Rate = 0.;

		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
		{
			for( ion=0; ion<nelem+2 ; ++ion )
			{
				dynamics.Source[nelem][ion] = 0.;
			}
		}

		for (mol=0;mol<N_H_MOLEC;mol++)
			dynamics.Molec[mol] = 0.;

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

	if( MAINPRINT )
	{
		fprintf(ioQQQ,"workwork\t%li\t%.3e\t%.3e\t%.3e\n",
			nzone,
			phycon.EnthalpyDensity,
			0.5*POW2(wind.windv)*dense.xMassDensity ,
			5./2.*pressure.PresGasCurr
			); /**/
	}

	/* net cooling due to advection */
	/* >>chng 01 aug 01, removed hden from dilution, put back into here */
	/* >>chng 01 sep 15, do not update this variable the very first time this routine
	 * is called at the new zone. */
	dynamics.Cool = phycon.EnthalpyDensity/timestep*dynamics.lgCoolHeat;
	dynamics.Heat = AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/timestep*dynamics.lgCoolHeat;
	/* temp deriv of cooling minus heating */
	dynamics.dCooldT = 5./2.*pressure.PresGasCurr/phycon.te/timestep*dynamics.lgCoolHeat;
	dynamics.dHeatdT = 0.*dynamics.lgCoolHeat;

#if 0
		/* printf("DynaCool %g %g %g\n",	dynamics.Cool,	phycon.EnthalpyDensity/timestep,AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/timestep); */
		if (dynamics.Cool > 0) {
			dynamics.Heat = 0.;
			/* temp deriv of cooling minus heating */
			dynamics.dCooldT = 5./2.*pressure.PresGasCurr/phycon.te/timestep*dynamics.lgCoolHeat;
			dynamics.dHeatdT = 0.*dynamics.lgCoolHeat;
		} else {
			dynamics.Heat = -dynamics.Cool;
			dynamics.Cool = 0.;
			/* temp deriv of cooling minus heating */
			dynamics.dCooldT = 0.*dynamics.lgCoolHeat;
			dynamics.dHeatdT = -5./2.*pressure.PresGasCurr/phycon.te/timestep*dynamics.lgCoolHeat;
		}
#endif

#	if 0
	if( MAINPRINT || nzone>17 && iteration == 10)
	{
		fprintf(ioQQQ,
			"dynamics cool-heat\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t %.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
			nzone,
			phycon.te, 
			dynamics.CoolHeat,
			thermal.htot,
			phycon.EnthalpyDensity/timestep,
			AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/timestep,
			phycon.EnthalpyDensity,
			AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN],
			dense.gas_phase[ipHYDROGEN],
			timestep);
	}
#	endif

	/* second or greater iteration, have advective terms */
	/* this will evaluate advective terms for current physical conditions */

	/* the rate (s^-1) atoms drift in from upstream, a source term for the ground */

	/* dynamics.Hatom/timestep is the source (cm^-3 s^-1) of neutrals,
		 normalized to (s^-1) by the next higher ionization state as for all 
		 other recombination terms */

	/* dynamics.xIonDense[ipHYDROGEN][1]/timestep is the sink (cm^-3 s^-1) of 
		 ions, normalized to (s^-1) by the ionization state as for all other
		 ionization terms */

	dynamics.Rate = 1./timestep;

	for (mol=0;mol<N_H_MOLEC;mol++)
	{
		dynamics.Molec[mol] = UpstreamMolec[mol]* dense.gas_phase[ipHYDROGEN]*dynamics.Rate;
	}

	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
	{
		if( dense.lgElmtOn[nelem] )
		{
			if (fabs(UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN] -dense.gas_phase[nelem])/dense.gas_phase[nelem]>=1e-3) {
				fprintf(ioQQQ,"Conservation error: elem %li upstream %g abund %g\n",
								nelem,UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN],dense.gas_phase[nelem]);
			}
			/* ASSERT( fabs(UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN] -dense.gas_phase[nelem])/dense.gas_phase[nelem]<1e-3); */
			for( ion=0; ion<dense.IonLow[nelem]; ++ion )
			{
				dynamics.Source[nelem][ion] = 0.;
			}
			for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
			{
				/* Normalize to next higher state in current zone, except at first iteration
				where upstream version may be a better estimate (required for
				convergence in the small timestep limit) */

				dynamics.Source[nelem][ion] = 
					/* UpstreamIon is ion number per unit hydrogen because dilution is 1/hden 
					 * this forms the ratio of upstream atom over current ion, per timestep,
					 * so Recomb has units s-1 */
					UpstreamIon[nelem][ion] * dense.gas_phase[ipHYDROGEN] / timestep;

			}
			for( ion=dense.IonHigh[nelem]+1;ion<nelem+2; ++ion )
			{
				dynamics.Source[nelem][ion] = 0.;
			}
		}
	}
#	if 0
	fprintf(ioQQQ,"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
		nzone,
		dynamics.Rate,
		dynamics.Source[ipHYDROGEN][0],
		dynamics.Rate,
		dynamics.Source[ipCARBON][3]);
#	endif
#if 0
	nelem = ipCARBON;
	ion = 3;
	/*if( nzone > 160 && iteration > 1 )*/
	fprintf(ioQQQ,"dynaionizeeee\t%li\t%i\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
		nzone, 
		ipUpstream,
		radius.depth ,
		Old_depth[ipUpstream],
		dense.xIonDense[nelem][ion], 
		UpstreamIon[nelem][ion]* dense.gas_phase[ipHYDROGEN],
		Old_xIonDense[ipUpstream][nelem][ion] ,
		dense.xIonDense[nelem][ion+1], 
		UpstreamIon[nelem][ion+1]* dense.gas_phase[ipHYDROGEN],
		Old_xIonDense[ipUpstream][nelem][ion+1] ,
		timestep,
		dynamics.Source[nelem][ion]
		);
#endif
	if( MAINPRINT )
	{
		fprintf(ioQQQ,"    DynaIonize, %4li photo=%.2e , H recom= %.2e \n",
			nzone,dynamics.Rate , dynamics.Source[0][0]  );
	}

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

/* ============================================================================== */
/* DynaStartZone called at start of zone calculation when advection is turned on */
void DynaStartZone(void)
{
	/* this routine is called at the start of a zone calculation, by ZoneStart:
	 *
	 * it sets deduced variables to zero if first iteration,
	 *
	 * if upstream depth is is outside the computed structure on previous iteration,
	 * return value at shielded face 
	 *
	 * Also calculates discretization_error, an estimate of the accuracy of the source terms.
	 *
	 * */

	/* this is index of upstream cell in structure stored from previous iteration */
	double upstream, dilution, dilutionleft, dilutionright, nextfrac;

	/* Properties for cell half as far upstream, used to converge timestep */
	double hupstream, hnextfrac=-BIGFLOAT, hion, hmol/*, hhden*/;

	/* Properties for cell at present depth, used to converge timestep */
	double ynextfrac=-BIGFLOAT, yion, ymol/*, yhden*/;

	long int nelem , ion, mol;

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

	/* do nothing on first iteration */
	if( iteration < 2 )
	{
		Upstream_hden = 0.;
		AdvecSpecificEnthalpy = 0.;
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
		{
			for( ion=0; ion<nelem+2 ; ++ion )
			{
				UpstreamIon[nelem][ion] = 0.;
			}
		}
		/* hydrogen molecules */
#if 1
		for (mol=0; mol<N_H_MOLEC;mol++) 
		{
			UpstreamMolec[mol] = 0;
		}
#else
		Upstream_htwo = 0;
		Upstream_htwo_star = 0;
		Upstream_hehp = 0;
		Upstream_h2plus = 0;
		Upstream_h3plus = 0;
#endif

		ipUpstream = 0;
		iphUpstream = 0;
		ipyUpstream = 0;

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

	/* radius.depth is distance from illuminated face of cloud to outer edge of
	 * current zone, which has thickness radius.drad */

	/* find where the down stream point is, in previous iteration,
	 * NB neg sign since velocity is negative, we are looking for point
	 * where current gas cell was, in previous iteration */

	/* true, both advection and wind solution */
	/* don't interpolate to the illuminated side of the first cell */
	upstream = MAX2(Old_depth[0] , radius.depth + Dyn_dr);
	hupstream = 0.5*(radius.depth + upstream);

	/* now locate upstream point in previous stored structure,
	 * will be at the same point or higher than we found previously */
	while( (Old_depth[ipUpstream+1] < upstream ) && 
		( ipUpstream < nOld_zone-1 ) )
	{
		++ipUpstream;
	}
	ASSERT( ipUpstream <= nOld_zone-1 );
	
	/* above loop will give ipUpstream == nOld_zone-1 if computed structure has been overrun */

	if( (ipUpstream != nOld_zone-1)&& ((Old_depth[ipUpstream+1] - Old_depth[ipUpstream])> SMALLFLOAT) )
	{
		/* we have not overrun radius scale of previous iteration */
		nextfrac = ( upstream - Old_depth[ipUpstream])/
			(Old_depth[ipUpstream+1] - Old_depth[ipUpstream]);
		Upstream_hden = Old_hden[ipUpstream] +
			(Old_hden[ipUpstream+1] - Old_hden[ipUpstream])*
			nextfrac;
		dilutionleft = 1./Old_hden[ipUpstream];
		dilutionright = 1./Old_hden[ipUpstream+1];

		/* fractional changes in density from passive advection */
		/* >>chng 01 May 02 rjrw: use hden for dilution */
		/* >>chng 01 aug 01, remove hden here, put back into vars when used in DynaIonize */
		dilution = 1./Upstream_hden;

		/* the advected enthalphy */
		AdvecSpecificEnthalpy = (Old_EnthalpyDensity[ipUpstream]*dilutionleft +
			(Old_EnthalpyDensity[ipUpstream+1]*dilutionright - Old_EnthalpyDensity[ipUpstream]*dilutionleft)*
			nextfrac);

		ASSERT( Old_depth[ipUpstream] <= upstream && upstream <= Old_depth[ipUpstream+1] );
		ASSERT( (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) *
			(AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright) > -SMALLFLOAT );
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
		{
			UpstreamElem[nelem] = 0.;
			for( ion=0; ion<nelem+2 ; ++ion )
			{
				/* Robin - I made several changes like the following - it seems easier to
				 * bring out the division by the old hydrogen density rather than putting in
				 * dilution and then looking for how diluation is defined.  I think the code
				 * is equivalent */
				/* UpstreamIon is ion number per unit hydrogen, both at the upstream position */
				UpstreamIon[nelem][ion] = 
					Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream] +
					(Old_xIonDense[ipUpstream+1][nelem][ion]/Old_hden[ipUpstream+1] - 
					Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream])*
					nextfrac;

				UpstreamElem[nelem] += UpstreamIon[nelem][ion];
			}
		}
#if 1
		for (mol=0;mol<N_H_MOLEC;mol++)
		{
			UpstreamMolec[mol] = Old_Molec[ipUpstream][mol]/Old_hden[ipUpstream] +
				(Old_Molec[ipUpstream+1][mol]/Old_hden[ipUpstream+1] - 
				 Old_Molec[ipUpstream][mol]/Old_hden[ipUpstream])*
				nextfrac;
			if (mol > 1)
				UpstreamElem[ipHYDROGEN] += UpstreamMolec[mol]*hmi.nProton[mol];
		}
#else
		Upstream_hminus = 
			Old_hminus[ipUpstream]/Old_hden[ipUpstream] +
			(Old_hminus[ipUpstream+1]/Old_hden[ipUpstream+1] - 
			Old_hminus[ipUpstream]/Old_hden[ipUpstream])*
			nextfrac;
		Upstream_htwo = 
			Old_htwo[ipUpstream]/Old_hden[ipUpstream] +
			(Old_htwo[ipUpstream+1]/Old_hden[ipUpstream+1] - 
			Old_htwo[ipUpstream]/Old_hden[ipUpstream])*
			nextfrac;
		Upstream_htwo_star = 
			Old_htwo_star[ipUpstream]/Old_hden[ipUpstream] +
			(Old_htwo_star[ipUpstream+1]/Old_hden[ipUpstream+1] - 
			Old_htwo_star[ipUpstream]/Old_hden[ipUpstream])*
			nextfrac;
		Upstream_hehp = 
			Old_hehp[ipUpstream]/Old_hden[ipUpstream] +
			(Old_hehp[ipUpstream+1]/Old_hden[ipUpstream+1] - 
			Old_hehp[ipUpstream]/Old_hden[ipUpstream])*
			nextfrac;
		Upstream_h2plus = 
			Old_h2plus[ipUpstream]/Old_hden[ipUpstream] +
			(Old_h2plus[ipUpstream+1]/Old_hden[ipUpstream] - 
			Old_h2plus[ipUpstream]/Old_hden[ipUpstream])*
			nextfrac;
		Upstream_h3plus = 
			Old_h3plus[ipUpstream]/Old_hden[ipUpstream] +
			(Old_h3plus[ipUpstream+1]/Old_hden[ipUpstream] - 
			Old_h3plus[ipUpstream]/Old_hden[ipUpstream])*
			nextfrac;
		UpstreamElem[ipHYDROGEN] += Upstream_hminus+2*Upstream_htwo+2*Upstream_htwo_star+Upstream_hehp+
			2*Upstream_h2plus+3*Upstream_h3plus;
#endif
	}
	else
	{
		/* SPECIAL CASE - we have overrun the previous iteration's radius */
		Upstream_hden = Old_hden[ipUpstream];
		/* fractional changes in density from passive advection */
		/* >>chng 01 aug 01, remove hden here, put back into vars when used in DynaIonize */
		dilution = 1./Upstream_hden;
		/*AdvecSpecificEnthalpy = Old_EnthalpyDensity[ipUpstream]*dilution;*/
		AdvecSpecificEnthalpy = Old_EnthalpyDensity[ipUpstream]/Old_hden[ipUpstream];
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
		{
			UpstreamElem[nelem] = 0.;
			for( ion=0; ion<nelem+2 ; ++ion )
			{
				/* UpstreamIon is ion number per unit hydrogen */
				UpstreamIon[nelem][ion] = 
					Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream];
				UpstreamElem[nelem] += UpstreamIon[nelem][ion];
			}
		}
#if 1
		for (mol=0;mol<N_H_MOLEC;mol++) 
		{
			UpstreamMolec[mol] = Old_Molec[ipUpstream][mol]/Old_hden[ipUpstream];
			if (mol > 1)
				UpstreamElem[ipHYDROGEN] += UpstreamMolec[mol]*hmi.nProton[mol];
		}
#else
		Upstream_hminus = 
			Old_hminus[ipUpstream]/Old_hden[ipUpstream];
		Upstream_htwo = 
			Old_htwo[ipUpstream]/Old_hden[ipUpstream];
		Upstream_htwo_star = 
			Old_htwo_star[ipUpstream]/Old_hden[ipUpstream];
		Upstream_hehp = 
			Old_hehp[ipUpstream]/Old_hden[ipUpstream];
		Upstream_h2plus = 
			Old_h2plus[ipUpstream]/Old_hden[ipUpstream];
		Upstream_h3plus = 
			Old_h3plus[ipUpstream]/Old_hden[ipUpstream];
		UpstreamElem[ipHYDROGEN] += Upstream_hminus+2*Upstream_htwo+2*Upstream_htwo_star+Upstream_hehp+
			2*Upstream_h2plus+3*Upstream_h3plus;
#endif
	}

	/* Repeat enough of the above for half-step and no-step to judge convergence:
		 the result of this code is the increment of discretization_error */

	while( (Old_depth[iphUpstream+1] < hupstream ) && 
		( iphUpstream < nOld_zone-1 ) )
	{
		++iphUpstream;
	}
	ASSERT( iphUpstream <= nOld_zone-1 );

	while( (Old_depth[ipyUpstream+1] < radius.depth ) && 
		( ipyUpstream < nOld_zone-1 ) )
	{
		++ipyUpstream;
	}
	ASSERT( ipyUpstream <= nOld_zone-1 );

	dynamics.dRad = BIGFLOAT;

	if (iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
	{
		hnextfrac = ( hupstream - Old_depth[iphUpstream])/
			(Old_depth[iphUpstream+1] - Old_depth[iphUpstream]);
		/* >>chng 02 nov 11, hhden never appears on rhs of = */
		/*hhden = Old_hden[iphUpstream] +
			(Old_hden[iphUpstream+1] - Old_hden[iphUpstream])*
			hnextfrac;*/
	}
	else
	{
		/* >>chng 02 nov 11, hhden never appears on rhs of = */
		/*hhden = Old_hden[iphUpstream];*/
	}
	if (ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
	{
		ynextfrac = ( radius.depth - Old_depth[ipyUpstream])/
			(Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]);
		/* >>chng 02 nov 11, yhden never appears on rhsof = */
		/*yhden = Old_hden[ipyUpstream] +
			(Old_hden[ipyUpstream+1] - Old_hden[ipyUpstream])*
			ynextfrac;*/
	}
	else
	{
		/* >>chng 02 nov 11, yhden never appears on rhsof = */
		/*yhden = Old_hden[ipyUpstream];*/
	}

	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
	{
		for( ion=0; ion<nelem+2 ; ++ion )
		{
			double f1;
			if (ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
			{
				yion = 
					Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream] +
					(Old_xIonDense[ipyUpstream+1][nelem][ion]/Old_hden[ipyUpstream+1] - 
					 Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream])*
					ynextfrac;
			}
			else
			{		
				yion = Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream];
			}
			if (iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
			{
				hion = 
					Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream] +
					(Old_xIonDense[iphUpstream+1][nelem][ion]/Old_hden[iphUpstream+1] - 
					 Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream])*
					hnextfrac;
			}
			else
			{		
				hion = Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream];
			}

			/* the proposed thickness of the next zone, there will be a scale factor, something
			 * like 1/500, that will be applied when this is used in nextdr */
			f1 = fabs(yion - UpstreamIon[nelem][ion] );
			dynamics.dRad = MIN2(dynamics.dRad,Dyn_dr * 
			/* don't pay attention to species with abundance relative to H less tghan 1e-6 */
				MAX2(yion + UpstreamIon[nelem][ion],1e-6 ) / 
				SDIV(f1));
				/*MAX2(fabs(yion - UpstreamIon[nelem][ion] ),SMALLFLOAT));*/
			/* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */

			/* Must be consistent with convergence_error below */
			/* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */
			/* >>chngf 02 aug 01, multiply by cell width */
			/* this error is error due to the advection length not being zero - a finite
			 * advection length.  no need to bring convergence error to below
			 * discretization error.  when convergece error is lower than a fraction of this
			 * errror we reduce the advection length. */
			dynamics.discretization_error += POW2(yion-2.*hion+UpstreamIon[nelem][ion]); 
			dynamics.error_scale2 += POW2(UpstreamIon[nelem][ion]-yion);
		}
	}

	for (mol=0; mol < N_H_MOLEC; mol++) 
	{
		double f1;
		if (ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
		{
			ymol = 
				Old_Molec[ipyUpstream][mol]/Old_hden[ipyUpstream] +
				(Old_Molec[ipyUpstream+1][mol]/Old_hden[ipyUpstream+1] - 
				 Old_Molec[ipyUpstream][mol]/Old_hden[ipyUpstream])*
				ynextfrac;
		}
		else
		{		
			ymol = Old_Molec[ipyUpstream][mol]/Old_hden[ipyUpstream];
		}
		if (iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
		{
			hmol = 
				Old_Molec[iphUpstream][mol]/Old_hden[iphUpstream] +
				(Old_Molec[iphUpstream+1][mol]/Old_hden[iphUpstream+1] - 
				 Old_Molec[iphUpstream][mol]/Old_hden[iphUpstream])*
				hnextfrac;
		}
		else
		{		
			hmol = Old_Molec[iphUpstream][mol]/Old_hden[iphUpstream];
		}

		/* the proposed thickness of the next zone, there will be a scale factor, something
		 * like 1/500, that will be applied when this is used in nextdr */
		f1 = fabs(ymol - UpstreamMolec[mol] );
		dynamics.dRad = MIN2(dynamics.dRad,Dyn_dr * 
			/* don't pay attention to species with abundance relative to H less tghan 1e-6 */
			MAX2(ymol + UpstreamMolec[mol],1e-6 ) / 
			SDIV(f1));
			/*MAX2(fabs(ymol - UpstreamMolec[mol] ),SMALLFLOAT));*/
		/* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */

		/* Must be consistent with convergence_error below */
		/* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */
		/* >>chngf 02 aug 01, multiply by cell width */
		dynamics.discretization_error += POW2(ymol-2.*hmol+UpstreamMolec[mol]); 
		dynamics.error_scale2 += POW2(UpstreamMolec[mol]-ymol);
	}

	if( MAINPRINT )
	{
		fprintf(ioQQQ," DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n",
			nzone,dynamics.Rate , dynamics.Source[0][0] , dilution*dense.gas_phase[ipHYDROGEN] );
	}

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

/* DynaEndZone called at end of zone calculation when advection is turned on */
void DynaEndZone(void)
{
#	ifdef DEBUG_FUN
	fputs( "<+>DynaEndZone()\n", debug_fp );
#	endif

	/* this routine is called at the end of a zone calculation, by ZoneEnd */

	DivergePresInteg += wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad)); 
	if(DIAG_PRINT)
		fprintf(ioQQQ,"Check dp: %g %g mom %g %g mas %g\n",
			wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad)),
			2*wind.windv*DynaFlux(radius.depth)*radius.drad/(1e16-radius.depth),
			wind.windv*DynaFlux(radius.depth),
			wind.windv*DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth),
			DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth));

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


/* ============================================================================== */
/* DynaEndIter called at end of iteration when advection is turned on */
void DynaEndIter(void)
{
	/* this routine is called by IterRestart at the end of an iteration 
	 * when advection is turned on.  currently it only derives a 
	 * timestep by looking at the spatial derivative of
	 * some stored quantities */

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

	DivergePresInteg = 0.;

	/* This routine is only called if advection is turned on at the end of
	 * each iteration.  The test 
	 * is to check whether wind velocity is also set by dynamics code */

	if (iteration == 2)
	{
	
		if( dynamics.AdvecLengthInit> 0. )
		{
			Dyn_dr =  dynamics.AdvecLengthInit;
		}
		else
		{
			/* -ve dynamics.adveclengthlimit sets length as fraction of first iter total depth */
			Dyn_dr = -dynamics.AdvecLengthInit*radius.depth;
		}

		if( MAINPRINT )
		{
			fprintf(ioQQQ," DynaEndIter, dr=%.2e \n",
				Dyn_dr );
		}
	}
	else
	{
		DynaNewStep();
	}

	ASSERT( Dyn_dr > 0. );

	/* reset the upstream counters */
	ipUpstream = iphUpstream = ipyUpstream = 0;
	dynamics.discretization_error = 0.;
	dynamics.error_scale2 = 0.;

	DynaSaveLast();

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

static void DynaNewStep(void)
{
	long int ilast = 0,
		i,
		nelem,
		ion,
		mol;

	double nextfrac=-BIGFLOAT,
		Oldi_hden,
		Oldi_ion,
		Oldi_mol;

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

	/*n = MIN2(nzone, NZLIM-1);*/
	dynamics.convergence_error = 0;
	dynamics.error_scale1 = 0.;

	ASSERT( nzone < struc.nzlim);
	for (i=0;i<nzone;++i) 
	{
		/* Interpolate for present position in previous solution */
		while( (Old_depth[ilast] < struc.depth[i] ) && 
			( ilast < nOld_zone-1 ) )
		{
			++ilast;
		}
		ASSERT( ilast <= nOld_zone-1 );
			
		if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
		{
			nextfrac = ( struc.depth[i] - Old_depth[ilast])/
				(Old_depth[ilast+1] - Old_depth[ilast]);
			Oldi_hden = Old_hden[ilast] +
				(Old_hden[ilast+1] - Old_hden[ilast])*
				nextfrac;
		}
		else
		{
			Oldi_hden = Old_hden[ilast];
		}
		/* Must be consistent with discretization_error above */
		/* >>chngf 02 aug 01, multiply by cell width */
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
		{
			for( ion=0; ion<nelem+2 ; ++ion )
			{
				if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
				{
					Oldi_ion = (Old_xIonDense[ilast][nelem][ion] +
						(Old_xIonDense[ilast+1][nelem][ion]-Old_xIonDense[ilast][nelem][ion])*
						nextfrac);
				}
				else
				{
					Oldi_ion = Old_xIonDense[ilast][nelem][ion];
				}
				dynamics.convergence_error += POW2(Oldi_ion/Oldi_hden-struc.xIonDense[nelem][ion][i]/struc.hden[i]) /* *struc.dr[i] */ ;  

				/* >>chng 02 nov 11, add first error scale estimate from Robin */
				dynamics.error_scale1 += POW2(struc.xIonDense[nelem][ion][i]/struc.hden[i]);			
			}
		}
		for( mol=0; mol < N_H_MOLEC; mol++)
		{
			if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
			{
				Oldi_mol = (Old_Molec[ilast][mol] +
										(Old_Molec[ilast+1][mol]-Old_Molec[ilast][mol])*
										nextfrac);
			}
			else
			{
				Oldi_mol = Old_Molec[ilast][mol];
			}
			dynamics.convergence_error += POW2(Oldi_mol/Oldi_hden-struc.Hmolec[mol][i]/struc.hden[i]) /* *struc.dr[i] */ ;  
			
			/* >>chng 02 nov 11, add first error scale estimate from Robin
			 * used to normalize the above convergence_error */
			dynamics.error_scale1 += POW2(struc.Hmolec[mol][i]/struc.hden[i]);			
		}
	}

	/* convergence_error is an estimate of the convergence of the solution from its change during the last iteration,
		 discretization_error is an estimate of the accuracy of the advective terms, calculated in DynaStartZone above:
		 if dominant error is from the advective terms, need to make them more accurate.
	*/

	/* report properties of previous iteration */
	fprintf(ioQQQ,"DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n",
		Dyn_dr, dynamics.convergence_error , dynamics.discretization_error ,
		dynamics.error_scale1 , dynamics.error_scale2 
		);

	/* >>chng 02 nov 29, dynamics.convergence_tolerance is now set to 0.1 in init routine */
	if ( dynamics.convergence_error < dynamics.convergence_tolerance*dynamics.discretization_error ) 
		Dyn_dr /= 1.5;

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

static void DynaSaveLast(void)
{
	long int i,
		ion,
		nelem,
		mol;

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

	/* Save results from previous iteration */
	nOld_zone = nzone;
	dynamics.oldFullDepth = struc.depth[nzone-1];
	ASSERT( nzone < struc.nzlim );
	for( i=0; i<nzone ; ++i )
	{
		Old_histr[i] = struc.histr[i];
		Old_depth[i] = struc.depth[i];
		Old_xLyman_depth[i] = struc.xLyman_depth[i];
		/* old n_p density from previous iteration */
		Old_hiistr[i] = struc.hiistr[i];
		/* old pressure from previous iteration */
		Old_pressure[i] = struc.pressure[i];
		/* old electron density from previous iteration */
		Old_ednstr[i] = struc.ednstr[i];
		/* energy term */
		Old_EnthalpyDensity[i] = EnthalpyDensity[i];
		/* >>chng 02 May 2001 rjrw: add hden for dilution */
		Old_hden[i] = struc.hden[i];
		Old_DenMass[i] = struc.DenMass[i];
#if 1
		for (mol=0;mol<N_H_MOLEC;mol++)
		{
			Old_Molec[i][mol] = struc.Hmolec[mol][i];
		}
#else
		Old_hminus[i] = struc.hminus[i];
		Old_htwo[i] = struc.htwo[i];
		Old_htwo_star[i] = struc.htwo_star[i];
		Old_hehp[i] = struc.hehp[i];
		Old_h2plus[i] = struc.h2plus[i];
		Old_h3plus[i] = struc.h3plus[i];
#endif

		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
		{
			Old_gas_phase[i][nelem] = dense.gas_phase[nelem];
			for( ion=0; ion<nelem+2 ; ++ion )
			{
				Old_xIonDense[i][nelem][ion] = struc.xIonDense[nelem][ion][i];
			}
		}
	}
#	ifdef DEBUG_FUN
	fputs( " <->DynaSaveLast()\n", debug_fp );
#	endif
	return;
}

float DynaFlux(double depth)

{
	float flux;

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

	if (dynamics.FluxIndex == 0) 
	{
		flux = (float)dynamics.FluxScale;
	}	
	else 
	{
		flux = (float)(dynamics.FluxScale*pow(fabs(depth-dynamics.FluxCenter),dynamics.FluxIndex));
		if (depth < dynamics.FluxCenter)
			flux = -flux;
	}
	if (dynamics.lgFluxDScale) 
	{
		/*flux *= struc.DenMass[0]; */
		/* WJH 21 may 04, changed to use dense.xMassDensity0, which should be strictly constant */
		flux *= dense.xMassDensity0; 
	}

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

	return flux;
}

/* ============================================================================== */
/*DynaZero zero some dynamics variables, called from zero.c, 
 * before parsing commands */
void DynaZero( void )
{
	int ipISO;

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

	/* the number of zones in the previous iteration */
	nOld_zone = 0;

	/* by default advection is turned off */
	dynamics.lgAdvection = FALSE;
	/*dynamics.Velocity = 0.;*/
	AdvecSpecificEnthalpy = 0.;
	dynamics.Cool = 0.;
	dynamics.Heat = 0.;
	dynamics.dCooldT = 0.;
	dynamics.dHeatdT = 0.;
	dynamics.HeatMax = 0.;
	dynamics.CoolMax = 0.;
	dynamics.Rate = 0.;

	/* set inital value of the advection length,
	 * neg => fraction of depth of init model, + length cm */
	dynamics.AdvecLengthInit = -0.1;

	/* this is a tolerance for determining whether dynamics has converged */
	dynamics.convergence_tolerance = 0.1;

	/* this says that set dynamics pressure mode was set */
	dynamics.lgSetPresMode = FALSE;

	/* set default values for uniform mass flux */
	dynamics.FluxScale = 0.;
	dynamics.lgFluxDScale = FALSE;
	dynamics.FluxCenter = 0.;
	dynamics.FluxIndex = 0.;
	dynamics.dRad = BIGFLOAT;

	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* factor to turn off advection for H-like iso seq, no advection h-like*/
		dynamics.lgISO[ipISO] = 1;
	}
	/* factor to turn off advection for rest of ions, treated in bidiag, no advection h-like*/
	dynamics.lgMETALS = 1;
	dynamics.lgCoolHeat = 1;
	DivergePresInteg = 0.;

	dynamics.discretization_error = 0.;
	dynamics.error_scale2 = 0.;

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


/* ============================================================================== */
/* DynaCreateArrays allocate some space needed to save the dynamics structure variables, 
 * called from DynaCreateArrays */
void DynaCreateArrays( void )
{
	long int nelem,
		ns,
		i,
		ion,
		mol;

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

	if( (dynamics.Source = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ))) == NULL )
		BadMalloc();
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		if( (dynamics.Source[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ))) == NULL )
			BadMalloc();
		for( ion=0; ion<nelem+2; ++ion )
		{
			dynamics.Source[nelem][ion] = 0.;
		}
	}
	dynamics.Rate = 0.;

	if( (Old_EnthalpyDensity = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_ednstr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (EnthalpyDensity = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_DenMass = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_hden = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_pressure = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_histr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_hiistr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_depth = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_xLyman_depth = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )))) == NULL )
		BadMalloc();

	if( (Old_xIonDense = (float ***)MALLOC(sizeof(float **)*(unsigned)(struc.nzlim) ) )==NULL )
		BadMalloc();

	if( (Old_gas_phase = (float **)MALLOC(sizeof(float *)*(unsigned)(struc.nzlim) ) )==NULL )
		BadMalloc();

	if( (Old_Molec = (float **)MALLOC(sizeof(float *)*(unsigned)(struc.nzlim) ) )==NULL )
		BadMalloc();

	/* now create diagonal of space for ionization arrays */
	for( ns=0; ns < struc.nzlim; ++ns )
	{
		if( (Old_xIonDense[ns] = 
			(float**)MALLOC(sizeof(float*)*(unsigned)(LIMELM+3) ))==NULL )
			BadMalloc();

		if( (Old_gas_phase[ns] = 
			(float*)MALLOC(sizeof(float*)*(unsigned)(LIMELM+3) ))==NULL )
			BadMalloc();

		if( (Old_Molec[ns] = 
 		  (float*)MALLOC(sizeof(float*)*(unsigned)(N_H_MOLEC) ))==NULL )
			BadMalloc();

		for( nelem=0; nelem< (LIMELM+3) ;++nelem )
		{
			if( (Old_xIonDense[ns][nelem] = 
				(float*)MALLOC(sizeof(float)*(unsigned)(LIMELM+1) ))==NULL )
				BadMalloc();
		}
	}

	for( i=0; i < struc.nzlim; i++ )
	{
		/* these are values if H0 and tau_912 from previous iteration */
		Old_histr[i] = 0.;
		Old_xLyman_depth[i] = 0.;
		Old_depth[i] = 0.;
		dynamics.oldFullDepth = 0.;
		/* old n_p density from previous iteration */
		Old_hiistr[i] = 0.;
		/* old pressure from previous iteration */
		Old_pressure[i] = 0.;
		/* old electron density from previous iteration */
		Old_ednstr[i] = 0.;
		Old_hden[i] = 0.;
		Old_DenMass[i] = 0.;
		Old_EnthalpyDensity[i] = 0.;
		for( nelem=0; nelem< (LIMELM+3) ;++nelem )
		{
			for( ion=0; ion<LIMELM+1; ++ion )
			{
				Old_xIonDense[i][nelem][ion] = 0.;
			}
		}
		for (mol=0;mol<N_H_MOLEC;mol++)
		{
			Old_Molec[i][mol] = 0.;
		}
	}

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

/* ============================================================================== */
/* ParseDynaWind parse the wind command, called from ParseCommands */
void ParseDynaWind( char *chCard )
{
	long int i;
	int iVelocity_Type;
	int lgEOL;
	/* compiler flagged possible paths where dfdr used but not set -
	 * this is for safety/keep it happy */
	double dfdr=-BIGDOUBLE;

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

	/* Flag for type of velocity law: 
	 * 1 is original, give initial velocity at illuminated face
	 * 2 is face flux gradient (useful if face velocity is zero),
	 * set to zero, but will be reset if vel specified */
	iVelocity_Type = 0;
	/* wind structure, parameters are initial velocity and optional mass
	 * v in km/sec, mass in solar masses */
	if( (i = lgMatch("VELO",chCard))>0 )
	{
		i += 5;
		wind.windv0 = (float)(FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)*1e5);
		wind.windv = wind.windv0;
		iVelocity_Type = 1;
	}

	if( (i = lgMatch("DFDR",chCard))>0 )
	{
		/* velocity not specified, rather mass flux gradient */
		i += 5;
		dfdr = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		iVelocity_Type = 2;
	}

	/* center option, gives xxx */
	if( (i = lgMatch("CENT",chCard))>0 )
	{
		/* physical length in cm, can be either sign */
		i += 5;
		dynamics.FluxCenter = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	}

	/* flux index */
	if( (i = lgMatch("INDE",chCard))>0 )
	{
		/* power law index */
		i += 5;
		dynamics.FluxIndex = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	}

	/* the case where velocity was set */
	if (iVelocity_Type == 1) 
	{
		/* was flux index also set? */
		if (dynamics.FluxIndex == 0)
		{
			dynamics.FluxScale = wind.windv0;
			dynamics.lgFluxDScale = TRUE;
			/* Center doesn't mean much in this case -- make sure it's
			 * in front of grid so DynaFlux doesn't swap signs where
			 * it shouldn't */
			dynamics.FluxCenter = -1.;
		}
		else
		{
			/*TODO	2	Need to include divergence terms in pressure balance 
			 * if flux index is != 0 */
			/* velocity was set but flux index was not set - estimate it */
			dynamics.FluxScale = wind.windv0*
				pow(fabs(dynamics.FluxCenter),-dynamics.FluxIndex);

			dynamics.lgFluxDScale = TRUE;
			if (dynamics.FluxCenter > 0)
			{
				dynamics.FluxScale = -dynamics.FluxScale;
			}
		}
	} 
	/* the case where flux gradient is set */
	else if (iVelocity_Type == 2)
	{
		if (dynamics.FluxIndex == 0)
		{
			fprintf(ioQQQ,"Can't specify gradient when flux is constant!\n");
			/* use this exit handler, which closes out MPI when multiprocessing */
			cdEXIT(EXIT_FAILURE);
		}
		/*TODO	2	Need to include divergence terms in pressure balance 
		 * if flux index is != 0 */
		/* Can't specify FluxScale from dvdr rather than dfdr, as
		 * d(rho)/dr != 0 */ 
		dynamics.FluxScale = dfdr/dynamics.FluxIndex*
			pow(fabs(dynamics.FluxCenter),1.-dynamics.FluxIndex);
		if (dynamics.FluxCenter > 0)
		{
			dynamics.FluxScale = -dynamics.FluxScale;
		}
		dynamics.lgFluxDScale = FALSE;

		/* put in bogus value simply as flag -- assume that surface velocity
		 * is small or we wouldn't be using this to specify. */
		wind.windv0 = -0.01f;
	}
	else
	{
		/* assume old-style velocity-only specification */ 
		/* wind structure, parameters are initial velocity and optional mass
		 *  v in km/sec, mass in solar masses */
		i = 5;
		wind.windv0 = (float)(FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)*1e5);
		dynamics.FluxScale = wind.windv0;
		dynamics.FluxIndex = 0.;
		dynamics.lgFluxDScale = TRUE;
		/* Center doesn't mean much in this case -- make sure it's
		 * in front of grid so DynaFlux doesn't swap signs where
		 * it shouldn't */
		dynamics.FluxCenter = -1.;
		if( lgEOL ) 
		{
			NoNumb(chCard);
		}
	}
	
	wind.windv = wind.windv0;

#ifdef FOO
	fprintf(ioQQQ,"Scale %g (*%c) Index %g Center %g\n",
					dynamics.FluxScale,(dynamics.lgFluxDScale)?'D':'1',dynamics.FluxIndex,dynamics.FluxCenter);
#endif

	/* wind < 0 is flow from blister */
	if( wind.windv0 < 0. )
	{
		/* gravity is not included in hydro flow*/
		wind.comass = 0.;

		/* option to include advection */
		if( lgMatch( "ADVE" , chCard ) )
		{
			/* turn on advection */
			dynamics.lgAdvection = TRUE;

			/* turn off prediction of next zone's temperature, as guessed in ZoneStart,
			 * also set with no tepredictor */
			thermal.lgPredNextTe = FALSE;

			/* turn off both CO and H2 networks since advection not included, also say physical
			 * conditions are not ok*/
			/* the heavy element molecules are turned off for the moment,
			 * until I am confident that the H-only models work robustly */
			/* co.lgNoH2Mole = TRUE; */
			/*TODO	1	get following to work */
			co.lgNoCOMole = TRUE; /* >>chng 04 apr 23, tried to rm this line - problems */
			phycon.lgPhysOK = FALSE;/* >>chng 04 apr 23, rm this line */

			/* use the new temperature solver */
			strcpy( conv.chSolverEden , "new" );

			/*  constant total pressure, gas+rad+incident continuum
			 *  turn on radiation pressure - off by default in all models - must be
			 * explicitly turned on  */
			pressure.lgRadPresON = TRUE;

			/* increse precision of solution */
			conv.EdenErrorAllowed = 1e-3;
			/* the actual relative error is relative to the total heating and cooling,
			 * which include the dynamics.heat and .cool, which are the advected heating/cooling.
			 * the two terms can be large and nearly cancel, what is written to the .heat and cool files
			 * by punch files has had the smaller of the two subtracted, leaving only the net advected 
			 * heating and cooling */
			conv.HeatCoolRelErrorAllowed = 3e-4f;
			conv.PressureErrorAllowed = 1e-3f;
		}
	}

	/* this is usual hypersonic outflow */
	else
	{
		if( wind.windv0 <= 1.e6 )
		{
			/* speed of sound roughly 10 km/s */
			fprintf( ioQQQ, " >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" );
			wind.lgWindOK = FALSE;
		}

		/* set the central object mass, in solar masses */
		wind.comass = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		if( lgEOL )
		{
			wind.comass = 1.;
		}

	}

	/* this is used in convpres to say wind solution - both cases use this*/
	strcpy( dense.chDenseLaw, "WIND" );

	/*  option to turn off continuum radiative acceleration */
	if( lgMatch("O CO",chCard) )
	{
		pressure.lgContRadPresOn = FALSE;
	}
	else
	{
		pressure.lgContRadPresOn = TRUE;
	}

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

/*DynaPrtZone called to print zone results */
void DynaPrtZone( void )
{

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

	ASSERT( nzone>0 && nzone<struc.nzlim );

	if( nzone > 0 )
	{
		fprintf(ioQQQ," Advection: Uad%.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n",
			timesc.sound_speed_adiabatic/1e5 ,
			wind.windv/1e5 ,
			dynamics.Cool/thermal.ctot,
			dynamics.Heat/thermal.ctot);
	}

	ASSERT( EnthalpyDensity[nzone-1] > 0. );

	fprintf(ioQQQ," Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n",
		phycon.EnergyExcitation,
		phycon.EnergyIonization,
		phycon.EnergyBinding,
  		0.5*POW2(wind.windv)*dense.xMassDensity,
		5./2.*pressure.PresGasCurr ,
		EnthalpyDensity[nzone-1]/dense.gas_phase[ipHYDROGEN] , AdvecSpecificEnthalpy
	);

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

/*DynaPunch punch info related to advection */
void DynaPunch(FILE* ipPnunit , char chJob )
{

	if( chJob=='a' )
	{
		/* this was the original punch advection, now punch dynamics advection */
		fprintf( ipPnunit , "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
			radius.depth,
			thermal.htot , 
			dynamics.Cool , 
			dynamics.Heat , 
			dynamics.dCooldT ,
			dynamics.Source[ipHYDROGEN][ipHYDROGEN],
			dynamics.Rate,
			phycon.EnthalpyDensity/dense.gas_phase[ipHYDROGEN] ,
			AdvecSpecificEnthalpy
			);
	}
	else
		TotalInsanity();
}
