#!/opt/local/bin/perl
#
#  Name:  get_rockvol.pl -- Get rock volumes for standard layers
#
#  Purpose:
#
#	This script uses soil layer and texture data from the STATSGO
#	COMP and LAYER tables to estimate the percentage by volume of
#	rock fragments for each of a set of standard-thickness soil
#	layers for each mapunit of a state.  See "Algorithmic details",
#	below.
#
#  Usage:
#	perl get_rockvol.pl layerdeftab	workspace outfile rptfile
#
#  Options and Arguments:
#	layerdeftab	File specification for table defining
#			thicknesses of standard soil layers.
#	workspace	Pathname of directory containing info/
#			subdirectory for COMP and LAYER tables.
#	outfile		File to receive output table.
#	rptfile		File to receive report of processing statistics,
#			problems, etc.
#
#  External associations:
#
#	This script is part of a package for processing the STATSGO soil
#	properties database.  Output from this script is in ASCII
#	format.  Script pack_rockvol.pl may be used to convert it to a
#	format suitable for ingest into an Arc/Info INFO table defined
#	by script def_rockvol_items.aml.
#
#  Algorithmic details:
#
#	For each component of each mapunit, the volume percent of rock
#	fragments is estimated for each component layer using the bulk
#	density (mean of BDL and BDH), percent rock fragments larger
#	than 10 inches (INCH10L/H), percent rock fragments between 3 and
#	10 inches (INCH3L/H), and percent of material less than 3 inches
#	which passes a No. 10 sieve (NO10L/H).  Each standard output
#	layer is assigned the volume percent of rock fragments value for
#	the component layer within which it is contained; if the
#	standard output layer overlaps more than one component layer,
#	it's value is computed by a weighted average of the values for
#	the overlapped component layers, with the weights proportional
#	to the amount of overlap.  The average volume percent of rock
#	fragments over the entire mapunit for a given standard output
#	layer is then computed by summing the products of the value for
#	each component times the component percent, COMPPCT.
#
#	The region from the bottom of the last component layer to the
#	bottom of the last standard layer, if any, is assumed to be the
#	same as the lowest component layer down to the mean bedrock
#	depth (mean of ROCKDEPL and ROCKDEPH), below which the rock
#	volume percent is set to 100.  Because the STATSGO data
#	frequently use a rock depth of 60 inches to indicate that no
#	soil information was available below this depth, rock volume
#	values for layer extending below this depth will often be
#	misleading.
#
#	When the rock volume for some, but not all, components of a map
#	unit cannot be calculated because of missing or contradictory
#	data values or because the component is specified to be water,
#	these components are ignored when computing the average percent.
#	If all non-water components of the map unit contain invalid or
#	contradictory information, the rock volume is assigned a flag
#	value of 101%.  If the map unit is entirely water, the rock
#	volume is set to zero.
#
#  External files accessed:
#	<layerdeftab>	READ to get thickness of each standard layer.
#			This file contains one record for each layer,
#			starting at the top, giving its thickness in cm.
#	<workspace>	READ to get values of INFO items MUID, SEQNUM,
#	 /info/COMP	COMPPCT, ROCKDEPL, and ROCKDEPH for each
#			component of each mapunit.
#	<workspace>	READ to get values of INFO items MUID, SEQNUM,
#	 /info/LAYER    LAYERNUM, LAYDEPL, LAYDEPH, INCH10L, INCH10H,
#			INCH3L, INCH3H, NO10L, NO10H, BDL, and BDH for
#			each layer of each component.
#	<outfile>	WRITE to record volumetric percent of rock for
#			each standard layer.  Each record contains data
#			for all layers of one mapunit.
#	<rptfile>	WRITE to record processing arguments, statistics
#			of bedrock/bottom-layer mismatches, and other
#			problems, etc., encountered during processing.
#
#  Environment variables referenced:
#	PERLINC		Libraries to be searched by require() ahead of
#			those specified by -I and the default libraries.
#                                  
#  External subroutines included:
#	openinfo.pl	Open an INFO table, returning unpack() template
#			for reading a specified list of items.
#
#  Internal subroutines defined:
#	bld_layers	Build layers for one mapunit component.
#	help		Output help text.
#	output_mu	Generate output values for each layer of one
#			mapunit, and reinitialize for next layer.
#	set_allrock	Build single rock layer when no soil.
#	usage		Output terse usage text.
#
#  Internal variables:
#	COMP		File handle for reading COMP table.
#	FLAGACTION	Associative array specifying actions for flags.
#	LAYDEF		Filehandle for reading layer definition table.
#	LAYER		Filehandle for readeing LAYER table.
#	OUT		File handle for writing computed rock fraction.
#	PD		Particle density of rock.
#	RPT		File handle for output report file.
#	UCONV		Conversion factor from STATSGO depth units to
#			standard depth units.
#	allnull		Number of layers with all null fragment values.
#	bdl		Lower and upper ends of range of values for
#	bdh		 moist bulk density of soil.
#	bdm		Mean moist bulk density of soil.
#	cl		Current component layer.
#	clbot		Array of component layer bottoms, in cm.
#	cldeph		Bottom of component layer, LAYDEPH item.
#	cldepl		Top of component layer, LAYDEPL item.
#	clflag		Array of component layer problem flags.
#	cltop		Array of component layer tops, in cm.
#	clvrf		Array of component layer volume fraction of rock
#			fragments.
#	cmd		Name of this script; used for help text.
#	comprec		Number of current COMP record.
#	comprecl	Record length for COMP table.
#	comptempl	Template for unpacking COMP table.
#	comppct		Percent of component in map unit.
#	curmuid		ID of mapunit currently being processed.
#	denom		Denominator for computing vol. fract. of fines.
#	errcnt		Count of severe errors.
#	inch3l		Lower and upper ends of range of percent by
#	inch3h		 weight of 3 to 10 inch rock fragments.
#	inch3m		Mean percent of 3 to 10 inch rock fragments.
#	inch10l		Lower and upper ends of range of percent by
#	inch10h		 weight of rock fragments larger than 10 inches.
#	inch10m		Mean percent of rock fragments > than 10 inches.
#	l2laynum	LAYNUM value for current l2rec.
#	l2muid		MUID value for current l2rec.
#	l2rec		Current record of second part of LAYER listing.
#	l2seqnum	SEQNUM value for current l2rec.
#	layerdeftab	Filename for table defining standard layers.
#	layer		Number of current output layer.
#	laynum		Value of LAYERNUM item.
#	layerrec	Number of current LAYER record.
#	layrecl		Record length for LAYER table.
#	laytempl	Template for unpacking LAYER table.
#	lmuid		Mapunit ID of current layer.
#	lowbd		Number of layers with 0.1 <= BDL/H <= .90.
#	lseqnum		Component sequence number of current layer.
#	lthick		Array of thicknesses of standard layers.
#	ltop		Array of depths to tops of standard layers.
#	lvrf		Array of rock fraction accumulators for each
#			standard output layer.
#	muid		Mapunit ID.
#	no10l		Lower and upper ends of range of percent of soil
#	no10h		 less than 3 inches which passes No. 10 filter.
#	no10m		Mean percent of soil < 3" passing No. 10 filter.
#	ncl		Number of layers for current component.
#	nextcl		Next component layer.
#	nlayer		Number of standard layers.
#	opt		Name of current option.
#	outfile		Name of file to receive relayered output.
#	overlap		Amount of overlap of component & output layers. 
#	rdeph		Bedrock depth for component, higher abs value.
#	rdepl		Bedrock depth for component, lower abs value.
#	rdepm		Bedrock depth for component, mean.
#	rptfile		Name of file to receive report.
#	savecl		Array to save data for next component layer.
#	seqnum		Component sequence number.
#	stats		Associative array of statistics for anomalies.
#	time		Array of values returned by localtime().
#	totpct		Sum of percents for all components in map unit.
#	vals		Array of soil data values for one layer.
#	valtotpct	Sum of percents of components permitting valid
#			(usable) rock volume computation.
#	value		Value of one array element.
#	varnames	Array of soil variable names.
#	vindex		Loop counter/index for soil variables.
#	watermu		Flag set if mapunit is all water.
#	workspace	Arc/Info workspace containing tables.
#
#  Script history:
#	10-95	Initial version, R. A. White
#	03-96	Added support for all-water components, used openinfo.pl
#		to access INFO tables.
#
############   End of Prolog for get_rockvol.pl   ############

#  Load external library modules, adding $PERLINC to search path

unshift (@INC, split (/ /,$ENV{"PERLINC"}));
require "openinfo.pl";

#  Define usage and help text

@cmd =  (split ('/', $0));
$cmd = pop (@cmd);

sub usage
    {
    print "\n";
    print "Usage: ", "perl ", $cmd, " layerdeftab workspace",
	" outfile rptfile\n";
    print "       ", "perl ", $cmd, " -u|h\n";

    }		#  End of subroutine usage

sub help
    {
    print "\n";
    print $cmd, " -- Get rock volume percents for standard layers\n";
    &usage;
    print "\n";
    print "Options and arguments are\n";
    print "     -h   	Write this help text to stdout and exit.\n";
    print "     -u   	Display terse usage message and exit.\n";
    print " layerdeftab	File specification for table defining",
	  " thicknesses\n",
	  "		of standard soil layers.\n";
    print " workspace	Pathname of directory containing info/\n",
	  "		subdirectory for COMP and LAYER tables.\n";
    print " outfile	File to receive output table.\n",
	  " rptfile	File to receive report of processing",
			" statistics,\n",
	  "		problems, etc.\n";

    }		#  End of subroutine help

#
#  Subroutine to build layers for one mapunit component
#

sub bld_layers
    {
    $stats{"totalcomp"}++;	# Total components
    $ncl = $nextcl;
    return 1 if (($ncl == 0) && ( ! &set_allrock));
    $totpct += $comppct;

#  Convert bedrock depths to cm

    if (($rdepl eq "") && ($rdeph eq ""))
	{
	if (($clflag[0] ne "water") && ($clflag[0] ne "allrock"))
	    {
	    print RPT "Mapunit $curmuid, component $seqnum:	no",
	      " ROCKDEP specified -- component ignored\n";
	    $stats{"norock"}++;
	    }
	$rdepl = $rdeph  = 0;
	}
    $rdepl *= $UCONV;
    $rdeph *= $UCONV;
    $rdepm = ($rdepl + $rdeph) / 2.0;

#  Determine whether component contains valid data

    if ($clflag[0] eq "water")
	{
	$watermu = 1 if ($comppct == 100);
	return 1;
	}
    for ($cl = 0; $cl < $ncl; $cl++)
	{
	return 1 if ($FLAGACTION{$clflag[$cl]} eq "skipcomp");
	if ($cltop[$cl] >= $rdepm)
	    {
	    if ($FLAGACTION{$clflag[$cl]} eq "checkbr")
		{
		$value = $cl + 1;
		if ($clflag[$cl] eq "rocksnofines")
		    {
#		    print RPT "Mapunit $curmuid, component $seqnum:",
#		      "  layer $value has INCH3/10 and\n  no NO10L/H",
#		      " data but below mean bedrock -- assuming all",
#		      " rock\n";
		    $stats{"rnfbr"}++;
		    }
		else
		    {
#		    print RPT "Mapunit $curmuid, component $seqnum:",
#		      "  layer $value has no NO10L/H data\n   but below",
#		      " mean bedrock -- assuming all rock\n";
		    $stats{"no10br"}++;
		    }
		$clvrf[$cl] = 1.0;
		}
	    else
		{
		$stats{"belowbr"}++ if ($cltop[$cl] < 1.0);
		}
	    }
	elsif ($clflag[$cl] eq "rocksnofines")
	    {
	    $stats{rnfcomp}++;
	    return 1;
	    }
	}

#  Check for layer overlaps, inversion, etc.

    if ($cltop[0] != 0.0)
	{
	print RPT "Layer 1 of $curmuid $seqnum starts at $cltop[0]",
	  " -- changed to 0.0\n";
	$cltop[0] = 0.0;
	$stats{"nonzero"}++;
	}
    for ($cl = 0; $cl < $ncl; $cl++)
	{
	if ($cltop[$cl] > $clbot[$cl])
	    {
	    print RPT "LAYDEPH = $clbot[$cl] less than LAYDEPL =",
	      " $cltop[$cl], $curmuid $seqnum, layer $cl\n";
	    $clbot[$cl] = $cltop[$cl];
	    die "Too many errors -- see report file $rptfile: "
	      if ($errcnt++ > 10);
	    }	    
	if (($cl > 0) && ($cltop[$cl] < $clbot[$cl - 1]))
	    {
	    if (($cltop[$cl] == $cltop[$cl -1])
	      && ($clbot[$cl] == $clbot[$cl - 1]))
		{
		$value = $cl + 1;
		print RPT "Layer $value duplicates $cl, $curmuid",
		  " $seqnum -- using first only\n";
		$cltop[$cl] = $clbot[$cl];
		$stats{"duplayer"}++;
		}
	    elsif ($cltop[$cl] > $cltop[$cl - 1])
		{
		$value = $cl + 1;
		print RPT "Layer $value overlaps layer $cl, $curmuid",
		  " $seqnum -- moving top of layer $value downwards\n";
		$cltop[$cl] = $clbot[$cl - 1];
		$clbot[$cl] = $cltop[$cl]
		  if ($clbot[$cl] < $cltop[$cl]); 
		$stats{"overlap"}++;
		}
	    else
		{
		$value = $cl + 1;
		print RPT "Layer $value starts above $cl, $curmuid",
		  " $seqnum -- moving layer $value downwards\n";
		$cltop[$cl] = $clbot[$cl - 1];
		$clbot[$cl] = $cltop[$cl]
		  if ($clbot[$cl] < $cltop[$cl]); 
		die "Too many errors -- see report file $rptfile: "
		  if ($errcnt++ > 10);
		}
	    }
	elsif (($cl > 0) && ($cltop[$cl] > $clbot[$cl - 1]))
	    {
	    $value = $cl + 1;
	    print RPT "Gap between layer $cl and layer $value,",
	      " $curmuid $seqnum -- extending layer $cl downwards\n";
	    $clbot[$cl - 1] = $cltop[$cl];
	    $stats{"layergap"}++;
	    }
	}

#  Extend component layer coverage to bottom of output layers

    if ($clbot[$ncl - 1] < $ltop[$nlayer])
	{
	$stats{"lextend"}++;	# Component layer extended downwards
	if ($clbot[$ncl - 1] >= $rdepm)
	    {
	    $cltop[$ncl] = $clbot[$ncl - 1];
	    $clbot[$ncl] = $ltop[$nlayer];
	    $clvrf[$ncl++] = 1.0;
	    }
	else
	    {
	    $clbot[$ncl - 1] = $rdepm;
	    if ($ltop[$nlayer] > $rdepm)
		{
		$cltop[$ncl] = $clbot[$ncl - 1];
		$clbot[$ncl] = $ltop[$nlayer];
		$clvrf[$ncl++] = 1.0;
		}
	    }
	}
    

#  Get rock fragment volume contributions to output layers

    $cl = 0;
    for ($layer = 0; $layer < $nlayer; $layer++)
	{
	$overlap = $lthick[$layer];
	for ( ; $clbot[$cl] < $ltop[$layer + 1]; $cl++)
	    {
	    $overlap -= ($ltop[$layer + 1] - $clbot[$cl]);
	    $lvrf[$layer]
	      += $comppct * $clvrf[$cl] * ($overlap / $lthick[$layer]); 
	    $overlap = $ltop[$layer + 1] - $clbot[$cl];
	    }
	$lvrf[$layer]
	  += $comppct * $clvrf[$cl] * ($overlap / $lthick[$layer]); 
	}
    $valtotpct += $comppct;

    return 1;
    }		# End of subroutine bld_layers

#
#  Subroutine to output layer data for one mapunit
#

sub output_mu
    {
    return $muid if ($curmuid eq "");

#  Report problems and anomalies

    print RPT "Total COMPPCT for MUID $curmuid = $totpct\n"
      if ((($totpct < 99) || ($totpct > 101))
	&& (++$stats{"badtotpct"}));
    if ($allnull)
	{
#	print RPT "$curmuid has $allnull layers with no size data -- ",
#	  "assumed rock fraction = 0\n";
	$stats{"allnull"} += $allnull;
	$allnull = 0;
	}
    if ($allzero)
	{
#	print RPT "$curmuid has $allzero layers with BD = NO10",
#	  " = INCH3/10 = 0 -- assumed rock fraction = 0\n";
	$stats{"allzero"} += $allzero;
	$allzero = 0;
	}
    if ($lowbd)
	{
#	print RPT "$curmuid has $lowbd layers with BDL >= .01 and",
#	  " BDH <= .90\n";
	$stats{"lowbd"} += $lowbd;
	$lowbd = 0;
	}
    if ($valtotpct <= 0)
	{
	if ($watermu)
	    {
	    $watermu = 0;
	    $value = 0.0;
	    }
	else
	    {
	    $value = 101.0;
	    print RPT "Mapunit $curmuid:  No components with usable",
	      " data -- setting rockvol = 101%\n";
	    $stats{"nousable"}++;
	    }
	$valtotpct = 100;
	for ($layer =1; $layer < $nlayer; $layer++)
	    {
	    $lvrf[$layer] = $value;
	    }
	}

#  Write record for this mapunit

    printf OUT "%-8s", $curmuid;
    for ($layer = 0; $layer < $nlayer; $layer++)
	{
	printf OUT "%4d", int (100.0 * ($lvrf[$layer] / $valtotpct) + 0.5);
	$lvrf[$layer] = 0.0;
	}
    print OUT "\n";

#  Reinitialize accumlators not done above

    $totpct = 0;
    $valtotpct = 0;

    return $muid;
    }			#  End of subroutine output_mu

#
#  Subroutine to build single rock or water layer when no soil
#

sub set_allrock
    {
    if (($compname =~ /^\s*WATER\s*$/) || ($curmuid =~ /^[A-Z]{2}W/))
	{
	$nextcl = $ncl = 1;
	$cltop[0] = 0;
	$clbot[0] = $ltop[$nlayer];
	$clvrf[0] = 0.0;
	$clflag[0] = "water";
	$comppct = 100
	  unless (($compname =~ /^\s*WATER\s*$/) && ($comppct > 0)); 
	}
    elsif ($rdeph == 0)
	{
	print RPT "No layers for $curmuid $seqnum -- assuming",
	  " solid rock\n";
	$nextcl = $ncl = 1;
	$cltop[0] = 0;
	$clbot[0] = $ltop[$nlayer];
	$clvrf[0] = 1.0;
	$clflag[0] = "allrock";
	}
    else
	{
	print RPT "MUID/SEQNUM mismatch -- ",
	  "compfile:  $curmuid $seqnum,  ",
	  "layerfile:  $lmuid $lseqnum\n";
	die "Too many errors -- see report file $rptfile: "
	  if ($errcnt++ > 10);
	return 0;
	}
    return 1;
    }			# End of subroutine set_allrock

#
#  Start of main routine
#

#  Get options

while ($ARGV[0] =~ /-(\w+)(\s+|$)/)
    {
    push (@opt, $1);
    shift;
    }	
if ($opt[0] eq 'u')
    {
    &usage;
    exit;
    }
if ($opt[0] eq 'h')
    {
    &help;
    exit;
    }

#  Preprocess arguments

if  ($#ARGV < 3)
    {
    &usage;
    exit 1;
    }
$layerdeftab = shift;
$workspace = shift;
$outfile = shift;
$rptfile = shift;

#  Set program constants

$UCONV = 2.54;		# Convert inches to centimeters
$PD = 2.65;		# Use typical density of silicate rocks
%FLAGACTION = ("allnull", "checkbr",
	       "allzero", "checkbr",
	       "bdnofines", "checkbr",
	       "toomuchcoarse", "skipcomp",
	       "valrev", "skipcomp",
#	       "rocksnofines", "skipcomp",
	       "rocksnofines", "checkbr",
	       "finesnobd", "skipcomp",
	       "bdnonzero", "skipcomp",
	       "baddenom", "skipcomp",
	       "bigrocks", "",
	       "water", "",
	       "allrock", "",
	       "", "");

#  Get layer depths


$nlayer = 0;
$ltop[0] = 0;
open (LAYDEF, $layerdeftab)
  || die "Error opening layer table, $layerdeftab,";
while (<LAYDEF>)
    {
    if (/^\s*(\d+)($|\s+)/)
	{
	$lthick[$nlayer++] = $1;
	$ltop[$nlayer] = $ltop[$nlayer - 1] + $1;
	}
    }
close LAYDEF;

#  Open component, layer, and output files

(($comprecl, $comptempl) = &openinfo ("COMP", $workspace, "COMP",
  "MUID", "SEQNUM", "COMPNAME", "COMPPCT", "ROCKDEPL", "ROCKDEPH"))
  || die "Error opening COMP table"; 
(($layrecl, $laytempl) = &openinfo ("LAYER", $workspace, "LAYER",
  "MUID", "SEQNUM", "LAYERNUM", "LAYDEPL", "LAYDEPH", "INCH10L",
  "INCH10H", "INCH3L", "INCH3H", "NO10L", "NO10H", "BDL", "BDH"))
  || die "Error opening LAYER table";
@varnames = ("INCH10", "INCH3", "NO10", "BD");
open (OUT, "> $outfile") || die "Error opening output file, $outfile,";
open (RPT, "> $rptfile") || die "Error opening report file, $rptfile,";

#  Create report header

@time = localtime (time);
$time[5] += 1900;
printf RPT "	%s, %02d/%02d/%02d %02d:%02d:%02d\n\n",
  "Summary of get_rockvol.pl processing", $time[5], $time[4],
  $time[3], $time[2], $time[1];
print RPT "Arguments are:\n";
printf RPT "    %-16s%s\n", "layerdeftab", $layerdeftab;
printf RPT "    %-16s%s\n", "workspace", $workspace;
printf RPT "    %-16s%s\n", "outfile", $outfile;
printf RPT "    %-16s%s\n\n", "rptfile", $rptfile;

#  Clear all layer rock fraction entries

for ($layer = 0; $layer < $nlayer; $layer++)
    {
    $lvrf[$layer] = 0.0;
    }

#  Process components of each mapunit in turn

$allnull = 0;
$allzero = 0;
$lowbd = 0;
$curmuid = "";
$totpct = 0;
$valtotpct = 0;
$watermu = 0;
$nextcl = 0;
$comprec = 0;
$layerrec = 0;

while (read (COMP, $_, $comprecl))
    {
    $comprec++;
    (($muid, $seqnum, $compname, $comppct, $rdepl, $rdeph)
      = unpack ($comptempl, $_)) || die "Bad COMP record $comprec";
    $curmuid = &output_mu if ($muid ne $curmuid);

#  Collect layer data for current component

    if (defined $laynum)
	{
	print RPT "Bad LAYERNUM \"$laynum\" for layer 1,",
	  " $lmuid $lseqnum\n"
	  if (($laynum != 1) && ( ! eof (LAYER)));
	if (($lmuid ne $curmuid) || ($lseqnum != $seqnum))
	    {
	    & bld_layers if (&set_allrock);
	    ($cltop[0], $clbot[0], $clvrf[0]) = @savecl;
	    next;
	    }
	}
    while (read (LAYER, $_, $layrecl))
	{
	$layerrec++;
	$clflag[$nextcl] = "";
	(($lmuid, $lseqnum, $laynum, $cldepl, $cldeph, @vals)
	  = unpack ($laytempl, $_)) || die "Bad LAYER record $layerrec";
	$cltop[$nextcl] = $cldepl * $UCONV;
	$clbot[$nextcl] = $cldeph * $UCONV;

#  Compute rock fragment volume for this component layer

#  First check for problems and compute means

	@missing = (0, 0, 0, 0);
	for ($vindex = 0; $vindex < 8; $vindex++)
	    {
	    if ($vals[$vindex] eq "")
		{
		$vals[$vindex] = 0;
		$missing[$vindex/2]++;
		}
	    }
	$stats{"nobd"}++ if ($missing[3]);
	$lowbd++ if (($vals[6] >= .01) && ($vals[7] <= .90));
	if ($vals[0] + $vals[2] > 100)
	    {
	    print RPT "Mapunit $lmuid, component $lseqnum:  layer",
	      " $laynum has INCH3L + INCH10L > 100%\n  -- ignoring",
	      " component\n";
	    $clflag[$nextcl] = "toomuchcoarse";
	    $stats{"toomuchcoarse"}++;
	    }
	for ($vindex = 1; $vindex < 8; $vindex += 2)
	    {
	    if ($vals[$vindex -1] > $vals[$vindex])
		{
		print RPT "Mapunit $lmuid, component $lseqnum:  ",
		  "layer $laynum has $varnames[$vindex/2]L",
		  " > $varnames[$vindex/2]H -- ignoring component\n";
		$clflag[$nextcl] = "valrev";
		$stats{"valrev"}++;
		die "Too many errors -- see report file $rptfile: "
		  if ($errcnt++ > 10);
		}
	    $vals[$vindex/2] = ($vals[$vindex-1] + $vals[$vindex])
	      / (($vindex < 7) ? 200.0 : 2.0);
	    }
	$#vals = 3;

#  Determine rock volume for layer

	($inch10m, $inch3m, $no10m, $bdm) = @vals;
	$rft = $no10m * (1.0 - $inch3m - $inch10m);
	$denom = $rft + (1.0 - $rft) * ($bdm / $PD);
	if ($FLAGACTION{$clflag[$nextcl]} eq "skipcomp")
	    {
			#  skip computation, since component ignored
	    }
	elsif ($missing[2])	# No entry for NO10L/H
	    {
	    if ($missing[0] && $missing[1] && $missing[3])
		{
		$clvrf[$nextcl] = 0.0;
		$clflag[$nextcl] = "allnull";
		$allnull++;
		}
	    elsif (( ! $missing[0]) || ( ! $missing[1]))
		{
#		print RPT "Mapunit $lmuid, component $lseqnum:  ",
#		  "layer $laynum has INCH3 or INCH10 specified,",
#		  " but not NO10 -- ignoring component\n";
		$clvrf[$nextcl] = 1.0;
		$clflag[$nextcl] = "rocksnofines";
		$stats{"rocksnofines"}++;
		}
	    elsif ( ! $missing[3])	# BDL/H specified
		{
		$clvrf[$nextcl] = 0.0;
		$clflag[$nextcl] = "bdnofines";
		$stats{"bdnofines"}++;
		}
	    }
	elsif ($no10m == 0.0)
	    {
	    if (($inch10m > 0.0) || ($inch3m > 0.0))
		{
		$clvrf[$nextcl] = 1.0;
		$clflag[$nextcl] = "bigrocks";
		$stats{"bigrocks"}++;
		}
	    elsif ($missing[3])	 # BDL/H not specified
		{
		print RPT "Mapunit $lmuid, component $lseqnum:  ",
		  "layer $laynum specifies NO10L/H = 0 but not BDL/H",
		  " -- ignoring component\n";
		$clvrf[$nextcl] = 1.0;
		$clflag[$nextcl] = "finesnobd";
		$stats{"finesnobd"}++;
		}
	    elsif ($bdm > 0)
		{
		print RPT "Mapunit $lmuid, component $lseqnum:  ",
		  "layer $laynum specifies NO10L/H = 0, BDL/H > 0",
		  " -- ignoring component\n";
		$clvrf[$nextcl] = 1.0;
		$clflag[$nextcl] = "bdnonzero";
		$stats{"bdnonzero"}++;
		}
	    else  # ($bdm == 0.0)
		{
		$clvrf[$nextcl] = 1.0;
		$clflag[$nextcl] = "allzero";
		$allzero++;
		}
	    }
	elsif ($inch10m + $inch3m >= 1.0)
	    {
	    if ($inch10m + $inch3m > 1.0)
		{
		print RPT "Mapunit $lmuid, component $lseqnum:  layer",
		  " $laynum has INCH3 + INCH10 > 100%\n  -- reducing",
		  " to 100%\n";
		}
	    else
		{
		print RPT "Mapunit $lmuid, component $lseqnum:  layer",
		  " $laynum all coarse rock > 3 inches\n";
		}
	    $clvrf[$nextcl] = 1.0;
	    $clvrf[$nextcl] = "allcoarse";
	    $stats{"allcoarse"}++;
	    }
	elsif ($denom <= 0.0)
	    {
	    print RPT "Mapunit $lmuid, component $lseqnum:  Bad denom",
	      " = $denom, layer $laynum\n  -- ignoring component\n";
	    $clvrf[$nextcl] = 0.0;
	    $clflag[$nextcl] = "baddenom";
	    die "Too many errors -- see report file $rptfile: "
	      if ($errcnt++ > 10);
	    }
	else
	    {
	    $clvrf[$nextcl] = 1.0 - ($rft / $denom);
	    }
	last if (($lseqnum != $seqnum) || ($lmuid ne $curmuid)); 
	print RPT "Bad LAYERNUM \"$laynum\" for layer",
	  " $nextcl, $lmuid $lseqnum\n"
	  if (++$nextcl != $laynum);
	}
    @savecl = ($cltop[$nextcl], $clbot[$nextcl], $clvrf[$nextcl],
      $clflag[$nextcl]);
    $nextcl = &bld_layers;
    ($cltop[0], $clbot[0], $clvrf[0], $clflag[0]) = @savecl;
    }

&output_mu;		# Output results for last mapunit

close COMP;
close LAYER;
close OUT;

#  Add statisitics to processing report.

print RPT "\nProcessing statistics\n\n";
$value = 0 unless ($value = $stats{"totalcomp"});
print RPT "    Total components in input compfile:            $value\n";
$value = 0 unless ($value = $stats{"lextend"});
print RPT "    Number of bottom layers extended downwards:    $value\n"; 
print RPT "    Number of layers with no size data:            $value\n",
	  "      (assumed rock fraction = 0)\n",
  if ($value = $stats{"allnull"});
print RPT "    Number of layers with BD and size data = 0:    $value\n",
	  "      (assumed rock fraction = 0 unless below bedrock)\n",
  if ($value = $stats{"allzero"});
print RPT "    Number of layers with 0.01 <= BDL/H <= .90:    $value\n"
  if ($value = $stats{"lowbd"});
print RPT "    Number of components w/o rock depth:           $value\n",
	  "      (component not used)\n"
  if ($value = $stats{"norock"});
print RPT "    Number of layers with high value < low value:  $value\n",
	  "      (layer ignored)\n"
  if ($value = $stats{"valrev"});
print RPT "    Number of layers with INCH3/10, no NO10:       $value\n",
	  "      (entire component not used unless below bedrock)\n"
  if ($value = $stats{"rocksnofines"});
print RPT "    Number of above layers below bedrock:          $value\n",
	  "      (assumed all rock)\n"
  if ($value = $stats{"rnfbr"});
print RPT "    Number of components w. such layers above BR:  $value\n"
  if ($value = $stats{"rnfcomp"});
print RPT "    Number of layers specifying BD, but no NO10:   $value\n",
	  "      (assumed all fines unless below bedrock)\n"
  if ($value = $stats{"bdnofines"});
print RPT "    Number of layers with INCH3/10 > 0, NO10 = 0:  $value\n"
  if ($value = $stats{"bigrocks"});
print RPT "    Number of layers with INCH3 + INCH10 > 100:    $value\n",
	  "      (entire component not used)\n"
  if ($value = $stats{"toomuchcoarse"});
print RPT "    Number of layers with INCH3 + INCH10 = 100:    $value\n"
  if ($value = $stats{"allcoarse"});
print RPT "    Number of layers with BD > 0, NO10 = 0:        $value\n",
	  "     (entire component not used)\n"
  if ($value = $stats{"bdnonzero"});
print RPT "    Number of layers below bedrock, NO10 = 0/null: $value\n",
	  "     (assumed all rock)\n"
  if ($value = $stats{"no10br"});
print RPT "    Number of layers below bedrock, rockvol < 1.0: $value\n",
	  "     (used computed value)\n"
  if ($value = $stats{"belowbr"});
print RPT "    Number of layers with BD not specified:        $value\n"
  if ($value = $stats{"nobd"});
print RPT "    Number of mapunits with no usable components:  $value\n"
  if ($value = $stats{"nousable"});

close RPT;

exit 0;		# From Script get_rockvol.pl
