#!/opt/local/bin/perl
#
#  Name:  get_fragtext.pl -- Get texture and rock fragment classes
#
#  Purpose:
#
#	This script uses soil layer and texture data from the STATSGO
#	COMP and LAYER tables to estimate the dominant soil texture and
#	soil fragment classes for each of a set of standard-thickness
#	soil layers for each mapunit of a state.  This process is
#	necessary because the number and depth of the layers vary widely
#	from one component to another within a given mapunit, and even
#	more widely between mapunits.  The script also determines the
#	mean depth to bedrock for each mapunit.
#
#  Usage:
#	perl get_fragtext.pl layerdeftab fragtextdeftab workspace
#		outfile rptfile
#
#  Options and Arguments:
#	layerdeftab	File specification for table defining
#			thicknesses of standard soil layers.  Each
#			thickness is specified by a separate record;
#			layers are specified from top to bottom, with
#			thicknesses in cm.
#	fragtextdeftab	File specification for table assigning standard
#			coarse rock fragment and soil texture classes to
#			each texture-plus-fragment class used in STATSGO
#			tables.
#	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 an ASCII
#	format; script pack_fragtext.pl may be used to convert it to a
#	format suitable for ingest into Arc/Info INFO tables defined by
#	scripts def_soiltext_items.aml and def_rockfrag_items.aml.
#
#  Algorithmic details:
#
#	The dominant texture and framgent class for each standard layer
#	is derived from the TEXTURE1 values specified in the LAYER table
#	entries for each component.  Although other texture/fragment
#	classes may be present in a given layer, as indicated by the
#	TEXTURE2 and TEXTURE3 variables, no information is provided as
#	to the relative proportions of these other classes.
#	Accordingly, only the TEXTURE1 value is used.
#
#	For each mapunit, the dominant category for each standard layer
#	is determined by weighting the contribution of each component to
#	that layer by the component percent, COMPPCT, value.  When a
#	standard output layer overlaps more than one layer of a given
#	component in the map unit, the dominant texture for each of the
#	overlapped layers is weighted by the amount of the overlap.
#
#	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 this depth, the
#	soil texture and rock fragment modifiers are flagged as BRT and
#	BRF, repectively, to indicate that bedrock may be present, but
#	that the ROCKDEPL/ROCKDEPH values may instead merely indicate
#	the depth below which no soil data are available.
#
#
#  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.
#	fragtextdeftab	READ to get assignment of each texture-plus-
#			fragment-modifier to standard coarse fragment
#			and soil texture classes.  Each record contains
#			one texture-plus-modifier code (ASCII) followed
#			by the standard abbreviations for the standard
#			fragment and texture classes.
#	<workspace>	READ to get values of INFO items MUID, SEQNUM,
#	 /info/COMP	COMPPCT, COMPNAME, SURFTEX, ROCKDEPL, and
#			ROCKDEPH for each component of each mapunit.
#	<workspace>	READ to get values of INFO items MUID, SEQNUM,
#	 /info/LAYER    LAYERNUM, LAYDEPL, LAYDEPH, and TEXTURE1 for
#			each layer of each component.
#	outfile		WRITE to record dominant texture class and
#			rock-fragment class for each standard layer.
#			Each record contains data for all layers of one
#			mapunit, beginning with the mapunit ID and the
#			mean depth to bedrock, followed by the texture
#			classes for all layers and then the coarse
#			fragment classes.
#	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.
#			This variable should be set if openinfo.pl is in
#			neither the directory from which get_fragtext.pl
#			is invoked nor in 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:
#	BRF		Fragment class for bedrock layers.
#	BRT		Texture class for bedrock layers.
#	COMP		File handle for reading COMP table.
#	IN		Filehandle for input file.
#	LAYER		File handle for reading LAYER table.
#	OTHER		Texture and fragment code for type OTHER.
#	OUT		File handle for relayered output.
#	RPT		File handle for output report file.
#	STDFRAG		List of standard rock fragment classes.
#	STDTEXT		List of standard soil texture classes.
#	UCONV		Conversion factor from STATSGO depth units to
#			standard depth units.
#	WATER		Texture and fragment class for surface water.
#	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.
#	clfclass	Array of dom. fragment class for comp. layers. 
#	cltclass	Array of dom. texture class for comp. layers. 
#	cltop		Array of component layer tops, in cm.
#	cltext1		Array of component layer dominant textures.
#	cmd		Name of this script; used for help text.
#	cnmuid		MUID associated with COMPNAME value.
#	cnrec		Current record from compname file.
#	cnseqnum	Component SEQNUM associated with COMPNAME value.
#	compname	Name of current component.
#	comppct		Percent of component in map unit.
#	comprec		Number of current COMP record.
#	comprecl	Record length for COMP table.
#	comptempl	Template for unpacking COMP table.
#	curmuid		ID of mapunit currently being processed.
#	curseq		Sequence number of current component.
#	domfrag		Array of dominant fragment class for each layer.
#	domfragwght	Array of weights for each layer's dominant frag.
#	domtext		Array of dominant texture class for each layer.
#	domtextwght	Array of weights for each layer's dominant text.
#	frag		Current output fragment class.
#	fragmap		Associative array of fragment assignments.
#	fragtextdeftab	Filename for fragment/texture definition table.
#	fragwght	Associative array of weights for each frag type.
#	layerdeftab	Filename for table defining standard layers.
#	layer		Number of current output layer.
#	layerrec	Number of current LAYER record.
#	laynum		Value of LAYERNUM item.
#	layrecl		Record length for LAYER table.
#	laytempl	Template for unpacking LAYER table.
#	lmuid		Mapunit ID of current layer.
#	lseqnum		Component sequence number of current layer.
#	lthick		Array of thicknesses of standard layers.
#	ltop		Array of depths to tops of standard layers.
#	murdep		Mean depth to bedrock for mapunit.
#	muid		Mapunit ID.
#	nextcl		Next component layer.
#	ncl		Number of layers for current component.
#	nlayer		Number of standard layers.
#	opt		Name of current option.
#	outfile		Name of file to receive relayered output.
#	outfrag		Array of output fragment types.
#	outtext		Array of output textures classes.
#	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.
#	surftex		Dominant surface texture of component.
#	text		Current output texture class.
#	text1		Dominant texture for component layer, TEXTURE1.
#	textmap		Associative array of texture assignments.
#	textwght	Associative array of weights for each texture.
#	time		Array of values returned by localtime().
#	totpct		Sum of percents for all components in map unit.
#	value		Value of one array element.
#	workspace	Arc/Info workspace containing tables.
#
#  Script history:
#	09-95	Initial version, R. A. White.
#	01-96	Added check of COMPNAME for "WATER", R. A. White.
#	04-96   Modified for direct read of INFO tables using subroutine
#		openinfo, R. A. White.
#	12-96	Modified to use consolidated fragment/texture definition
#		table, R. A. White.
#
############   End of Prolog for get_fragtext.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 fragtextdeftab",
	" workspace\n",
	"		outfile rptfile\n"; 
    print "       ", "perl ", $cmd, " -u|h\n";

    }		#  End of subroutine usage

sub help
    {
    print "\n";
    print $cmd, " -- Get dominant texture and rock fragment classes\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\n",
	  "		thicknesses of standard soil layers.\n";
    print "fragtextdeftab	File specification for table assigning",
			" standard\n", 
	  "		coarse rock fragment and soil texture classes",
			" to\n",
	  "		each texture-plus-fragment class used in",
			" STATSGO\n",
	  "		tables.\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
    {
    $ncl = $nextcl;
    return 1 if (($ncl == 0) && ( ! &set_allrock));

#  Get texture and fragment class for each component layer

    print RPT "Surface texture mismatch, $curmuid $seqnum --\n",
      "    compfile:  $surftex, layerfile:  $cltext1[0]\n"
      if (($surftex ne $cltext1[0])
        && (($ncl > 1) || ($surftex ne " ") || ($rdeph > 0)));
    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 (! ($cltclass[$cl] = $textmap{$cltext1[$cl]}))
	    {
	    print RPT "No texture class defined for $cltext1[$cl],",
	      " $curmuid $seqnum, layer ";
	    printf RPT "%d\n   -- assigning \"OTHER\"\n", $cl + 1;
	    $cltclass[$cl] = $OTHER;
	    }
	if (! ($clfclass[$cl] = $fragmap{$cltext1[$cl]}))
	    {
	    print RPT "No fragment class defined for $cltext1[$cl],",
	      " $curmuid $seqnum, layer ";
	    printf RPT "%d\n   -- assigning \"OTHER\"\n", $cl + 1;
	    $clfclass[$cl] = $OTHER;
	    }

#  Check for layer overlaps, inversions, etc.

	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"}++;
	    }
	}

#  Convert bedrock depths to cm

    $totpct += $comppct;
    $rdepl = 0 if ($rdepl eq "");
    $rdeph = 0 if ($rdeph eq "");
    $rdepl *= $UCONV;
    $rdeph *= $UCONV;
    $rdepm = ($rdepl + $rdeph) / 2.0;
    $murdep += $comppct * $rdepm;

#  Extend component layer coverage to bottom of output layers

    $stats{"totalcomp"}++;	# Total components
    if ($clbot[$ncl - 1] < $ltop[$nlayer])
	{
	$stats{"lextend"}++;	# Component layer extended downwards
	if ($clbot[$ncl - 1] >= $rdepm)
	    {
	    if (($clbot[$ncl - 1] > $rdepm)
	      && (grep (/^$cltclass[$ncl - 1]$/, @STDTEXT)))
		{
		$stats{"belowbr"}++;	# Bedrock truncates c.l.
		$stats{"belowbrh"}++
		  if ($clbot[$ncl - 1] > $rdeph);
		}
	    $cltop[$ncl] = $clbot[$ncl - 1];
	    $clbot[$ncl] = $ltop[$nlayer];
	    $cltclass[$ncl] = $BRT;
	    $clfclass[$ncl++] = $BRF;
	    }
	else
	    {
	    $stats{"abovebr"}++;
	    $stats{"abovebrl"}++ if ($clbot[$ncl - 1] < $rdepl);
	    $clbot[$ncl - 1] = $rdepm;
	    if ($ltop[$nlayer] > $rdepm)
		{
		$cltop[$ncl] = $clbot[$ncl - 1];
		$clbot[$ncl] = $ltop[$nlayer];
		$cltclass[$ncl] = $BRT;
		$clfclass[$ncl++] = $BRF;
		}
	    }
	}
    

#  Get texture and fragment class 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]);
	    $textwght{$cltclass[$cl], $layer}
	      += ($comppct * $overlap) / $lthick[$layer];
	    $fragwght{$clfclass[$cl], $layer}
	      += ($comppct * $overlap) / $lthick[$layer];
	    $overlap = $ltop[$layer + 1] - $clbot[$cl];
	    }
	$textwght{$cltclass[$cl], $layer}
	  += ($comppct * $overlap) / $lthick[$layer];
	$fragwght{$clfclass[$cl], $layer}
	  += ($comppct * $overlap) / $lthick[$layer];
	}

    return 1;
    }		# End of subroutine bld_layers

#
#  Subroutine to output layer data for one mapunit
#

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

#  Compute mean depth to bedrock

    print RPT "Total COMPPCT for MUID $curmuid = $totpct\n"
      if ((($totpct < 99) || ($totpct > 101))
	&& (++$stats{"badtotpct"}));
    $murdep /= $totpct unless ($totpct == 0);

#  Find dominant texture for each output layer

    for ($layer = 0; $layer < $nlayer; $layer++)
	{
	$domtext[$layer] = "";
	$domtextwght[$layer] = 0.0;
	foreach $text (@outtext)
	    {
	    if ( ($textwght{$text, $layer} > $domtextwght[$layer])
	      && ( ($text ne $WATER) || ($layer == 0)
	      || ($domtext[$layer - 1] eq $WATER) ) )
		{
		$domtextwght[$layer] = $textwght{$text, $layer};
		$domtext[$layer] = $text;
		}
	    $textwght{$text, $layer} = 0.0;
	    }

#  Find dominant rock fragment class for each output layer

	$domfrag[$layer] = "";
	$domfragwght[$layer] = 0.0;
	foreach $frag (@outfrag)
	    {
	    if ( ($fragwght{$frag, $layer} > $domfragwght[$layer])
	      && ( ($frag ne $WATER) || ($layer == 0)
	      || ($domfrag[$layer - 1] eq $WATER) ) )
		{
		$domfragwght[$layer] = $fragwght{$frag, $layer};
		$domfrag[$layer] = $frag;
		}
	    $fragwght{$frag, $layer} = 0.0;
	    }
	}

#  Write record for this mapunit

    printf OUT "%8s%4d", $curmuid, $murdep;
    for ($layer = 0; $layer < $nlayer; $layer++)
	{
	printf OUT "%6s", $domtext[$layer];
	}
    for ($layer = 0; $layer < $nlayer; $layer++)
	{
	printf OUT "%6s", $domfrag[$layer];
	}
    print OUT "\n";

#  Reinitialize accumlators not done above

    $murdep = 0.0;
    $totpct = 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];
	$cltext1[0] = $surftex = $WATER;
	$comppct = 100
	  unless (($compname =~ /^\s*WATER\s*$/) && ($comppct > 0)); 
	}
    elsif (($surftex eq " ") || ($rdeph == 0))
	{
	print RPT "No layers for $curmuid $seqnum -- assuming",
	  " solid rock\n";
	$nextcl = $ncl = 1;
	$cltop[0] = 0;
	$clbot[0] = $ltop[$nlayer];
	$cltext1[0] = ($surftex ne " ") ? $surftex : $BRT;
	}
    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 < 4)
    {
    &usage;
    exit 1;
    }
$layerdeftab = shift;
$fragtextdeftab = shift;
$workspace = shift;
$outfile = shift;
$rptfile = shift;

#  Get layer depths

$UCONV = 2.54;		# Convert inches to centimeters
$BRF = "BR";		# Default classes for bedrock
$BRT = "BR";
$OTHER = "O";
$WATER = "W";
@STDTEXT = ("SL", "C", "CL", "S", "L", "LS", "SC", "SCL", "SI", "SIC",
  "SICL", "SIL");

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

#  Build associative array of fragment and texture definitions

open (IN, $fragtextdeftab)
 || die "Error opening texture definition table, $textdeftab,";
while (<IN>)
    {
    if (/^\s*(\w\S*)\s+(\S+)\s+(\S+)($|\s+)/)
	{
	$textmap{$1} = $2;
	push (@outtext, $2) if ! grep (/^$2$/, @outtext);
	$fragmap{$1} = $3;
	push (@outfrag, $3) if ! grep (/^$2$/, @outfrag);
	}
    }
close IN;
$textmap{$BRT} = $BRT if ( ! defined ($textmap{$BRT}));
push (@outtext, $BRT) if ! grep (/^$BRT$/, @outtext);
$textmap{$WATER} = $WATER if ( ! defined ($textmap{$WATER}));
push (@outtext, $WATER) if ! grep (/^$WATER$/, @outtext);
$fragmap{$BRT} = $BRF if ( ! defined ($fragmap{$BRT}));
push (@outfrag, $BRF) if ! grep (/^$BRF$/, @outfrag);
$fragmap{$WATER} = $WATER if ( ! defined ($fragmap{$WATER}));
push (@outfrag, $WATER) if ! grep (/^$WATER$/, @outfrag);

#  Open component, layer, and output files

(($comprecl, $comptempl) = &openinfo ("COMP", $workspace, "COMP",
  "MUID", "SEQNUM", "COMPNAME", "COMPPCT", "SURFTEX", "ROCKDEPL",
  "ROCKDEPH")) || die "Error opening COMP table"; 
(($layrecl, $laytempl) = &openinfo ("LAYER", $workspace, "LAYER",
  "MUID", "SEQNUM", "LAYERNUM", "LAYDEPL", "LAYDEPH", "TEXTURE1"))
  || die "Error opening LAYER table";
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_fragtext.pl processing", $time[5], $time[4] + 1,
  $time[3], $time[2], $time[1];
print RPT "Arguments are:\n";
printf RPT "    %-16s%s\n", "layerdeftab", $layerdeftab;
printf RPT "    %-16s%s\n", "fragtextdeftab", $fragtextdeftab;
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 weight entries

for ($layer = 0; $layer < $nlayer; $layer++)
    {
    foreach $text (@outtext)
	{
	$textwght{$text,$layer} = 0.0;
	}
    foreach $frag (@outfrag)
	{
	$fragwght{$frag,$layer} = 0.0;
	}
    }

#  Process components of each mapunit in turn

$curmuid = "";
$murdep = 0.0;
$totpct = 0;
$nextcl = 0;

while (read (COMP, $_, $comprecl))
    {
    $comprec++;
    (($muid, $seqnum, $compname, $comppct, $surftex, $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], $cltext1[0]) = @savecl;
	    next;
	    }
	}
    while (read (LAYER, $_, $layrecl))
	{
	$layerrec++;
	(($lmuid, $lseqnum, $laynum, $cldepl, $cldeph, $text1)
	  = unpack ($laytempl, $_)) || die "Bad LAYER record $layerrec";
	$cltop[$nextcl] = $cldepl * $UCONV;
	$clbot[$nextcl] = $cldeph * $UCONV;
	$cltext1[$nextcl] = $text1;
	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], $cltext1[$nextcl]);
    $nextcl = &bld_layers;
    ($cltop[0], $clbot[0], $cltext1[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";
$value = 0 unless ($value = $stats{"belowbr"});
print RPT "    Number of non-rock layers beyond mean bedrock:  $value\n";
$value = 0 unless ($value = $stats{"belowbrh"});
print RPT "    Number of non-rock layers beyond max bedrock:   $value\n";
$value = 0 unless ($value = $stats{"abovebr"});
print RPT "    Number of bottom layers above mean bedrock:     $value\n";
$value = 0 unless ($value = $stats{"abovebrl"});
print RPT "    Number of bottom layers above min bedrock:      $value\n";

close RPT;

exit 0;		# From Script get_fragtext.pl