/*GrnVryDpth set grains abundance as a function of depth into cloud*/
STATIC double GrnVryDpth(

/* nd is the number of the grain bin. The values are listed in the Cloudy output,
 * under "Average Grain Properties", and can easily be obtained by doing a trial
 * run without varying the grain abundance and setting stop zone to 1 */

	long int nd)
{
	double GrnVryDpth_v;

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

	/* set grains abundance as a function of depth into cloud
	 * NB most quantities are undefined for first calls to this sub */
	/* nd is the pointer to the grain bin.  This routine must return
	 * a scale factor for the abundance at this position. */

	if( gv.bin[nd]->lgDustVary )
	{
#if 0
		/* sample code - the scale factor will be the hydrogen ionization fraction
		 * following can be used for any other element - first dimension
		 * is atomic number (ipHYDROGEN, ipHELIUM, etc.), and second is stage
		 * of ionization (0=neutral, 1=singly ionized, etc.) - this function is
		 * probably appropriate for PAHs, which do not exist in ionized gas but
		 * are present beyond the ionization front */
		if( dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHYDROGEN][0] > 0. )
			GrnVryDpth_v = dense.xIonDense[ipHYDROGEN][0]/
				(dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHYDROGEN][0]);
		else
			GrnVryDpth_v = 1.;
#else
		if( phycon.te >= 4000. )
			GrnVryDpth_v = 1.e-10;
		else if( phycon.te >= 3000. )
			GrnVryDpth_v = 1.e-10 - 1./PI*atan2((phycon.te-3000.)*(phycon.te-4000.),
							    (phycon.te-3000.)*2000.+(phycon.te-4000.)*2000.);
		else
			GrnVryDpth_v = 1.;

		ASSERT( GrnVryDpth_v > 0. && GrnVryDpth_v <= 1.000001 );
#endif
	}
	else
	{
		GrnVryDpth_v = 1.;
	}


#	ifdef DEBUG_FUN
	fputs( " <->GrnVryDpth()\n", debug_fp );
#	endif
	return GrnVryDpth_v;
}
/*lint +e662 creation of out of bounds pointer */