#!/usr/bin/perl

#======================================================================

#                    M E R G E C T D + L A D C P 

#                    doc: Wed Apr 28 21:50:01 2010

#                    dlm: Fri Dec 17 16:51:33 2010

#                    (c) 2010 A.M. Thurnherr

#                    uE-Info: 25 44 NIL 0 0 72 2 2 4 NIL ofnI

#======================================================================

# TODO:

#	- make sure range data in ADCP header are really along-beam

$antsSummary = 'match and merge LADCP with CTD time series for caclulating vertical velocities';

# HISTORY:

#	Apr 28, 2010: - created early version

#	May 12, 2010: - time-shift each ensemble separately (V2)

#	May 21, 2010: - read raw LADCP files (V3)

#	Jun  5, 2010: - first "fully debugged" version (V4)

#	Oct 29, 2010: - hastily adapted to Earth coordinates

#	Oct 30: 2010: - -w => -l, max -g)ap removed, added -w)ater depth

#	Dec 14, 2010: - BUG: percentages were wrong (forgot to round)

#	Dec 16, 2010: - removed hard-coded PERL_TOOLS directory

#				  - added -m)ax vertical velocity

#	Dec 17, 2010: - added -a)djust CTD depth

# ASSUMPTIONS:

#	- LADCP elapsed field lagged approximately at beginning of record 

#	- LADCP data collected with nominal soundspeed of 1500m/s

#	- correct layout of old time-series file

# NOTES:

#	- by default, bin 1 is not used

#	- the initial lag estimate can be calculated by creating a 1Hz version

#	  of the CTD file; the number of seconds that are added to the CTD elapsed

#	  field to achieve a visually satisfying match is the initial guess for the lag

#	- if there are a sequence of +ve jumps at the beginning of the profile, the

#	  initial lag estimate should be decreased

#	- in a clean data set, decreasing -m should not be required; however,

#	  in some cases (ambiguity velocity too small?) there are large LADCP

#	  ref-lr w spikes that must be ignored; this can be accomplished by

#	  setting -m to a value like 1m/s; the symptom of the problem is usually

#	  that the time-lagging does not work

# POTENTIAL IMPROVEMENTS:

#	- correct for sound-speed variations within instrument range

#		- probably not important except perhaps near sea surface

#		- cannot be done for data below max CTD depth

#	- allow specifying vertical offset between CTD & ADCP

#		- not important, given size of output bins

#	- add edits for vertical & horizontal rotation velocities

#		- may be important, as it is implemented in veloity inversion

#	- add user-definable edit functions

#	- estimate first guess of CTD_timoff

#	- output correct ensemble number during warnings

($ANTS) 	  = ($0 =~ m{^(.*)/[^/]*$});

($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});

require "$ANTS/ants.pl";

require "$ANTS/libstats.pl";

require "$ANTS/mergeCTD+LADCP.bestLag";

require "$ANTS/mergeCTD+LADCP.editData";

require "$ANTS/mergeCTD+LADCP.timeMatchingProblems";

require "$ANTS/mergeCTD+LADCP.LADCPtimeSeries";

require "$ANTS/mergeCTD+LADCP.seabed";

require "$PERL_TOOLS/RDI_BB_Read.pl";

require "$PERL_TOOLS/RDI_Coords.pl";

$antsParseHeader = 0;

&antsUsage('4a:b:c:d:e:fl:m:o:r:w:p:s:t:z:',2,

	'[-w)ater <depth>]',

	'[-l)ag-search <window[60s]>]',

	'[max lag-s)earch <time[2s]> between LADCP ensembles]',

	'[ref-l)ayer <bin[2],bin[6]>] [valid -b)ins for final w <bin[2],bin[*]>',

	'[-c)orrelation <min[70]>] [-p)itch/roll <max[10deg]>',

	'[-e)rr-vel <max[0.1m/s]>]',

	'[-m)ax vertical <velocity[10m/s]>',

    '[require -4)-beam solutions]',

	'[read/write -t)imes-eries <file> [-f)orce re-make]]',

	'[-o)uput bin <resolution[10m]>',

	'[-d)ump depth-bins to <basename>]',

	'[-a)djust CTD depth <by>] [CTD time re-z) <val[1/24s]>]',

	'[--] <initial CTD time offset> <LADCP-file> [CTD-file]');

&antsFloatOpt(\$opt_a,0);					# CTD depth adjustment

&antsFloatOpt(\$opt_z,1/24);				# CTD time resolution

&antsFloatOpt(\$opt_l,60);					# time-lagging window size

&antsFloatOpt(\$opt_s,2);					# max possible CTD time jump btw LADCP pings

&antsFloatOpt(\$opt_c,70);					# min correlation

&antsFloatOpt(\$opt_p,5);					# max pitch/roll

&antsFloatOpt(\$opt_e,0.1);					# max err vel

&antsFloatOpt(\$opt_m,1.0);					# max allowed vertical velocity

&antsFloatOpt(\$opt_o,10);					# max err vel

$opt_b = '2,*' unless defined($opt_b);		# bins to use

$opt_r = '2,6' unless defined($opt_r);		# reference layer bins for w for time matching

$RDI_Coords::minValidVels = 4 if ($opt_4);	# suppress 3-beam solutions

($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b);				# decode -b

croak("$0: cannot decode -b $opt_b\n")

	unless (numberp($LADCP_firstBin) &&

			($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));

($refLr_firstBin,$refLr_lastBin) = split(',',$opt_r);				# decode -r

croak("$0: cannot decode -r $opt_r\n")

	unless numberp($refLr_firstBin) && numberp($refLr_lastBin);

$CTD_timoff = &antsFloatArg();										# initial (good) guess 

$LADCP_file = &antsFileArg();

{ my($warned) = 0;													# handle warnings

  sub warning(@)

  {

	  print(STDERR "\n");

	  &antsInfo(@_);

	  $warned++;

	  croak("$0: too many warnings\n")

	  	if ($warned > 10);

  }

}

#======================================================================

# STAGE 1: match CTD to LADCP time series, find seabed, edit sidelobes

#

# KEY VARIABLES:

#	$LADCP_bufStart		index of first element in LADCP "buffer" [0]

#	$ens				index of current target elapsed time of LADCP

#	$LADCP_bufEnd		index of first element in LADCP "buffer" [$#LADCP]

#	$CTD_trg			index of record closest to target elapsed time in CTD buf

#	$CTD_timoff			current best CTD time offset

#======================================================================

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

# S1 STEP: Read LADCP data

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

print(STDERR "Reading LADCP data ($LADCP_file)...");

readData($LADCP_file,\%LADCP);

printf(STDERR "\n\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));

croak("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -l\n")

	unless ($LADCP{N_BINS} >= $refLr_lastBin);

#croak("$LADCP_file: require beam coordinates\n")

#	unless ($LADCP{BEAM_COORDINATES});

croak("$0: first reference-layer bin outside valid range\n")

	unless ($refLr_firstBin>=1 && $refLr_firstBin<=$LADCP{N_BINS});

croak("$0: last reference-layer bin outside valid range\n")

	unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS});

croak("$0: first reference-layer bin > last reference-layer bin\n")

	unless ($refLr_firstBin <= $refLr_lastBin);

$LADCP_lastBin = $LADCP{N_BINS}-1

	if ($LADCP_lastBin eq '*');

croak("$0: first valid LADCP bin outside valid range\n")

	unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS});

croak("$0: last valid LADCP bin outside valid range\n")

	unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS});

croak("$0: first valid LADCP bin > last valid LADCP bin\n")

	unless ($LADCP_firstBin <= $LADCP_lastBin);

warning("WARNING: first reference-layer bin < first valid LADCP bin")

	unless ($refLr_firstBin >= $LADCP_firstBin);

warning("WARNING: last reference-layer bin > last valid LADCP bin")

	unless ($refLr_lastBin <= $LADCP_lastBin);

warning("WARNING: if at all, bin 1 should not be used for short blank-after-transmit values")

	if ($LADCP{BLANKING_DISTANCE} < $LADCP{BIN_LENGTH} &&

		($LADCP_firstBin == 1 || $refLr_firstBin == 1));

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

# S1 STEP: Edit beam-velocity data

#	1) correlation threshold

#	2) pitch/roll threshold

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

if ($LADCP{BEAM_COORDINATES}) {

	print(STDERR "Editing beam-velocity data...");

	$nvv = $cte = 0;

	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

		$nvv += countValidBeamVels($ens);

		$cte += editCorr($ens,$opt_c);

		$pte += editPitchRoll($ens,$opt_p);

	}

	print(STDERR "\n");

	printf(STDERR "\tcorrelation threshold (-c %d): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));

	printf(STDERR "\tattitude threshold (-p %d): %d velocites removed (%d%% of total)\n",$opt_p,$pte,round(100*$pte/$nvv));

} else {

	print(STDERR "Editing velocity data...");

	$nvv = $cte = 0;

	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

		$nvv += countValidBeamVels($ens);

		$cte += editCorr_Earthcoords($ens,$opt_c);

		$pte += editPitchRoll($ens,$opt_p);

	}

	print(STDERR "\n");

	printf(STDERR "\tcorrelation threshold (-c %d): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));

	printf(STDERR "\tattitude threshold (-p %d): %d velocites removed (%d%% of total)\n",$opt_p,$pte,round(100*$pte/$nvv));

}

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

# S1 STEP: Calculate earth velocities

#	- this is done for all bins (not just valid ones), to allow

#	  useless possibility that invalid bins are used for reflr calcs

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

if ($LADCP{BEAM_COORDINATES}) {

	print(STDERR "Calculating earth-coordinate velocities...");

	$nvw = 0;

	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

		for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

			($LADCP{ENSEMBLE}[$ens]->{U}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{V}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{W}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) =

				velInstrumentToEarth(\%LADCP,$ens,

					velBeamToInstrument(\%LADCP,

										@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));

			$nvw += defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);

		}

	}

	print(STDERR "\n\t$nvw valid velocities\n");

	printf(STDERR "\t3-beam solutions : $RDI_Coords::threeBeam_1 " .

									   "$RDI_Coords::threeBeam_2 " .

									   "$RDI_Coords::threeBeam_3 " .

									   "$RDI_Coords::threeBeam_4\n")

	    unless ($opt_4);

} else {

	print(STDERR "Counting valid vertical velocities...");

	$nvw = 0;

	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

		for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

			($LADCP{ENSEMBLE}[$ens]->{U}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{V}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{W}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) = @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]};

			$nvw += defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);

		}

	}

	print(STDERR "\n\t$nvw valid velocities\n");

}

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

# S1 STEP: Edit earth-coordinate -velocity data

#	1) error-velocity threshold

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

print(STDERR "Editing earth-coordinate velocity data...");

$ete = 0;

for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

	$ete += editErrVel($ens,$opt_e);

}

printf(STDERR "\n\t error-velocity threshold (-e %.1f): %d velocites removed (%d%% of total)\n",$opt_e,$ete,round(100*$ete/$nvw));

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

# S1 BRANCH POINT: If up-to-date time-series file is found, read it

#				   and skip to STAGE 2

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

if ($opt_t && -f $opt_t) {

	print(STDERR "Time-series file <$opt_t> found:\n");

	my(@s,$ct,$Lct,$Cct);

	my($ctf) = 10;

	die unless (@s = stat($opt_t)); $ct = $s[$ctf];				# ctime of time-series file

	die unless (@s = stat($LADCP_file)); my($Lct) = $s[$ctf]; 	# ctime of LADCP file

	printf(STDERR "\t%s wrt $LADCP_file\n",$Lct>$ct ? 'stale' : 'up-to-date');

	if ($ARGV[0] eq '-') {

		print(STDERR "\tCTD data on stdin assumed up-to-date\n");

	} else {

		die("$ARGV[0]: $!") unless (@s = stat($ARGV[0])); $Cct = $s[$ctf];			# ctime of CTD file

	    printf(STDERR "\t%s wrt $ARGV[0]\n",$Cct>$ct ? 'stale' : 'up-to-date');

	}

	if ($ct >= $LCt && $ct >= $Cct && !$opt_f) {				# read time-series file

		print(STDERR "\t-> reading time-series data from <$opt_t>...");

		croak("$0: need CTD data\n")							# read 1st CTD rec

			unless (&antsIn());									# to get %PARAMs

		open(F,$opt_t) || croak("$opt_t: $!\n");

		$firstGoodEns 	= &antsFileScanParam(F,'LADCP_prof_start');

		$LADCP_atbottom = &antsFileScanParam(F,'LADCP_at_bottom');

		$lastGoodEns 	= &antsFileScanParam(F,'LADCP_prof_end');

		$water_depth 	= &antsFileScanParam(F,'water_depth');

		$sig_water_depth= &antsFileScanParam(F,'water_depth.sig');

		&antsAddParams('rms_w_reflr_err',&antsFileScanParam(F,'rms_w_reflr_err'),

					   'rms_w_reflr_err_interior',&antsFileScanParam(F,'rms_w_reflr_err_interior'),

					   'LADCP_prof_start',&antsFileScanParam(F,'LADCP_prof_start'),

					   'LADCP_at_bottom',&antsFileScanParam(F,'LADCP_at_bottom'),

					   'LADCP_prof_end',&antsFileScanParam(F,'LADCP_prof_end'));

		for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {

			my(@r) = &antsFileIn(F);

			die("$0: assertion failed") unless(@r);

			my($f) = 0;

			croak("Incompatible time-series file (expected ens $ens, got $r[$f])?\n")

				unless ($r[$f] == $ens);

			$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}      = $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{ELAPSED_MISMATCH} 	= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{CTD_TIMOFF} 		= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{PAST_LAG}			= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{PAST_ERR}			= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{FUTURE_LAG}		= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{FUTURE_ERR}		= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{DEPTH}				= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{SOUND_SPEED}		= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{PITCH}				= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{ROLL}				= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{HEADING}			= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{CTD_W}				= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{CTD_W_T}			= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{REFLR_W}			= $r[++$f];

			$LADCP{ENSEMBLE}[$ens]->{REFLR_W_ERR}		= $r[++$f];

		}

		close(F);

		print(STDERR "done\n");

		goto STAGE_2;

	}

}

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

# S1 STEP: Calculated LADCP time series

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

print(STDERR "Calculating LADCP time-series...");

($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) =

	calcLADCPts(\%LADCP,$refLr_firstBin,$refLr_lastBin,$opt_g);

croak("$LADCP_file: no good ensembles\n")

	unless defined($firstGoodEns);

printf(STDERR "\n\tStart of cast : %s (#%5d)\n",

					$LADCP{ENSEMBLE}[$firstGoodEns]->{TIME},

					$LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER});

printf(STDERR "\tBottom of cast: %s (#%5d)\n",

					$LADCP{ENSEMBLE}[$LADCP_atbottom]->{TIME},

					$LADCP{ENSEMBLE}[$LADCP_atbottom]->{NUMBER});

printf(STDERR "\tEnd of cast   : %s (#%5d)\n",

					$LADCP{ENSEMBLE}[$lastGoodEns]->{TIME},

					$LADCP{ENSEMBLE}[$lastGoodEns]->{NUMBER});

$cast_duration = int($LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED_TIME} -

				     $LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED_TIME});

printf(STDERR "\tCast Duration : %.1f hours (pinging for %.1f hours)\n",

					$cast_duration / 3600,

					($LADCP{ENSEMBLE}[$#{$LADCP{ENSEMBLE}}]->{UNIX_TIME} -

						$LADCP{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);

croak("$0: no valid LADCP data\n")

	unless ($lastGoodEns-$firstGoodEns > 0);

croak("$0: implausibly short cast ($cast_duration seconds)\n")

	unless ($cast_duration > 600);

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

# S1 STEP: time-lag CTD time series to each ensemble

#

# NOTES:

#	- time lag done in window of size (in seconds) $opt_l

#	- after this step, reflr w is sound-speed corrected!!!

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

print(STDERR "Matching CTD to LADCP data");

$LADCP_bufStart = $ens = $LADCP_bufEnd = $firstGoodEns;					# init pointers

sub growLADCPbuf()														# adjust $LADCP_bufEnd

{																		# can safely be called often :-)

	$LADCP_bufEnd++ while ($LADCP_bufEnd<=$lastGoodEns && &Le($LADCP_bufEnd)-&Le($ens) < $opt_l);

	return $LADCP_bufEnd<=$lastGoodEns;

}

$sumWsq = $nWsq = $sumWsqI= $nWsqI = 0;									# for rms(w_ocean_ref_lr)

$antsBufSkip = 0;														# don't remove recs from buf

croak("$0: need CTD data\n")

	unless (&antsIn());													# read 1st CTD rec

($CTD_elapsed,$CTD_depth,$CTD_w,$CTD_w_t,$CTD_svel) =					# decode CTD Layout

	&fnr('elapsed','depth','w','w_t','ss');								

&antsInstallBufFull('&Ce($#ants_)-&Le($ens) >= $opt_l+$opt_s');			# CTD buffer full criterion

&antsIn();																# fill CTD buffer for 1st LADCP ens

&growLADCPbuf();														# set up LADCP buffer pointers

if (&Ce(0) > &Le($LADCP_bufStart)) {									# CTD data begin with CTD in water

	print(STDERR "\n");

	&antsInfo("CTD data begin with LADCP in water (Ce = %f, Le = %f)",&Ce(0),&Le($LADCP_bufStart));

	$ens++,$LADCP_bufStart++ while (&Ce(0) > &Le($LADCP_bufStart));

	&growLADCPbuf();

	&antsIn();

}

die(sprintf("$0: assertion failed (ens=$ens, Ce($#ants_) = %f, Le($ens) = %f)",&Ce($#ants_),&Le($ens)))

	unless (&Ce($#ants_) >= &Le($ens));

$CTD_trg = $#ants_ - int(($opt_l+$opt_s)/$opt_z);						# nearest CTD record is -w + -s before end

$CTD_trg++ if (abs(&Ce($CTD_trg)-&Le($ens)) > abs(&Ce($CTD_trg+1)-&Le($ens))); # determine nearer of two bracketing recs

$CTD_trg-- if (abs(&Ce($CTD_trg)-&Le($ens)) > abs(&Ce($CTD_trg-1)-&Le($ens))); # can happen when CTD data begin with LADCP in water

while ($ens <= $lastGoodEns && &growLADCPbuf()) {						# 2nd test required on EOF(CTD)

#	printf(STDERR "target = %f\n",&Le($ens));

	die("$0: assertion failed (CTD_trg=$CTD_trg, ens=$ens)") unless ($CTD_trg>=0 && $CTD_trg<=$#ants_);

	die("$0: assertion failed (ens=$ens)") unless ($ens>=0 && $ens<=$LADCP_bufEnd);

	die(sprintf("$0: assertion failed (ens=$ens, abs(%f-%f) > $opt_z/2)",&Ce($CTD_trg-1),&Le($ens)))

		unless (abs(&Ce($CTD_trg)-&Le($ens)) <= $opt_z/2);

	$LADCP_bufStart++ while (&Le($ens)-&Le($LADCP_bufStart+1) >= $opt_l);		# shift beginning of LADCP buffer

	die("$0: assertion failed (ens=$ens)") unless ($ens>=$LADCP_bufStart);

	die("$0: assertion failed (ens=$ens)") unless (abs(&Ce($CTD_trg)-&Le($ens)) <= $opt_z/2);

	if (&Le($ens)-&Le($LADCP_bufStart) >= $opt_l) {	 					# find best lag using past data

		($pCts,$pErr) = &bestLag($LADCP_bufStart,$ens,$ens,$CTD_trg);

		$pLag = $pCts - $CTD_trg;

	} else {

		$pCts = $pErr = $pLag = undef;

	}

	if (&Le($LADCP_bufEnd)-&Le($ens) >= $opt_l &&	 					# find best lag using future data

		&Ce($#ants_)-&Le($ens) >= $opt_l+$opt_s) {

		($fCts,$fErr) = &bestLag($ens,$LADCP_bufEnd,$ens,$CTD_trg);

		$fLag = $fCts - $CTD_trg;

	} else {

		$fCts = $fErr = $fLag = undef;

	}

	if ((defined($fErr) && defined($pErr) && $pErr <= $fErr) ||			# past data are better correlated

	    (defined($pErr) && !defined($fErr))) {

		    $CTD_timoff -= ($pCts-$CTD_trg) * $opt_z;

		    $CTD_trg = $pCts;

    } elsif ((defined($fErr) && defined($pErr) && $fErr <= $pErr) ||	# return w based on future lagcorr

			 (defined($fErr) && !defined($pErr))) {

		    $CTD_timoff -= ($fCts-$CTD_trg) * $opt_z;

		    $CTD_trg = $fCts;

	}

	die("$0: assertion failed (ens=$ens)") unless ($CTD_trg>=0 && $CTD_trg<=$#ants_);

	die("$0: assertion failed (ens=$ens)") unless (abs(&Ce($CTD_trg)-&Le($ens)) <= $opt_z/2);

	croak("$0: CTD velocity spike >5m/s at elapsed time = $ants_[$CTD_trg][$CTD_elapsed]\n")

		if (abs($ants_[$CTD_trg][$CTD_w] > 5));

	&detectTimeMatchingProblems($ens,$ants_[$CTD_trg][$CTD_depth]+$opt_a,$oCo,$CTD_timoff)

		if defined($oCo);

	$oCo = $CTD_timoff;

	if (defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W}) &&							# not a gap

		(defined($fErr) || defined($pErr))) {									# at least one correlation is defined

	    $LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $ants_[$CTD_trg][$CTD_svel]/1500; 	# correct for sound-speed variations

		$LADCP{ENSEMBLE}[$ens]->{ELAPSED_MISMATCH} = &Ce($CTD_trg) - &Le($ens);

		$LADCP{ENSEMBLE}[$ens]->{CTD_TIMOFF} = $CTD_timoff;

		$LADCP{ENSEMBLE}[$ens]->{PAST_LAG} = $pLag;

		$LADCP{ENSEMBLE}[$ens]->{PAST_ERR} = $pErr;

		$LADCP{ENSEMBLE}[$ens]->{FUTURE_LAG} = $fLag;

		$LADCP{ENSEMBLE}[$ens]->{FUTURE_ERR} = $fErr;

		croak(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed = %.1fs (use -a?)\n",

			$ants_[$CTD_trg][$CTD_depth]+$opt_a,$ants_[$CTD_trg][$CTD_elapsed]))

			unless ($ants_[$CTD_trg][$CTD_depth]+$opt_a >= 0);

		$LADCP{ENSEMBLE}[$ens]->{DEPTH} = $ants_[$CTD_trg][$CTD_depth]+$opt_a;

		$LADCP{ENSEMBLE}[$ens]->{SOUND_SPEED} = $ants_[$CTD_trg][$CTD_svel];

		$LADCP{ENSEMBLE}[$ens]->{CTD_W} = $ants_[$CTD_trg][$CTD_w];

		$LADCP{ENSEMBLE}[$ens]->{CTD_W_T} = $ants_[$CTD_trg][$CTD_w_t];

		my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $LADCP{ENSEMBLE}[$ens]->{CTD_W};

		if (abs($reflr_ocean_w) <= $opt_m) {

			$sumWsq += &SQR($reflr_ocean_w);

			$nWsq++;

			if ($LADCP{ENSEMBLE}[$ens]->{DEPTH} > 100 &&

				$LADCP{ENSEMBLE}[$ens]->{DEPTH} < $LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH}-100) {

					$sumWsqI += &SQR($reflr_ocean_w);

					$nWsqI++;

	        }

	    } else {

	    	undef($LADCP{ENSEMBLE}[$ens]->{DEPTH});						# DON'T USE THIS ENSEMBLE LATER

	    }

	} else{

	    undef($LADCP{ENSEMBLE}[$ens]->{REFLR_W});						# don't output in time-series file

	    undef($LADCP{ENSEMBLE}[$ens]->{DEPTH});							# old DEPTH from calcLADCPts()

	}

	$ens++;																# update target indices

	$CTD_trg += int((&Le($ens)-&Le($ens-1)) / $opt_z);					# NB: int() rounds down

	unless ($CTD_trg <= $#ants_) {

		print(STDERR "\n");

		&antsInfo("CTD data end with LADCP in water");

		last;

	}

	$CTD_trg++ if (abs(&Le($ens)-&Ce($CTD_trg+1)) < abs(&Le($ens)-&Ce($CTD_trg)));

	die("$0: assertion failed (ens=$ens)") unless ($CTD_trg>=0 && $CTD_trg<=$#ants_);

	die("$0: assertion failed (ens=$ens)") unless (abs(&Ce($CTD_trg)-&Le($ens)) <= $opt_z/2);

	while (&Le($ens)-&Ce(1) >= $opt_l+$opt_s) {							# remove no-longer-needed CTD records

		splice(@ants_,0,1);

		$CTD_trg--;

	}

	die("$0: assertion failed (ens=$ens)") unless ($CTD_trg>=0 && $CTD_trg<=$#ants_);

	die("$0: assertion failed (ens=$ens)") unless (abs(&Ce($CTD_trg)-&Le($ens)) <= $opt_z/2);

	$eof = !&antsIn() unless ($eof);									# don't antsIn after EOF

    print(STDERR '.')

		if ($ens % 100 == 0);

}

if ($nWsq > 0 && $nWsqI > 0) {

	&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI));

	printf(STDERR "\n\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n",

						100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI));

	croak("$0: rms reference-layer w_ocean is too large\n")

		unless (sqrt($sumWsqI/$nWsqI) < 0.05);

} elsif ($nWsq > 0) {

	&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',nan);

	printf(STDERR "\n\t%.2f cm/s rms reference-layer w_ocean\n",100*sqrt($sumWsq/$nWsq));

} else {

	croak("$0: no valid vertical velocities\n");

}

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

# S1 STEP: Find seabed

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

if (defined($opt_w)) {

	print(STDERR "Using specified water depth...");

	$water_depth = $opt_w; $sig_water_depth = nan;