#!/usr/bin/perl
#======================================================================
#                    L A D C P _ W _ C T D 
#                    doc: Mon Nov  3 17:34:19 2014
#                    dlm: Mon Jan 25 09:32:32 2016
#                    (c) 2014 A.M. Thurnherr
#                    uE-Info: 51 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';

# HISTORY:
#	Nov  3, 2014: - created
#	Nov  4, 2014: - improved
#	Nov  6, 2014: - BUG: sound speed was not calculated correctly
#				  - added -a
#				  - added conductivity & temperature editing
#	Nov  7, 2014: - loosened outlier editing
#				  - added no-valid-data error message
#				  - modified binning criterion to allow any sampling
#					frequency (not just divisors of 24)
#	Apr 17, 2015: - added in-situ temperature to output
#	Apr 23, 2015: - had to losen unrealistic-soundspeed criterion
#					because of DUK2#45
#	May 13, 2015: - BUG: did not work with casts beginning at depth
#				  - added errors when editing edits too much
#	May 14, 2015: - BUG: depth was wrong for partial-depth casts
#				  - BUG: $badval had not been handled correctly
#				  - BUG: editing errors were too tight
#	Jun 16, 2015: - added STDOUT redirection (default output file) on tty stdout
#				  - added -p)lot option
#	Jun 17, 2015: - added cond to std output
#				  - added conductivity spike filter (different from pressure)
#	Jun 18, 2015: - renamed from SBE_w
#	Jun 19, 2015: - added press & salin to make output suitable for u/v processing
#				  - removed a couple of data assertions
#	Jul 20, 2015: - allow for non-numeric -i
#	Aug 17, 2015: - modified plots
#	Sep 21, 2015: - added lat/lon to output if available in input
#	Sep 22, 2015: - added error on non-numeric station information
#				  - changed plot output name
#	Sep 23, 2015: - added -o
#				  - improved pressure editing
#				  - modified editing algorithm, primarily reducing temperature resolution
#					from 0.1deg to 1deg and finding contiguous range up to single gap in histogram
#					(requring halving of resolution)
#	Oct 12, 2015: - require ANTSlibs V6.2 for release
#	Oct 13, 2015: - added code to deal with CNV files with wrong number of scans in header
#   Oct 13, 2015: - adapted to [version.pl]
#   Jan 25, 2016: - added software version %PARAM

($ANTS)  = (`which ANTSlib`   =~ m{^(.*)/[^/]*$});
($WCALC) = ($0                =~ m{^(.*)/[^/]*$});
$WCALC   = '.' if ($WCALC eq '');

require "$WCALC/version.pl";
require "$ANTS/ants.pl";
require "$ANTS/fft.pl";
require "$ANTS/libGMT.pl";
require "$ANTS/libSBE.pl";
require "$ANTS/libEOS83.pl";
&antsAddParams('LADCP_w_CTD',"Version $VERSION");

$antsParseHeader = 0;											# usage
&antsUsage('ai:l:orp:s:v:',1,
	'[-v)erbosity <level[0]>]',
	'[use -a)lternate sensor pair]',
	'[-r)etain all data (no editing)] [allow infinite -o)utliers]',
	'[-s)ampling <rate[6Hz]>]',
	'[-l)owpass w <cutoff[2s]>]',
	'[profile -i)d <id>]',
	'[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]',
	'<SBE CNV file>');

&antsFloatOpt(\$opt_l,2);										# defaults
&antsCardOpt(\$opt_s,6);
&antsCardOpt(\$opt_v,0);

$CNVfile = $ARGV[0];											# open CNV file
open(F,&antsFileArg());
&antsActivateOut();												# activate ANTS file

#----------------------------------------------------------------------
# Read Data
#----------------------------------------------------------------------

print(STDERR "Reading $CNVfile...") if ($opt_v);

($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) =			# decode SBE header	
    SBE_parseHeader(F,0,0);										# SBE field names, no time check

croak("$CNVfile: unexpected sampling interval $sampint\n")
	unless (abs($sampint-1/24) < 1e-5);
croak("$CNVfile: no latitude in header\n")
	unless numberp($lat);

$pressF = fnr('prDM');

if ($opt_a) {													# temp/cond alternate sensor pair
	$tempF	= fnr('t190C');
	$condF	= fnrNoErr('c1S/m');
	if (defined($condF)) {
		$condHistRes = 20;										# 0.2 S/m bins
	} else {
		$condF  = fnr('c1mS/cm');
		$condHistRes = 2;										# 2.0 mS/cm bins
	}
} else {														# primary sensor pair
	$tempF	= fnr('t090C');
	$condF	= fnrNoErr('c0S/m');
	if (defined($condF)) {
		$condHistRes = 20;
	} else {
		$condF  = fnr('c0mS/cm');
		$condHistRes = 2;
	}
}

$latF = &fnrNoErr('lat');										# GPS data if available (to make files useful for u/v processing)
$lonF = &fnrNoErr('lon');

&antsInstallBufFull(0);											# read entire CNV file
&SBEin(F,$ftype,$nfields,$nscans,$badval);

if (@ants_ != $nscans) {
	if ($opt_v > 1) {
		printf(STDERR "\n\nWARNING: $CNVfile has wrong number of scans in header\n");
	} else {
		printf(STDERR "WARNING: $CNVfile has wrong number of scans in header\n");
	}
	$nscans = @ants_;
}

printf(STDERR "\n\t%d scans",$nscans) if ($opt_v > 1);
printf(STDERR "\n") if ($opt_v);

#----------------------------------------------------------------------
# Redirect STDOUT to %.6Hz & create %_w_CTD.ps,%_sspd.ps if STDOUT is a tty
#----------------------------------------------------------------------

$id = defined($opt_i) ? $opt_i : &antsParam('station');
croak("$CNVfile: no station information in header => -i required\n")
	unless defined($id);
croak("$CNVfile: non-numeric station information <$id> in header => -i required\n")
	unless numberp($id);
	
if (-t STDOUT) {
	if (numberp($id)) {
		$opt_p = '%03d_w_CTD.ps,%03d_sspd.ps'
			unless defined($opt_p);
		$outfile = sprintf('%03d.%dHz',$id,$opt_s);
	} else {
		$opt_p = '%s_w_CTD.ps,%s_sspd.ps'
			unless defined($opt_p);
		$outfile = sprintf('%s.%dHz',$id,$opt_s);
	}
    open(STDOUT,">$outfile") || die("$outfile: $!\n");
}

#----------------------------------------------------------------------
# Edit Data
#	- pressure outliers & spikes
#	- conductivity outliers & spikes
#	- max_bin+1 is returned to ensure full range
#----------------------------------------------------------------------

sub validRange($)
{
	my($guess_bin) = @_;
	my($min_bin,$max_bin);

	croak("assertion failed") unless ($hist[$guess_bin]);
	for ($max_bin=$guess_bin; $hist[$max_bin]; $max_bin++) { }
	for ($min_bin=$guess_bin; $hist[$min_bin]; $min_bin--) { }
	return ($min_bin,$max_bin+1);
}

unless ($opt_r) {
	print(STDERR "Editing Data...") if ($opt_v);

	#----------------------------------------
	# trim initial records with nan pressures
	#----------------------------------------
	my($trimmed) = 0;												# trim leading nan pressures
	shift(@ants_),$trimmed++
		until numberp($ants_[0][$pressF]) && numberp($ants_[0][$condF]);
	printf(STDERR "\n\t%d initial records with nan pressure and/or conductivity trimmed",$trimmed) if ($opt_v > 1);
	my($lvp) = $ants_[0][$pressF];
	my($lvc) = $ants_[0][$condF];
	
	#------------------------------------------------
	# edit pressure outliers outside contiguous range
	#	- 2dbar resolution
	#	- histogram shifted by 100dbar to allow for negative values
	#------------------------------------------------
	my($outliers) = my($modeSamp) = 0; my($modeBin,$min,$max); local(@hist);
	for (my($s)=0; $s<$nscans; $s++) {
		next unless ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);
		my($b) = ($ants_[$s][$pressF]+100) / 2;
		$hist[$b]++;
		next unless ($hist[$b] > $modeSamp);
		$modeSamp = $hist[$b]; $modeBin = $b;
	}
	printf(STDERR "\n\tvalid pressure guess: %d dbar (%d samples)",2*$modeBin-100,$modeSamp)
		if ($opt_v > 1);
	($min,$max) = validRange($modeBin);
	$min = 2*$min-100; $max = 2*$max-100;
	for (my($s)=0; $s<$nscans; $s++) {
		if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }
		if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }
	}
	&antsAddParams("pressure_outliers",sprintf("%d",$outliers));
	printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
	croak("$CNVfile: pressure editing removed too many 'outliers'\n")
		unless ($opt_o || $outliers < 100);
	
	#----------------------------------------------------
	# edit conductivity outliers outside contiguous range
	#----------------------------------------------------
	$outliers = $modeSamp = 0;
	undef(@hist);
	for (my($s)=0; $s<$nscans; $s++) {
		next unless ($ants_[$s][$condF] > 0);
		my($b) = $ants_[$s][$condF]*$condHistRes;					# 1/10 S/m histogram resolution (1 mS/cm)
		$hist[$b]++;
		next unless ($hist[$b] > $modeSamp);
		$modeSamp = $hist[$b]; $modeBin = $b;
	}
	($min,$max) = validRange($modeBin);
	$min /= $condHistRes; $max /= $condHistRes; 
	for (my($s)=0; $s<$nscans; $s++) {
		if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; }
		if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; }
	}
	&antsAddParams("conductivity_outliers",sprintf("%d",$outliers));
	printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
	croak("$CNVfile: conductivity editing removed too many 'outliers'\n")
		unless ($opt_o || $outliers/$nscans < 0.4);

	#----------------------------------------------------
	# edit temperature outliers outside contiguous range
	#	- Stan's NBP0901 profiles require resolution of 1deg
	#	- otherwise 0.2deg seems to be fine
	#----------------------------------------------------
	$outliers = $modeSamp = 0;
	undef(@hist);
	for (my($s)=0; $s<$nscans; $s++) {
		next unless ($ants_[$s][$tempF] > 0);
		my($b) = $ants_[$s][$tempF]*1;								# 10th of a degree histogram resolution
		$hist[$b]++;
		next unless ($hist[$b] > $modeSamp);
		$modeSamp = $hist[$b]; $modeBin = $b;
	}
	printf(STDERR "\n\ttemperature mode: %.1f degC (%d samples)",$modeBin/10,$modeSamp)
		if ($opt_v > 1);
	($min,$max) = validRange($modeBin);
	$min /= 1; $max /= 1; 
	for (my($s)=0; $s<$nscans; $s++) {
		if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; }
		if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; }
	}
	&antsAddParams("temperature_outliers",sprintf("%d",$outliers));
	printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers)
		if ($opt_v > 1);
	croak("$CNVfile: temperature editing removed too many 'outliers'\n")
		unless ($opt_o || $outliers/$nscans < 0.4);

	#----------------------------------------
	# edit pressure spikes based on gradients
	#----------------------------------------
	
	for (my($s)=1; $s<$nscans; $s++) {								# calculate pressure gradients (across gaps)
		if (numberp($ants_[$s][$pressF])) {
			$dp[$s-1] = $ants_[$s][$pressF] - $lvp;
			$lvp = $ants_[$s][$pressF];
		} else {
			$dp[$s-1] = nan;
		}
	}

	my($ns1,$ns2) = (0,0);
	for (my($s)=0; $s<$nscans-2; $s++) {							# consecutive large pressure gradients of opposite sign
		if (($dp[$s]*$dp[$s+1] < 0) &&								# tests return false if either of the dps is not defined
			(abs($dp[$s]) > 10) &&
			(abs($dp[$s+1]) > 10)) {
				$ants_[$s+1][$pressF] = nan;
				$dp[$s] = $dp[$s+1] = undef;
				$ns1++;
		}
	}
	for (my($s)=0; $s<$nscans-3; $s++) {							# 3 consecutive large pressure gradients of opposite sign
		if (($dp[$s]>2	&& $dp[$s+1]<-4 && $dp[$s+2]>2) ||
			($dp[$s]<-2 && $dp[$s+1]>4	&& $dp[$s+2]<-2)) {
				$ants_[$s+1][$pressF] = $ants_[$s+2][$pressF] = nan;
				$dp[$s] = $dp[$s+1] = $dp[$s+2] = undef;
				$ns2+=2;
		}
	}
	&antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2));
	printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v > 1);

	#--------------------------------------------------
	# edit conductivity spikes based on large gradients
	#	- $lvc = $ants_[0] is guaranteed numeric here
	#--------------------------------------------------
	
	my($nv) = my($ns) = my($last_dcond) = my($lvcs) = 0;
	my($run_start) = my($run_dcond) = 0;
	for (my($s)=1; $s<$nscans; $s++) {								# calculate conductivity gradients (across gaps)
		next unless numberp($ants_[$s][$condF]);

		my($dcond) = $ants_[$s][$condF] - $lvc;						# integrate gradient across runs
#		print(STDERR "ldc/dc/rdc: $last_dcond/$dcond/$run_dcond\n");
		if ($last_dcond*$dcond >= 0) {								# run is continuing (gradient does not change sign)
			$run_dcond += $dcond;
		} else {													# run has ended
#			print(STDERR "new run at $lvcs (run_dcond = $last_dcond)\n");
			$run_start = $lvcs;
			$run_dcond = $dcond;
		}

		if (abs($run_dcond) <= 0.02) {								# small integrated gradient => okay
			$lvc = $ants_[$s][$condF];								# update stored previous values
			$lvcs = $s;
			$last_dcond = $dcond;
			next;													# process next scan
		}

#		print(STDERR "large gradient ($ants_[$run_start][$condF]-$ants_[$s][$condF], $run_dcond) run $run_start-$s\n");
		my($i);
		my($max_spike_length) = 24;
		for ($i=$s; $i<=$run_start+$max_spike_length && $i<$nscans-1; $i++) {	# large gradient => check whether values return within 10s
			next unless (numberp($ants_[$i][$condF]) &&
						 numberp($ants_[$i+1][$condF]));
			last if ((abs($ants_[$i][$condF]-$ants_[$run_start][$condF]) <= 0.005) &&	# 2 vals to avoid large gradients
					 (abs($ants_[$i+1][$condF]-$ants_[$run_start][$condF]) <= 0.005));
		}

		if ($i>$run_start+$max_spike_length || $i==$nscans-1) {		# values don't return => leave data alone
#			print(STDERR "values don't return => new run at $s\n");
			$run_start = $s; $run_dcond = 0;
			$lvc = $ants_[$run_start][$condF];						# start new run
			$lvcs = $run_start;
			$last_dcond = 0;
			next;													# process next scan
		}

		$ns++;
		for (my($j)=$run_start+1; $j<$i; $j++) {					# values return => remove bad spike
			$ants_[$j][$condF] = nan;
			$nv++;
		}
#		print(STDERR "values return at $i ($ants_[$i][$condF]) => deleted $run_start+1-$i; new run\n");
		$run_start = $i; $run_dcond = 0;
		$lvc = $ants_[$run_start][$condF];
		$lvcs = $run_start;
		$last_dcond = 0;
	}
	&antsAddParams("conductivity_spikes_removed",$ns);
	printf(STDERR "\n\t%d conductivity values removed from %d spikes",$nv,$ns) if ($opt_v > 1);

	printf(STDERR "\n") if ($opt_v);

} # if $opt_r

#----------------------------------------------------------------------
# Correcting for pressure bias
#----------------------------------------------------------------------

print(STDERR "Correcting for pressure bias...") if ($opt_v);

my($minP) = 9e99;
for (my($s)=0; $s<$nscans; $s++) {
	$minP = $ants_[$s][$pressF]
		if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP);
}
croak("$CNVfile: no valid CTD pressure data below 25dbar\n")
	unless ($minP < 9e99);

if ($minP < 25) {
	&antsAddParams('pressure_bias',$minP);
	printf(STDERR "\n\tsubtracting %.1f dbar",$minP) if ($opt_v > 1);
	for (my($s)=0; $s<$nscans; $s++) {
		$ants_[$s][$pressF] -= $minP
			if numberp($ants_[$s][$pressF]);
	}
} else {
	printf(STDERR "\n\tpartial-depth cast below %.1f dbar (no correction applied)",$minP) if ($opt_v > 1);
}
printf(STDERR "\n") if ($opt_v);

#----------------------------------------------------------------------
# Binning data
#----------------------------------------------------------------------

my($sps) = round(1 / $sampint / $opt_s);
print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v);
&antsAddParams('sampling_frequency',1/$opt_s);
&antsAddParams('sampling_rate',$opt_s);

my(@press,@temp,@cond);
my($sp,$np,$st,$nt,$sc,$nc);

$sp = $st = $sc = $np = $nt = $nc = 0;
$slat = $slon = $nGPS = 0;
for (my($rec)=1,my($s)=0; $s<$nscans; $s++) {
	if ($s*$sampint > $rec/$opt_s) {
		$rec++;
		push(@press,$np>0?$sp/$np:nan);
		push(@temp, $nt>0?$st/$nt:nan);
		push(@cond, $nc>0?$sc/$nc:nan);
		$sp = $st = $sc = $np = $nt = $nc = 0;
		if (defined($latF)) {
			push(@lat,$nGPS>0?$slat/$nGPS:nan);
			push(@lon,$nGPS>0?$slon/$nGPS:nan);
			$slat = $slon = $nGPS = 0;
		}
	}
	$sp+=$ants_[$s][$pressF],$np++ if numberp($ants_[$s][$pressF]);
	$st+=$ants_[$s][$tempF],$nt++ if numberp($ants_[$s][$tempF]);
	$sc+=$ants_[$s][$condF],$nc++ if numberp($ants_[$s][$condF]);
	if (defined($latF) && numberp($ants_[$s][$latF])) {
		$slat += $ants_[$s][$latF];
		$slon += $ants_[$s][$lonF];
		$nGPS++;
	}
}

printf(STDERR "\n") if ($opt_v);

#----------------------------------------------------------------------
# Calculating derived quantities
#----------------------------------------------------------------------

print(STDERR "Calculating vertical package velocity & sound speed...") if ($opt_v);

my($maxP) = -9e99; my($atBtm);
my($min_sspd) = 9e99; my($max_sspd) = -9e99;
for (my($r)=0; $r<@press; $r++) {
	$maxP=$press[$r],$atBtm=$r if ($press[$r] > $maxP);
	$elapsed[$r] = $r/$opt_s;
	$depth[$r] = &depth($press[$r],$lat);
	$salin[$r] = &salin($cond[$r],$temp[$r],$press[$r]);
	$sspd[$r]  = &sVel($salin[$r],$temp[$r],$press[$r]);
	$min_sspd = $sspd[$r] if ($sspd[$r] < $min_sspd);
	$max_sspd = $sspd[$r] if ($sspd[$r] > $max_sspd);
}

$w[0] = nan;
for (my($r)=1; $r<@depth-1; $r++) {
	$w[$r] = numbersp($depth[$r-1],$depth[$r+1])
		   ? ($depth[$r+1] - $depth[$r-1]) * $opt_s
		   : nan;
}
push(@w,nan);

printf(STDERR "\n") if ($opt_v);

#----------------------------------------------------------------------
# Low-pass filter velocity data
#	- interpolate missing vertical velocities first
#----------------------------------------------------------------------

if ($opt_l > 0) {
	print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v);
	&antsAddParams('w_lowpass_cutoff',$opt_l);
	
	my($trimmed) = 0;
	shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++
		until numberp($w[0]);
	my($interpolated) = 0;
	for ($r=1; $r<@w; $r++) {
		next if numberp($w[$r]);
		my($lv) = $r-1;
		for ($nv=$r+1; $nv<@depth && !numberp($w[$nv]); $nv++) {}
		if ($nv < @depth) {
			while ($r < $nv) {
				$w[$r] = $w[$lv] + ($r-$lv)/($nv-$lv) * ($w[$nv]-$w[$lv]);
				$interpolated++;
				$r++;
			}
		    
		} else {
			$trimmed += @w-$r;
			splice(@w,$r); splice(@depth,$r);
			splice(@elapsed,$r); splice(@sspd,$r);
		}
	}
	&antsAddParams('w_interpolated',$interpolated);
	printf(STDERR "\n\t%d/%d vertical velocities trimmed/interpolated",$trimmed,$interpolated) if ($opt_v > 1);


#--------------------
# Zero Pad Data
#--------------------
	
	for ($pot=1; $pot<@w; $pot<<=1) {}									# determine power of two
	
	for ($r=0; $r<@w; $r++) {											# copy data
		$fftbuf[2*$r] = $w[$r];
		$fftbuf[2*$r+1] = 0;
	}
	printf(STDERR "\n\t%d zero records added",$pot-$r) if ($opt_v > 1);
	while ($r < $pot) { 												# pad with zeroes
		$fftbuf[2*$r] = $fftbuf[2*$r+1] = 0;
		$r++;
	}
	
#--------------------
# Low-Pass Filter
#--------------------
	
	@fco = &FOUR1(-1,@fftbuf);											# forward FFT
	$n = @fco/2;
	for (my($ip)=2; $ip<=$n; $ip+=2) {									# +ve freq fco
		my($in) = 2*$n-$ip; 											# -ve freq fco
		my($f)	= $ip/2/$n*$opt_s; 										# frequency
		$fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0
			if ($f > 1/$opt_l); 										# low-pass filter
	}
	@w_lp = &FOUR1(1,@fco); 											# inverse FFT
	
	printf(STDERR "\n") if ($opt_v);
} else {
	@w_lp = @w;
}

#----------------------------------------------------------------------

if (defined($opt_p)) {
	print(STDERR "Plotting data...\n") if ($opt_v);
	my(@pfmt) = split(',',$opt_p);
	croak("$0: cannot decode -p $opt_p\n")
		unless (@pfmt == 2);

	my($xmin) = $elapsed[0]/60;
	my($xmax) = $elapsed[$#elapsed]/60;
	my($ymin) = -3; my($ymax) = 3;
	my($plotsize) = '13c';
	GMT_begin(sprintf($pfmt[0],$id),"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
	GMT_psxy('-W2,255/127/80');
	for ($r=0; $r<@w; $r++) {
		printf(GMT "%f %f\n",$elapsed[$r]/60,$w_lp[2*$r]/@w_lp);
		GMT_psxy('-W2,46/139/87') if ($r == $atBtm);
	}
	GMT_psbasemap('-Bg60a30f5:"Elapsed Time [min]":/g1a1f0.1:"Vertical Package Velocity [ms@+-1@+]":WeSn');
	GMT_unitcoords();
	GMT_pstext('-G255/127/80'); print(GMT "0.98 0.98 14 0 0 TR dc\n");
	GMT_pstext('-G46/139/87');  print(GMT "0.98 0.94 14 0 0 TR uc\n");
	GMT_pstext('-Gblue');
	if (defined($outfile)) { printf(GMT "0.02 0.02 14 0 0 BL $outfile\n",$id); }
	else 				   { printf(GMT "0.02 0.02 14 0 0 BL %03d\n",$id); }
	GMT_end();	                

	my($xmin) = round($min_sspd-3,5);
	my($xmax) = round($max_sspd+3,5);
	my($ymin) = 0; my($ymax) = round($depth[$atBtm]+70,100);
	my($plotsize) = '13c';
	GMT_begin(sprintf($pfmt[1],$id),"-JX${plotsize}/-${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
	GMT_psbasemap('-Bg10a10f1:"Speed of Sound [m/s]":/g1000a500f100:"Depth [m]":WeSn');
	GMT_psxy('-W8,255/127/80');
	for ($r=0; $r<@w; $r++) {
		printf(GMT "%f %f\n",$sspd[$r],$depth[$r]);
		GMT_psxy('-W6,46/139/87') if ($r == $atBtm);
	}
	GMT_unitcoords();
	GMT_pstext('-G255/127/80'); print(GMT "0.98 0.02 14 0 0 TR dc\n");
	GMT_pstext('-G46/139/87');  print(GMT "0.98 0.06 14 0 0 TR uc\n");
	GMT_pstext('-Gblue');
	if (defined($outfile)) { printf(GMT "0.02 0.98 14 0 0 BL $outfile\n",$id); }
	else 				   { printf(GMT "0.02 0.98 14 0 0 BL %03d\n",$id); }
	GMT_end();	                
}

#----------------------------------------------------------------------

print(STDERR "Writing output...\n") if ($opt_v);

@antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w.raw','w');
push(@antsNewLayout,'lat','lon') if (defined($latF));
for ($r=0; $r<@w; $r++) {
	&antsOut($elapsed[$r],$press[$r],$temp[$r],$cond[$r],$depth[$r],$salin[$r],
			 $sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp,$lat[$r],$lon[$r]);
}

exit(0);															# don't flush @ants_