/*HydroLevel solve for ionization balance level populations of model hydrogen atom */
#include <string.h>
#include <float.h>
#include "cddefines.h"
#include "taulines.h"
#include "showme.h"
#include "zonecnt.h"
#include "hcolrt.h"
#include "secondaries.h"
#include "hphoto.h"
#include "charexc.h"
#include "trace.h"
#include "elementnames.h"
#include "phycon.h"
#include "hcbr.h"
#include "hstat.h"
#include "thlo.h"
#include "hcaseb.h"
#include "hlteok.h"
#include "destcrt.h"
#include "rcota.h"
#include "hionfraccom.h"
#include "elecden.h"
#include "ionrange.h"
#include "printe82.h"
#include "negcon.h"
#include "hydrogenic.h"

void HydroLevel(long int ipZ)
{
	char chType[5];
	long int i, 
	  ipHi, 
	  ipLo, 
	  level;
	double BigError ,
	  bottom,
	  colfrc, 
	  phtfrc, 
	  rfac, 
	  secfrc;

	/* this block of variables will be obtained and freed here */
	long int *ipiv ; /* malloc out to [nhlevel+1] */
	double *totcap;/* malloc out to [nhlevel+1]*/
	double **SaveZ/*[nhlevel+2][nhlevel+2]*/, 
	  *bvec/*[nhlevel+2]*/, 
	  *error/*[nhlevel+2]*/, 
	  *work/*[nhlevel+2]*/, 
	  **z/*[nhlevel+2][nhlevel+2]*/;

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

	/* check that we were called with valid charge */
	assert( ipZ >= 0);
	assert( ipZ < LIMELM );

	/* 
	 * following electron density has approximate correction for neutrals
	 * corr of hi*1.7e-4 accounts for col ion by HI; Drawin Zs Phys 225, 483.
	 * used EdenHCorr instead
	 * edhi = eden + hi * 1.7e-4
	 */

	/* now malloc some scratch space */
	ipiv = (long int *)malloc(sizeof(long int)*(unsigned)(nhlevel+1) ) ;
	totcap = (double *)malloc(sizeof(double)*(unsigned)(nhlevel+1) ) ;
	bvec = (double *)malloc(sizeof(double)*(unsigned)(nhlevel+2) ) ;
	error = (double *)malloc(sizeof(double)*(unsigned)(nhlevel+2) ) ;
	work = (double *)malloc(sizeof(double)*(unsigned)(nhlevel+2) ) ;
	if( ipiv==NULL || totcap==NULL || bvec==NULL || error==NULL || work==NULL )
	{
		fprintf( ioQQQ, " could not malloc 1D arrays\n" );
		puts( "[Stop in hydrolevel]" );
		exit(1);
	}

	/* now do the 2D arrays */
	if( (SaveZ = (double **)malloc(sizeof(double *)*(unsigned)(nhlevel+2) ) ) ==NULL )
	{
		fprintf( ioQQQ, " could not malloc SaveZ array in 1D\n" );
		puts( "[Stop in hydrolevel]" );
		exit(1);
	}

	if(  (z = (double **)malloc(sizeof(double *)*(unsigned)(nhlevel+2) ))==NULL  )
	{
		fprintf( ioQQQ, " could not malloc z arrays in 1D\n" );
		puts( "[Stop in hydrolevel]" );
		exit(1);
	}

	/* now do the second dimension */
	for( i=0; i<(nhlevel+2); ++i )
	{
		SaveZ[i] = (double *)malloc(sizeof(double)*(unsigned)(nhlevel+2) ) ;
		if( SaveZ[i]==NULL )
		{
			fprintf( ioQQQ, " could not malloc 2D SaveZ\n" );
			puts( "[Stop in hydrolevel]" );
			exit(1);
		}

		z[i] = (double *)malloc(sizeof(double)*(unsigned)(nhlevel+2) ) ;
		if( z[i]==NULL )
		{
			fprintf( ioQQQ, " could not malloc 2D z\n" );
			puts( "[Stop in hydrolevel]" );
			exit(1);
		}
	}

	/* remember creation and destruction rates */
	DestCrt.destroy[ipZ][ipZ] = HPhoto.hgamnc[ipZ][IP1S] + 
	  Secondaries.csupra + hydro.hcolnc[ipZ][IP1S]*phycon.eden;

	DestCrt.create[ipZ][ipZ] = phycon.eden*hcaseb.htotrc[ipZ];

	/* check how low temperature is; departure coef tend to infinity
	 * as TE goes to zero; matrix formalism won't work on a 32-bit machine
	 * if EXP(H NU/KT) > 1E38 for 13.6EV.
	 * the recombination coefficient is effect radiative (only) rec
	 * coef (optical depths in) */

	/* get simple estimate of ionization balance, first get denominator,
	 * which can be zero during search phase */
	bottom = phycon.eden*(hcaseb.htotrc[ipZ] + 
	  rcota.CotaRate[ipZ]) + CharExc.CTHrec[MIN2(ipZ,1)][1];

	if( bottom > 0. )
	{
		HIonFraccom.HIonSimple[ipZ] = (HPhoto.hgamnc[ipZ][IP1S] + 
		  Secondaries.csupra + hydro.hcolnc[ipZ][IP1S]*ElecDen.EdenHCorr + 
		  CharExc.CTHion[MIN2(ipZ+1,2)-1][1])/bottom;
	}
	else
	{
		HIonFraccom.HIonSimple[ipZ] = 0.;
	}

	if( trace.lgHBug && trace.lgTrace )
	{
		/* identify how atom is ionized for full trace */
		if( HIonFraccom.HIonSimple[ipZ] > 0. )
		{
			/* fraction of ionization due to photoionization */
			phtfrc = HPhoto.hgamnc[ipZ][IP1S]/((phycon.eden*(hcaseb.htotrc[ipZ] + 
			  rcota.CotaRate[ipZ]) + CharExc.CTHrec[MIN2(ipZ+1,2)-1][1])*
			  HIonFraccom.HIonSimple[ipZ]);

			/* fraction of ionization due to collisional ionization */
			colfrc = (hydro.hcolnc[ipZ][IP1S]*ElecDen.EdenHCorr + 
			  CharExc.CTHion[MIN2(ipZ+1,2)-1][1])/((phycon.eden*(hcaseb.htotrc[ipZ] + 
			  rcota.CotaRate[0]) + CharExc.CTHrec[MIN2(ipZ+1,2)-1][1])*
			  HIonFraccom.HIonSimple[ipZ]);

			/* fraction of ionization due to secondary ionization */
			secfrc = Secondaries.csupra/((phycon.eden*(hcaseb.htotrc[ipZ] + 
			  rcota.CotaRate[0]) + CharExc.CTHrec[MIN2(ipZ+1,2)-1][1])*
			  HIonFraccom.HIonSimple[ipZ]);
		}
		else
		{
			phtfrc = 0.;
			colfrc = 0.;
			secfrc = 0.;
		}

		fprintf( ioQQQ, "     HydroLevel Z=%2ld called, simple II/I=",ipZ);
		PrintE93( ioQQQ, HIonFraccom.HIonSimple[ipZ]);
		fprintf( ioQQQ," PhotFrc:");
		PrintE93( ioQQQ,phtfrc);
		fprintf(ioQQQ," ColFrc:");
		PrintE93( ioQQQ,colfrc);
		fprintf( ioQQQ," SecFrc");
		PrintE93( ioQQQ, secfrc);
		fprintf( ioQQQ,"  Te:");
		PrintE93( ioQQQ,phycon.te);
		fprintf( ioQQQ," eden:");
		PrintE93( ioQQQ,phycon.eden);
		fprintf( ioQQQ,"\n"); 
	}

	/* for dense models close to lte all ionizations can be from
	 * excited states, and simple and actual pops are very different.
	 * code used simple all the time, caused catastrophe for Kingdon nova model.
	 * now use actual pops if we are into model */
	if( ZoneCnt.nzone > 2 )
	{
		/* HIonFrac saves ion to neutral, used in bidiag to set hydrogenic ratio */
		rfac = HIonFraccom.HIonFrac[ipZ];
	}
	else
	{
		rfac = HIonFraccom.HIonSimple[ipZ];
	}

	/* end of prelimaries - start model atom */

	/* which case atom to solve??? */
	if( HIonFraccom.HIonSimple[ipZ] < 1e-30 )
	{
		/* don't bother if no ionizing radiation */
		strcpy( chType, "zero" );
		if( trace.lgHBug && trace.lgTrace )
		{
			fprintf( ioQQQ, "     HydroLevel Z=%2ld simple II/I=%10.2e so not doing equilibrium.\n", 
			  ipZ, HIonFraccom.HIonSimple[ipZ] );
		}

		for( i=IP1S; i <= nhlevel; i++ )
		{
			hydro.hbn[ipZ][i] = 1.;
			hydro.hn[ipZ][i] = 0.;
			error[i] = 0.;
		}

		hcaseb.HRecEffec[ipZ] = 0.;
		HIonFraccom.HIonFrac[ipZ] = 0.;
		HIonFraccom.HIonSimple[ipZ] = 0.;

	}
	/* chng 99 nov 23, very dense model close to lte had very low simple
	 * ionization fraction but moderate ionization since all pops in
	 * excited states - added test for density 
	 * second change - for all cases use this routine if ionization is
	 * very low indeed */
	else if( (HIonFraccom.HIonSimple[ipZ] < 1e-15) ||
		(HIonFraccom.HIonSimple[ipZ] < 1e-7 && phycon.hden < 1e6 ) )
	{
		strcpy( chType, "TLow" );
		/* this avoids departure coefficients, 
		 * by doing simple cascade, should always work, 
		 * but not very well at high densities */
		HydroT2Low(ipZ,rfac);

		/* error not set in HydroT2Low so set to zero here */
		for( i=IP1S; i <= nhlevel; i++ )
		{
			error[i] = 0.;
		}
	}
	/* >>chng 99 Jun 08 rfac from 1e-7 to 1e-5 due to orion nebula crash using
	 * costar stellar atmosphere, in He+ which was barely present */
	else if(
		(phycon.te <= thlo.HydTempLimit* (ipZ+1)) || 
		(rfac < 1e-5) || 
	   !HlteOk.lgHlteOk[ipZ] )
	{
		/* this is the low temperature branch, does not use dep coef */
		strcpy( chType, "Pops" );
		HydroLevelPop(ipZ ,
		  SaveZ/*[nhlevel+2][nhlevel+2]*/, 
		  bvec/*[nhlevel+2]*/, 
		  error/*[nhlevel+2]*/, 
		  work/*[nhlevel+2]*/, 
		  z/*[nhlevel+2][nhlevel+2]*/ ,
		  ipiv , /* malloc out to [nhlevel+1] */
		  totcap/* malloc out to [nhlevel+1]*/);
	}
	else
	{
		/* the main solver for most conditions */
		strcpy( chType, "Depr" );
		/* this is the high temperature branch, with dep coef */
		HydroLevelDep(ipZ ,
		  SaveZ/*[nhlevel+2][nhlevel+2]*/, 
		  bvec/*[nhlevel+2]*/, 
		  error/*[nhlevel+2]*/, 
		  work/*[nhlevel+2]*/, 
		  z/*[nhlevel+2][nhlevel+2]*/ ,
		  ipiv , /* malloc out to [nhlevel+1] */
		  totcap/* malloc out to [nhlevel+1]*/);
	}
	/* end branch for main solve */

	/* remember largest residual in matrix inversion */
	BigError = 0.;
	/* all three cases end up down here */
	for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ )
	{
		double abserror;
		abserror = fabs( error[ipLo]) ;
		/* this will be the largest residual in the matrix inversion */
		BigError = MAX2( abserror , BigError );

		for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
		{
			/* population of lower level rel to ion */
			HydroLines[ipZ][ipHi][ipLo].PopLo = 
				hydro.hn[ipZ][ipLo];

			/* population of upper level rel to ion */
			HydroLines[ipZ][ipHi][ipLo].PopHi = 
				hydro.hn[ipZ][ipHi];

			/* population of lower level rel to ion, corrected for stim em */
			HydroLines[ipZ][ipHi][ipLo].PopOpc = 
				(hydro.hn[ipZ][ipLo] - 
				HydroLines[ipZ][ipHi][ipLo].PopHi*
				hstat.HStatWght[ipLo]/hstat.HStatWght[ipHi]);
		}
	}

	/* matrix inversion should be nearly as good as the accuracy of a double,
	 * but demand that it is better than epsilon for a float */
	if( BigError > FLT_EPSILON ) 
	{
		fprintf(ioQQQ,
			" warning - largest residual in hydrogenic %s ipZ=%li matrix inversion is %g "
			"matrix type was %s \n",
			elementnames.chElementName[ipZ],ipZ , BigError , chType );
		ShowMe();
		puts( "[Stop in HydroLevel]" );
		exit(1);
	}

	/* now get real effective rec coef */
	hcaseb.HRecNet[ipZ] = 0.;
	hcaseb.HRecCol[ipZ] = 0.;
	hcaseb.HRecInd[ipZ] = 0.;

	for( i=IP2S; i <= nhlevel; i++ )
	{
		hcaseb.HRecInd[ipZ] += (float)HPhoto.rinduc[ipZ][i];

		hcaseb.HRecCol[ipZ] += (float)(hydro.hcolnc[ipZ][i]*phycon.eden*
		  hcolrt.hlte[ipZ][i]*(1. - hydro.hbn[ipZ][i]));

		hcaseb.HRecNet[ipZ] += (float)(hydro.hrec[ipZ][i][ipRecRad]*
		  hydro.hrec[ipZ][i][ipRecNetEsc] + HPhoto.rinduc[ipZ][i] + 
		  hydro.hcolnc[ipZ][i]*phycon.eden*hcolrt.hlte[ipZ][i]*
		  (1. - hydro.hbn[ipZ][i]));
	}

	/* this will eventually become the ratio of ionized to neutral hydrogenic 
	 * species, create sum of level pops per ion first */
	HIonFraccom.HIonFrac[ipZ] = 0.;
	for( level=IP1S; level <= nhlevel; level++ )
	{
		HIonFraccom.HIonFrac[ipZ] += hydro.hn[ipZ][level];
	}

	if( HIonFraccom.HIonFrac[ipZ] < 0. )
	{
		fprintf( ioQQQ, 
			" HydroLevel %s finds negative hydrogen ion fraction for ipZ=%2ld val= %10.3e simple=%10.3e TE=%10.3e ZONE=%4ld\n", 
		  chType , ipZ, HIonFraccom.HIonFrac[ipZ], 
		  HIonFraccom.HIonSimple[ipZ], 
		  phycon.te, ZoneCnt.nzone );
		fprintf( ioQQQ, " level pop are:" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%8.1e", hydro.hn[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );
		negcon();
		ShowMe();
		puts( "[Stop in HydroLevel]" );
		exit(1);
	}

	/* get level populations, two cases, 
	 * first, may be zero since all cases drop down to here, 
	 * this branch of if is for trivial abundance, so zero out species*/
	else if( HIonFraccom.HIonFrac[ipZ] >= 0. && HIonFraccom.HIonFrac[ipZ] < 1e-30 )
	{
		HIonFraccom.HIonFrac[ipZ] = 0.;

		/* reset pointer to one lower stage of ionization so this not
		 * considered again, hydrogenic considered if IonHigh is ipZ+2 */
		IonRange.IonHigh[ipZ] = ipZ+1;

		/* now zero this out */
		for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
			{
				/* population of lower level rel to ion */
				HydroLines[ipZ][ipHi][ipLo].PopLo = 0.;

				/* population of upper level rel to ion */
				HydroLines[ipZ][ipHi][ipLo].PopHi =  0.;

				/* population of lower level rel to ion, corrected for stim em */
				HydroLines[ipZ][ipHi][ipLo].PopOpc =  0.;

				/* local ots destruction, cm^-3 s^-1 */
				HydroLines[ipZ][ipHi][ipLo].ots =  0.;
			}
		}
	}
	/* this is the main branch that applies when we have non-trivial abundance */
	else
	{
		/* level populations 
		 * this brance we have significant population as sum of levels per ion,
		 * store inverse, ion per atom, here.  For non-H species this is ratio
		 * of highest to next highest stage of ionization */
		HIonFraccom.HIonFrac[ipZ] = 1./HIonFraccom.HIonFrac[ipZ];
	}

	if( (trace.lgHBugFull && trace.lgTrace) && (ipZ == trace.ipZTrace) )
	{
		fprintf( ioQQQ, "       HLEV HGAMNC" );
		PrintE93( ioQQQ, HPhoto.hgamnc[ipZ][IP1S] );
		for( i=IP2S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", HPhoto.hgamnc[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV TOTCAP" );
		for( i=1; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", totcap[i] ));
		}
		fprintf( ioQQQ,PrintEfmt("%10.2e", hcaseb.HRecNet[ipZ] ) );
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV IND Rc" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", HPhoto.rinduc[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		/* incuded recombination rate coefficients */
		fprintf( ioQQQ, "       IND Rc LTE " );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e",
				HPhoto.hgamnc[ipZ][i]*hcolrt.hlte[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		/* induced rates down */
		fprintf( ioQQQ, "       HLEV IND As" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", hcbr.hbul[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		/* LTE level populations */
		fprintf( ioQQQ, "       HLEV   HLTE" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", hcolrt.hlte[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		/* fraction of total ionization due to photo from given level */
		fprintf( ioQQQ, "       HLEVfr cion" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", 
				hydro.hcolnc[ipZ][i]*
			  hydro.hn[ipZ][i]*phycon.eden/(hcaseb.HRecEffec[ipZ]*phycon.eden) ) );
		}
		fprintf( ioQQQ, "\n" );

		/* fraction of total ionization due to photoionization from given level */
		if( hcaseb.HRecEffec[ipZ]> 0. )
		{
			fprintf( ioQQQ, "       HLEVfrPhIon" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", 
					HPhoto.hgamnc[ipZ][i]*hydro.hn[0][i]/(hcaseb.HRecEffec[ipZ]*phycon.eden)) );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "       HLEV     HN" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.hn[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV   b(n)" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.hbn[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV   X12tot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", Secondaries.x12tot ));
		fprintf( ioQQQ," Grn dest:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",
		  DestCrt.destroy[ipZ][ipZ] ));
		fprintf(ioQQQ, "\n"); 
	}

	/* check for non-positive level populations */
	if( HIonFraccom.HIonFrac[ipZ] > 0. )
	{
		for( level=IP1S; level <= nhlevel; level++ )
		{
			if( hydro.hn[ipZ][level] <= 0. )
			{
				fprintf( ioQQQ, 
					" HydroLevel finds negative hydrogen level population for %s ipZ=%ld level %ld value=%10.3e simple=%10.3e TE=%10.3e ZONE=%4ld\n", 
				  elementnames.chElementName[ipZ], ipZ, level, hydro.hn[ipZ][level], 
				  HIonFraccom.HIonSimple[ipZ], 
				  phycon.te, ZoneCnt.nzone );
				fprintf( ioQQQ, " level pop are:" );
				for( i=IP1S; i <= nhlevel; i++ )
				{
					fprintf( ioQQQ,PrintEfmt("%8.1e", hydro.hn[ipZ][i] ));
				}
				fprintf( ioQQQ, "\n" );
				negcon();
				ShowMe();
				puts( "[Stop in HydroLevel]" );
				exit(1);
			}
		}
	}
	else if( HIonFraccom.HIonSimple[ipZ] > 0 && 
		HIonFraccom.HIonFrac[ipZ] <= 0.)
	{
		/* this is case where we expected ionization, but found none, or negative */
		fprintf( ioQQQ, 
			" HydroLevel finds non-positive hydrogen ion frac for ipZ=%2ld value=%10.3e simple=%10.3e TE=%10.3e ZONE=%4ld\n", 
		  ipZ, hydro.hn[ipZ][level], 
		  HIonFraccom.HIonSimple[ipZ], 
		  phycon.te, ZoneCnt.nzone );
		fprintf( ioQQQ, " level pop are:" );
		for( i=IP1S; i <= nhlevel; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%8.1e", hydro.hn[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );
		negcon();
		ShowMe();
		puts( "[Stop in HydroLevel]" );
		exit(1);
	}

	if( trace.lgTrace )
	{
		/* HRecEffec(ipZ) is gross rec coef, computed here while filling in matrix
		 * elements, all physical processes included. 
		 * htotrc(ipZ) is radiavive only */
		fprintf( ioQQQ, "       HydroLevel%3ld retrn %4.4s te=",
			ipZ, 
			chType );
		PrintE93( ioQQQ,phycon.te);
		fprintf( ioQQQ," HII/HI=");
		PrintE93( ioQQQ,HIonFraccom.HIonFrac[ipZ]);
		fprintf( ioQQQ," simple=");
		PrintE93( ioQQQ,HIonFraccom.HIonSimple[ipZ]);
		fprintf( ioQQQ," b1=");
		PrintE82( ioQQQ,hydro.hbn[ipZ][IP1S]);
		fprintf( ioQQQ," ion rate=");
		PrintE82( ioQQQ,HPhoto.hgamnc[ipZ][IP1S] + Secondaries.csupra);
		fprintf( ioQQQ," TotRec");
		PrintE82( ioQQQ,hcaseb.HRecEffec[ipZ]);
		fprintf( ioQQQ," RadRec");
		PrintE82( ioQQQ,hcaseb.htotrc[ipZ]);
		fprintf( ioQQQ, "\n");
	}

	/* update destruction rate */
	DestCrt.destroy[ipZ][ipZ] = HPhoto.hgamnc[ipZ][IP1S] + 
	  Secondaries.csupra + hydro.hcolnc[ipZ][IP1S]*phycon.eden;
	DestCrt.create[ipZ][ipZ] = phycon.eden*hcaseb.htotrc[ipZ];

	/* now free up the arrays */
	free( ipiv );
	free( totcap);
	free( bvec );
	free( error );
	free( work );

	/* now free up the 2D arrays */
	for( i=0; i<nhlevel+2; ++i )
	{
		free( SaveZ[i] ) ;
		free( z[i] ) ;
	}
	free( SaveZ ) ;
	free( z ) ;

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