#!/usr/bin/perl

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

#                    L A D C P _ W 

#                    doc: Fri Dec 17 18:11:13 2010

#                    dlm: Sun Dec 26 20:14:32 2010

#                    (c) 2010 A.M. Thurnherr

#                    uE-Info: 34 86 NIL 0 0 72 2 2 4 NIL ofnI

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

# TODO:

#	make timelagging work for short casts (make sure 10% is not more than half window size)

#	own seabed detection (P403)

#	instrument tilt in sidelobe editing?

#	read ADCP soundspeed from data file

$antsSummary = 'calculate vertical velocities from LADCP & CTD time series';

# HISTORY:

#	Dec 17, 2010: - created from [mergeCTD+LADCP]

#	Dec 18, 2010: - made to work

#	Dec 19, 2010: - improved considerably

#	Dec 20, 2010: - onward

#				  - BUG: depth-binning was off by 1 bin?!

#				  - added binning correction for instrument tilt

#	Dec 21, 2010: - added -h (seafloor depth)

#	Dec 22, 2010: - BUG: had not applied soundspeed-correction to w

#				  - debugged opt_d

#	Dec 23, 2010: - continued implementation of soundspeed corrections

#	Dec 24, 2010: - added winch_w, wave_w

#				  - removed beampair velocities from code

#	Dec 25, 2010: - adapted for surface-wave correction in terms of acceleration (CTD_w_t)

#				  - removed elapsed_mismatch

#				  - removed winch_w, wave_w

#	Dec 26, 2010: - made -p output layout independent of -x to avoid Makefile problems

# CTD REQUIREMENTS

#	- elapsed		elapsed seconds; see note below

#	- depth

#	- ss			sound speed

#	- w				ddepth/dt

#	- w_t			dw/dt

# NUMERICAL OPTIONS

#	- the first option in the list cannot be numerical!

#	- if need be, use -v 1 as a dummy option

# ELAPSED TIMES

#	- there are 2 different elapsed times used in this program:

#		1) elapsed based on firstgoodens in the LADCP time series

#		2) CTD elapsed time

#	- CTD elapsed time does not have to start with zero!

#	- do not use the Seabird elapsed field, which is only reported to

#	  3 significant digits, causing significant jitter in dt; however,

#	  at least up to 2010 Seabird simply calculates elapsed time by 

#	  assuming a 24Hz sampling rate and no record drop; therefore,

#	  it is best to calculate elapsed time as %RECNO/24

#	- the elapsed field of the output is the elapsed time from the CTD

#	  file; this is required in order to be able to compare the times

#	  from the uplooker and downlooker-derived vertical velocity

#	  profiles

#	- as a result, a profile only starts with elapsed==0 if the CTD

#	  is turned on when the LADCP is already in the water

# TIME LAGGING

#	- occasionally, the time lagging algorithm fails, in particular if the CTD is turned on

#	  some time after the package enters the water

#	- in this case, an initial guess can be provided with -i

#	- e.g. plot 'CTD/24Hz/054.1Hz elapsed w','LADCP/raw/054UL.prof =$elapsed-850 w' => -i -850

# VELOCITY AMBIGUITY ERRORS

#	- quite extensive tests with DIMES US2 station 146, which has a lot of

#	  ambiguity velocity errors, reveal that -m 1 catches those errors

#	  quite nicely

#	- even when the errors are not filtered with -m 1, they do not

#	  affect the w profiles, as long as the median bin values are used

# SCREEN LOGGING

#	- there are 4 verbosity levels, selected by -v

#		0:	only print errors

#		1:	default, UNIX-like (warnings and info messages that are not produced for every cast)

#		2:	progress messages and useful information

#		>2:	debug messges

#	- the most useful ones of these are 1 & 2

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

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

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

require "$ANTS/ants.pl";

require "$ANTS/libstats.pl";

require "$WCALC/edit_LADCP_data.pl";

require "$WCALC/calc_LADCP_tis.pl";

require "$WCALC/calc_timelags.pl";

require "$WCALC/find_seabed.pl";

require "$PERL_TOOLS/RDI_BB_Read.pl";

require "$PERL_TOOLS/RDI_Coords.pl";

$antsParseHeader = 0;

&antsUsage('3:4a:b:c:d:e:f:g:h:i:m:n:o:p:r:s:t:v:w:x:',1,

	'[-v)erbosity <level[1]>]',

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

	'[-r)ef-layer <bin[2],bin[6]>]',

	'[-c)orrelation <min[70]>] [-t)ilt <max[10deg]> [-e)rr-vel <max[0.1m/s]>]',

	'[-h water <depth>]',

	'[max LADCP time-series -g)ap <length[60s]>]',

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

	'[-a)djust CTD depth <by[0m]>]',

	'[-i)nitial CTD time offset <guestimate>]',

	'[calculate -n) <lags,lags[10,100]>] [lag -w)indow <sz,sz[240s,20s]>]',

	'[require top-3) lags to account for <frac[0.6]> of all]',

	'[-x <dc_correction_params/uc_correction_param>]',

	'[use LADCP -b)ins for w profiles <bin[2],bin[*]>',

	'[-o)utput bin <resolution[10m]>] [require -s)amples <min[20]>]',

	'[-f) write time-series <file>] [-d)ump depth-bins to <basename>] [-p)ackage effect <file>]',

	'<LADCP-file> [CTD-file]');

&antsCardOpt(\$opt_v,1);											# suppress regular info

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

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

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

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

&antsFloatOpt(\$opt_g,60);											# max LADCP gap length

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

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

$opt_n = '10,100' unless defined($opt_n);							# number of time-lags to carry out

$opt_w = '240,20' unless defined($opt_w);							# time-lag search window (full width)

&antsFloatOpt(\$opt_3,0.6);

&antsFloatOpt(\$opt_o,50);											# output bin size

&antsCardOpt(\$opt_s,20);											# min samples required

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

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

@n_lags = split(',',$opt_n);										# decode -n

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

	unless numberp($n_lags[0]) && numberp($n_lags[1]);

@w_size = split(',',$opt_w);										# decode -w

croak("$0: cannot decode -w $opt_w\n")

	unless numberp($w_size[0]) && numberp($w_size[1]);

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

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

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

($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)));

if (defined($opt_x)) {												# decode corrections

	my($dccps,$uccps) = split('/',$opt_x);

	(@dc_corr_poly) = split(',',$dccps);

	(@uc_corr_poly) = split(',',$uccps);

	croak("$0: cannot decode -x $opt_x\n")

		unless @dc_corr_poly>0 && @uc_corr_poly>0;

	&antsAddParams('dc_corr_intercept',$dc_corr_poly[0],'dc_corr_slope',$dc_corr_poly[1]);

	&antsAddParams('uc_corr_intercept',$uc_corr_poly[0],'uc_corr_slope',$uc_corr_poly[1]);

}

if (defined($opt_d)) {												# make sure output directory is clean

	croak("$0: old depth-bin files <${opt_d}[0-9][0-9][0-9].dncast> found --- remove before creating new ones!\n")

		if (glob("${opt_d}[0-9][0-9][0-9].dncast"));

	croak("$0: old depth-bin files <${opt_d}[0-9][0-9][0-9].upcast> found --- remove before creating new ones!\n")

		if (glob("${opt_d}[0-9][0-9][0-9].upcast"));

}

$LADCP_file = &antsFileArg();

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

# Screen Logging

#	- warning levels:

#		0 probably unimportant, e.g. nonsensical parameters that probably won't affect solution

#		1 may be somewhat important

#		2 important

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

sub progress(@)

{  printf(STDERR @_) if ($opt_v > 1); }

sub info(@)

{

	print(STDERR ($opt_v > 1) ? "\t" : "$LADCP_file: ");

	printf(STDERR @_)  if ($opt_v > 0);

}

sub warning(@)

{

	my($lvl,@msg) = @_;

	return if ($opt_v == 0);

	print(STDERR "\n") if ($opt_v > 1);

	print(STDERR "WARNING (level $lvl): ");

	printf(STDERR @msg);

	print(STDERR "\n") if ($opt_v > 1);

}

sub debugmsg(@)

{ printf(STDERR @_) if ($opt_v > 2); }

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

# Read LADCP data

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

progress("Reading LADCP data ($LADCP_file)...\n");

readData($LADCP_file,\%LADCP);

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

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

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

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(0,"first reference-layer bin < first valid LADCP bin\n")

	unless ($refLr_firstBin >= $LADCP_firstBin);

warning(0,"last reference-layer bin > last valid LADCP bin\n")

	unless ($refLr_lastBin <= $LADCP_lastBin);

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

	if ($LADCP{BLANKING_DISTANCE}<$LADCP{BIN_LENGTH} && $refLr_firstBin==1);

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

# Edit beam-velocity data

#	1) correlation threshold

#	2) tilt threshold (also sets TILT field in all ensembles)

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

if ($LADCP{BEAM_COORDINATES}) {

	progress("Editing beam-velocity data...\n");

	$nvv = $cte = 0;

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

		$nvv += countValidBeamVels($ens);

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

		$pte += editTilt($ens,$opt_t);

	}

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

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

} else {

	progress("Editing velocity data...\n");

	$nvv = $cte = 0;

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

		$nvv += countValidBeamVels($ens);

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

		$pte += editTilt($ens,$opt_t);

	}

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

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

}

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

# 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

#	- also calculate separate beam-pair velocities

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

if ($LADCP{BEAM_COORDINATES}) {

	progress("Calculating earth-coordinate velocities...\n");

	$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]);

			$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =

			 	gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});

		}

	}

	progress("\t$nvw valid velocities\n");

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

								  "$RDI_Coords::threeBeam_2 " .

								  "$RDI_Coords::threeBeam_3 " .

								  "$RDI_Coords::threeBeam_4\n")

	    unless ($opt_4);

} else {

	progress("Counting valid vertical velocities...\n");

	$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]);

		}

	}

	progress("\t$nvw valid velocities\n");

}

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

# S1 STEP: Edit earth-coordinate -velocity data

#	1) error-velocity threshold

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

progress("Editing earth-coordinate velocity data...\n");

$ete = 0;

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

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

}

progress("\t error-velocity threshold (-e %.1f): %d velocites removed (%d%% of total)\n",

	$opt_e,$ete,round(100*$ete/$nvw));

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

# Calculate LADCP time series

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

progress("Calculating LADCP time-series...\n");

($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) && ($lastGoodEns-$firstGoodEns > 0);

my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} -

				     $LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED};

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

	unless ($cast_duration > 600);

$LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1);

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

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

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

progress("\tBottom of cast    : %s (#%5d) @ dz~%.1fm\n",

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

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

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

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

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

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

progress("\tCast duration     : %.1f hours (pinging for %.1f hours)\n",

					$cast_duration / 3600,

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

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

progress("\tMean ping interval: %.1f seconds\n",$LADCP{MEAN_DT});

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

# Read CTD data

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

progress("Reading CTD data...\n");

croak("$0: no CTD data\n") unless (&antsIn());

($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w,$CTD_w_t) =

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

$CTD_maxdepth = -1;

do {

	croak("$0: cannot deal with non-numeric CTD fields\n")

		unless &antsNumbers($CTD_elapsed,$CTD_depth,$CTD_w);

	push(@{$CTD{ELAPSED}},$ants_[0][$CTD_elapsed]);

	push(@{$CTD{DEPTH}},  $ants_[0][$CTD_depth]+$opt_a);

	push(@{$CTD{SVEL}},   $ants_[0][$CTD_svel]);

	push(@{$CTD{W}},      $ants_[0][$CTD_w]);

	push(@{$CTD{W_T}},	  $ants_[0][$CTD_w_t]);

	if ($ants_[0][$CTD_depth]+$opt_a > $CTD_maxdepth) {

		$CTD_maxdepth = $ants_[0][$CTD_depth]+$opt_a;

		$CTD_atbottom = $#{$CTD{DEPTH}};

	}

} while (&antsIn());

$CTD{DT} = ($CTD{ELAPSED}[$#{$CTD{ELAPSED}}] - $CTD{ELAPSED}[0]) / $#{$CTD{ELAPSED}};

progress("\t%d scans at %.1fHz\n",scalar(@{$CTD{DEPTH}}),1/$CTD{DT});

progress("\tstart depth = %.1fm\n",$CTD{DEPTH}[0]);

progress("\tmax depth   = %.1fm (# $CTD_atbottom)\n",$CTD_maxdepth);

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

# Construct sound-speed correction profile from CTD 1Hz downcast data

#	very simple algorithm that stores the last value found

#	in each 1m bin

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

progress("Constructing sound-speed correction profile\n");

my($scans_per_sec) = int(1/$CTD{DT}+0.5);

for (my($s)=0; $s<=$CTD_atbottom; $s+=$scans_per_sec) {

	next unless ($CTD{DEPTH}[$s] >= 0);

	$sVelProf[int($CTD{DEPTH}[$s])] = $CTD{SVEL}[$s];

}

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

# Determine time lag

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

if (defined($opt_i)) {

	progress("Setting initial time lag...\n");

	$CTD{TIME_LAG} = $opt_i;

	progress("\t-i => elapsed(CTD) ~ elapsed(LADCP) + %.1fs\n",$CTD{TIME_LAG});

} else {

	progress("Guestimating time lag...\n");

	my($CTD_10pct_down) = my($LADCP_10pct_down) = 0;

	$CTD_10pct_down++

		until ($CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0] >= 0.1*($CTD_maxdepth-$CTD{DEPTH}[0]));

	$LADCP_10pct_down++

		until ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} >= 0.1*$LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH});

	$CTD{TIME_LAG} = $CTD{ELAPSED}[$CTD_10pct_down] - $LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED};

	progress("\telapsed(dz(CTD)=%.1fm) ~ elapsed(dz(LADCP)=%.1fm) + %.1fs\n",

	    $CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH},$CTD{TIME_LAG});

}

$CTD{TIME_LAG} = calc_lag($n_lags[0],$w_size[0],int(1/$CTD{DT}+0.5));

progress("\telapsed(CTD) ~ elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG});

$CTD{TIME_LAG} = calc_lag($n_lags[1],$w_size[1],1);

progress("\telapsed(CTD) = elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG});

&antsAddParams('CTD_time_lag',$CTD{TIME_LAG});

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

# Merge CTD with LADCP data

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

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

progress("Merging CTD with LADCP data...\n");

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

	my($scan) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);

	if ($scan <= 0) {	# NB: must be <=, rather than <, or assertion below sometimes fails

		info("CTD data begin with LADCP in water\n") unless ($warned++);

		$firstGoodEns = $ens+1;

		next;

	}

	unless ($scan <= $#{$CTD{ELAPSED}}) {

		info("CTD data end with LADCP in water\n");

		$lastGoodEns = $ens-1;

		last;

	}

	die("assertion failed!\n" .

		"\ttest: abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2\n" .

		"\tens = $ens, scan = $scan\n" .

		sprintf("\tadjusted LADCP time = %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG}) .

		sprintf("\tCTD($scan) time = %f\n",$CTD{ELAPSED}[$scan]) .

		"=> Did you use SeaBird elapsed time? Don't!"

	) unless (abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2);

	$LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED} = $CTD{ELAPSED}[$scan];		# elapsed field for output

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

	    $LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $CTD{SVEL}[$scan]/1500; 	# correct (inadequately) for sound-speed variations

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

			$CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan]))

				unless ($CTD{DEPTH}[$scan] >= 0);

		$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan];

		$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan;

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

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

			$sumWsq += &SQR($reflr_ocean_w);

			$nWsq++;

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

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

					$sumWsqI += &SQR($reflr_ocean_w);

					$nWsqI++;

	        }

	    } else {

	    	undef($LADCP{ENSEMBLE}[$ens]->{CTD_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]->{CTD_DEPTH});							# old DEPTH from calcLADCPts()

	}

}

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

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

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

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

	warning(0,"%.2f cm/s reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI))

		if (sqrt($sumWsqI/$nWsqI) > 0.05);

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

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

} elsif ($nWsq > 0) {

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

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

} else {

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

}

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

# Find or set seabed

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

if (numberp($opt_h)) {

	progress("Setting water depth (-h)\n");

	$water_depth = $opt_h;

	$sig_water_depth = 0;

} else {

	progress("Finding seabed...\n");

	($water_depth,$sig_water_depth) =

	    find_seabed(\%LADCP,$LADCP_atbottom,$LADCP{BEAM_COORDINATES});

}

&antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth);

progress("\t%d(%d) m water depth\n",$water_depth,$sig_water_depth)

	if defined($water_depth);

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

# Remove data contaminated by sidelobe reflection from seabed and sea surface

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

if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {

	&antsAddParams('ADCP_orientation','downlooker');