#!/opt/local/bin/perl
#
#  Name: get_awc.pl -- Determine the AWC for each STATSGO Mapunit
#
#  Purpose:
#
#	This script uses available water capacity (AWC) and supporting
#	data from the STATSGO COMP and LAYER tables to determine
#	weighted full profile AWC values for each mapunit in the STATSGO
#	MAPUNIT table.  The output is in cm (more properly, cubic cm of
#	water per square cm of column cross section).
#
#	Available water capacity is defined in the NRSC Soil Survey
#	Manual as "the volume of water that should be available to
#	plants if the soil, inclusive of rock fragments, were at field
#	capacity".
#
#	Note that the STATSGO AWCL and AWCH values are volumetric
#	fraction, in the range 0.00 to < 1.00.  This script computes
#	total water in columns whose lengths are specified in  file
#	<coldepthtab>.  The values are rounded to the nearest cm to
#	permit packing the results as two-byte integers; for a column
#	length of 100 cm this matches the precision (.01) of the
#	AWCL/AWCH values.
#
#  Usage:
#	perl get_awc.pl coldepthtab workspace outfile rptfile
#
#  Options and Arguments:
#	coldepthtab	File specification for table defining lengths
#			of columns for which AWC values are wanted.
#	workspace	Arc/Info Workspace location of INFO 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. In particular, it was written to access the
#	Layer and Component STATSGO tables to determine the
#	available water capacity on a mapunit basis.  Output from this
#	script is in an ASCII format; script pack_awc.pl may be	used to
#	convert it to a format suitable for ingest into an Arc/Info INFO
#	table defined by script def_awc_items.aml.
#
#  Algorithmic details:
#
#	The STATSGO LAYER table gives a range of available water
#	capacity values (AWCL and AWCH) for each layer of each
#	component.  This script computes the average of these values for
#	each layer, and then sums the values to each of the depths
#	specified by the table pointed to by argument coldepthtab.  If a
#	specified depth is in the middle of a component layer, as
#	specified by LAYER table items LAYDEPH and LAYDEPL, the AWC for
#	the layer is multiplied by the fraction of the layer which is
#	within the depth.  The values for all components of a given
#	mapunit are then averaged, using the COMPPCT values in the COMP
#	table as weights, to determine the weighted average AWC values
#	for the mapunit.
#
#	If a mapunit has one or more components which are all water, and
#	also has other non-water components, the average AWC is computed
#	for the non-water components only.  If the mapunit is entirely
#	water, the AWC is set equal to the column length (i.e., the
#	volumetric fraction of water is set to 100%).  If a component is
#	all rock and there are no layers specified for it, AWC is set to
#	zero.
#
#  External files accessed:
#	<coldepthtab>	READ to get lengths of columns for which AWC
#			values are to be determined.  This file contains
#			one record for each depth; the depths must be in
#			cm and arranged in increasing order.
#	<workspace>	Arc/Info Workspace location of INFO tables being 
#			accessed for processing.
#	<outfile>	WRITE to record AWC for each mapunit.  Each
#			record contains data for all components of one 
#			mapunit.
#	<rptfile>	WRITE to record processing arguments, statistics
#			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_awc.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:
#	help		Output help text.
#	bld_layers	Subroutine to compute AWC for one mapunit
#			component for the layers defined by successive
#			segments of the specified column lengths.
#	output_mu	Generate output values for each component of one
#			mapunit, and reinitialize for next component.
#	set_allrock	Subroutine to build single rock layer when no soil.
#	usage		Output terse usage text.
#
#  Internal variables:
#	COMP		File handle for reading COMP table listing.
#	COLDEP		File handle for reading column depth table.
#	LAYER		File handle for reading LAYER table listing.	 
#	RPT		File handle for output report file.
#	UCONV		Conversion factor from STATSGO depth units to standard
#			depth units.
#	bufc		Array of values from COMP table.
#	bufl		Array of values from LAYER table.
#	cl		Current component layer.
#	clawc		Array of component layer available water capacities.
#	clbot		Array of component layer bottoms, in cm.
#	cldeph		Bottom of component layer, LAYDEPH item.
#	cldepl		Top of component layer, LAYDEPL item.
#	cltop		Array of component layer tops, in cm.
#	cmd		Name of this script; used for help text.
#	coldepthtab	Filename for table defining column depths.
#	compname	Component name from COMP table.
#	comppct		Percent of component in map unit from COMP table.
#	comprec		Number of current COMP record.
#	curmuid		ID of mapunit currently being processed.
#	errcnt		Count of severe errors.
#	lawc		Array of awc accumulators for each output layer.
#	layer		Number of current output layer.
#	laynum		Value of LAYERNUM item.
#	layerrec	Number of current LAYER record
#	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.
#	muid		Mapunit ID.
#	ncl		Number of layers for current component.
#	nextcl		Next component layer.
#	nlayer		Number of output layers (column depths).
#	nonwtotpct	Sum of percents for non water components in map
#			unit.
#	opt		Name of current option.
#	outfile		Name of file to receive relayered output.
#	overlap		Amount of overlap of component and output layers.
#	rdeph		Bedrock depth for component, higher abs. value.
#	rdepl		Bedrock depth for component, lower abs. value.
#	repm		Bedrock depth for component, mean.
#	reclen_comp	Length of COMP table record, in bytes.
#	reclen_layer	Length of LAYER table record, in bytes.
#	savecl		Array to save data for next component layer.
#	seqnum		Component sequence number.
#	stats		Associative array of statistics for anomalies.
#	template_comp	Template for reading COMP table, returned by openinfo
#	template_lay	Template for reading LAYER table, returned by openinfo
#	time		Array of values returned by localtime().
#	totpct		Sum of percents for all components in map unit.
#	value		Value of one array element.
#
#  Script history:
#
#	 4-96	Initial version, D. A. Miller.
#	12-96	Added support for multiple column depths, R. A. White.
#
############   End of Prolog for get_awc.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 coldepthtab workspace outfile rptfile\n"; 
    print "        perl $cmd -u|h\n";

    }		#  End of subroutine usage

sub help
    {
    print "\n";
    print $cmd, " -- Get available water capacity for each STATSGO",
      " mapunit\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 "  coldepthtab	File specification for table defining",
			" lengths\n", 
	  "		of columns for which AWC values are wanted.\n";
    print "   workspace	Arc/Info Workspace location of INFO 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));

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

#  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;

#  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];
	    $clawc[$ncl++] = 0.0;
	    }
	else
	    {
	    $clbot[$ncl - 1] = $rdepm;
	    if ($ltop[$nlayer] > $rdepm)
		{
		$cltop[$ncl] = $clbot[$ncl - 1];
		$clbot[$ncl] = $ltop[$nlayer];
		$clawc[$ncl++] = 0.0;
		}
	    }
	}
    
#  Get available water capacity 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]);
	    $lawc[$layer]
	      += $comppct * $clawc[$cl] * ($overlap / $lthick[$layer]); 
	    $overlap = $ltop[$layer + 1] - $clbot[$cl];
	    }
	$lawc[$layer]
	  += $comppct * $clawc[$cl] * ($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 "");

#  Report problems and anomalies

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

#  Compute total water for each column depth

    for ($layer = 0; $layer < $nlayer; $layer++)
	{
	if ($nonwtotpct > 0)
	    {
	    $lawc[$layer] *= ($lthick[$layer] / $nonwtotpct);
	    $lawc[$layer] += $lawc[$layer - 1] if ($layer > 0);
	    }
	else
	    {
	    print RPT "Non-water total COMPPCT for MUID $curmuid",
	      " = $nonwtotpct -- setting AWC = 0\n" if ($layer == 0);
	    $lawc[$layer] = 0.0;
	    }
	}

#  Write record for this mapunit

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

#  Reinitialize accumlators not done above

    $totpct = 0;
    $nonwtotpct = 100;

    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];
	if (($compname =~ /^\s*WATER\s*$/) && ($comppct > 0) &&
	  ($comppct < 100))
	    {
	    $clawc[0] = 0.0;
	    $nonwtotpct -= $comppct;
	    }
	else
	    {
	    $clawc[0] = 1.0;
	    $comppct = 100;
	    }
	}
    elsif ($rdeph == 0)
	{
	print RPT "No layers for $curmuid $seqnum -- assuming",
	  " solid rock\n";
	$nextcl = $ncl = 1;
	$cltop[0] = 0;
	$clbot[0] = $ltop[$nlayer];
	$clawc[0] = 0.0;
	}
    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;
    }

$coldepthtab = shift;
$workspace = shift;
$outfile = shift;
$rptfile = shift;

# Set program constants

$UCONV = 2.54;	 # Convert inches to centimeters

#  Open column depth file; get layer depths and thicknesses

$nlayer = 0;
$ltop[0] = 0;

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

#  Open	 output files

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_awc.pl processing", $time[5], $time[4]+1,
  $time[3], $time[2], $time[1];

print RPT "Arguments are:\n";
printf RPT "    %-16s%s\n", "coldepthtab", $coldepthtab;
printf RPT "    %-16s%s\n", "workspace", $workspace;
printf RPT "    %-16s%s\n", "outfile", $outfile;
printf RPT "    %-16s%s\n\n", "rptfile", $rptfile;

#  Call openinfo to open and unpack the INFO formatted COMP table

(($reclen_comp,$template_comp) = &openinfo("COMP",$workspace,
     "COMP","MUID","SEQNUM","COMPNAME","COMPPCT","ROCKDEPL","ROCKDEPH"))
     || die "Openinfo returned error";

#  Call openinfo to open and unpack the INFO formatted LAYER table

(($reclen_layer,$template_lay) = &openinfo("LAYER",$workspace,
     "LAYER","MUID","SEQNUM","LAYERNUM","LAYDEPL","LAYDEPH","AWCL","AWCH"))
     || die "Openinfo returned error";

#  Clear all awc entries

for ($layer = 0; $layer < $nlayer; $layer++)
    {
    $lawc[$layer] = 0.0;
    }
   
#  Process components of each mapunit in turn

$curmuid = "";
$totpct = 0;
$nonwtotpct = 100;
$nextcl = 0;
$comprec = 0;
$layerrec = 0;

while (read (COMP, $bufc, $reclen_comp))
    {
    $comprec++;
    (($muid, $seqnum, $compname, $comppct, $rdeph, $rdepl) = 
       unpack ($template_comp, $bufc))
      || 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], $clawc[0]) = @savecl;
	    next;
	    }
	}
    while (read (LAYER, $bufl, $reclen_layer))
	{
	$layerrec++;
	(($lmuid, $lseqnum, $laynum, $cldepl, $cldeph, $awcl, $awch)
	    = unpack ($template_lay, $bufl))
	  || die "Bad LAYER record $layerrec";
	$cltop[$nextcl] = $cldepl * $UCONV;
	$clbot[$nextcl] = $cldeph * $UCONV;

#  Compute available water capacity for this component layer

	if (($awcl eq "") || ($awch eq ""))
	    {
	    print RPT "AWCL or AWCH missing, $lmuid $lseqnum, layer",
	      "$laynum\n";
	    $stats{"noawc"}++;
	    $awcl = 0 if ($awcl eq "");
	    $awch = 0 if ($awch eq "");
	    }
	if ($awcl > $awch)
	    {
	    print RPT "AWCL > AWCH, $lmuid $lseqnum, layer $laynum",
	      " -- layer ignored\n";
	    $stats{"valrev"}++;
	    die "Too many errors -- see report file $rptfile: "
	      if ($errcnt++ > 10);
	    next;
	    }

	$clawc[$nextcl] = (($awcl + $awch)/2.0);

	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], $clawc[$nextcl]);
    $nextcl = &bld_layers;
    ($cltop[0], $clbot[0], $clawc[0]) = @savecl;
    }

&output_mu;		# Output results for last mapunit

close COMP;
close LAYER;
close OUT;

#  Add statistics to processing report.

print RPT "\nProcessing statistics\n\n";
$value = 0 unless ($value = $stats{"totalcomp"});
print RPT "    Total components in input compfile:          $value\n"; 
print RPT "    Number of layers with AWCL > AWCH:           $value\n"
  if (($value = $stats{"valrev"}) > 0);
print RPT "    Number of layers with zero or missing AWC:   $value\n"
  if (($value = $stats{"noawc"}) > 0);

close RPT;

exit 0;		# From Script get_awc.pl
