#!/usr/bin/perl

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

#                    L A D C P I N T S H 

#                    doc: Thu Oct 14 21:22:50 2010

#                    (c) 2010 A.M. Thurnherr & E. Firing

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

$antsSummary = 'integrate LADCP shear';

# NOTES:

#	- the core of this code is a simplified version of avg_sh.m and

#	  int_sh.m written by Eric Firing

#	- comments beginning with ## are taken from Eric's code

#	- cubic velocity interpolation across PPI gap from Eric's code has

#	  not been implemented (yet?)

#	- low-pass-filtered shear code has not yet been implemented

#	- currently, shear gaps are assumed to have vanishing shear;

#	  better solutions are possible

#	- elapsed time is simply copied from shear elapsed time (i.e.

#	  it is not interpolated onto the new depth grid)

# WEIRDNESSES IN Eric's CODE:

#	- in Eric's [avg_sh.m] the calculation of output shear stddev incorrectly assumes that the 4th column

#	  in the shear profile is stddev, rather than variance. However, as far as I can tell, this output

#	  is not used anywhere in Eric's code

# HISTORY:

#	Oct 14, 2010: - created

#	Oct 20, 2010: - first working version

#	Oct 23, 2010: - added support for -b)

#	Oct 24, 2010: - fix spuriously small variances that can occur for BT velocities based on very small

#					samples (i.e. primarily when chosing a small -r)

#	Dec  9, 2010: - allowed for empty BT file

#	Feb 16, 2011: - BUG: gaps in shear data were not handled correctly in baroclinic solution

#	Jul  7, 2011: - added -m

#	Jul 19, 2011: - added shear sigma output for shear inversion

#				  - BUG: dc/uc v component was not correctly referenced

#	Jul 24, 2011: - BUG: BT constraint was erroneously assumed to be available for dc,uc only, if it

#						 was available for <dc,uc> combo profile

#	Jul 25, 2011: - BUG: nan in either of ul/dl u_z.sig caused combined u_z.sig to be nan, too

#	Jul 27, 2011: - removed code related to smoothed shear

#				  - replaced -r by -n

#				  - removed shear sigma output

#				  - replaced ndata by nsamp

#				  - removed -w (Eric's way of dealing with dc/uc temporal variability)

#	Feb 19, 2012: - added processing of elapsed time

#				  - adapted to new shear file layout (nsamp instead of nshear)

#				  - BUG: uplooker data was not used for downcasts and and only partially for combo data

#   May 16, 2012: - adapted to ANTSlib V5.0

#	May 25, 2012: - added code to read LDEO_IX bottom-track data

#	Jun 14, 2012: - noticed that -b now works also with LDEO SADCP files :-); renamed -b option to -r

#	Jun  5, 2013: - BUG: code bombed when either UC or DC was missing

#	Jun 28, 2013: - adapated to new :: convention

#				  - make sure LADCP DUL metadata are dealt with correctly

#	Jul 12, 2013: - clarified -u usage with better messages

#	Mar 20, 2014: - fiddled while debugging [LADCPproc]

#	Jun  7, 2014: - improved error messages

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

require "$ANTS/ants.pl";

require "$ANTS/libstats.pl";

&antsUsage('dm:n:r:s:w:u:',0,

	'[-d)ebug]',

	'[-r)eference with <BT or SADCP file> [-m)in <samp[10]>]]',

	'[secondary -u)plooker <shear file>]',

	'[min -n) <shear samp[10]>]',

	'[output -s)hear-pro <file>]',

	'[LADCP shear file]');

croak("$0: -m meaningless without -r\n")

	if defined($opt_m) && !defined($opt_r);

&antsCardOpt(\$opt_n,10);				# minimum number of samples for shear

$minBTsamp = &antsCardOpt($opt_m,10);	# minimum number of samples for BT data

&antsFileOpt($opt_r);		# reference velocity file

&antsFileOpt($opt_u);		# UL shear file

if (defined($opt_u)) {

	open(ULF,$opt_u) || croak("$opt_u: $!\n");

	%UL_P = &antsFileParams(ULF);

}

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

# Step 1: Read and Average Shear Data

#	- depth bins with less than $opt_n values are blanked out

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

sub wavg_sig(@)

{

	my($sumSq) = my($n) = 0;

	for (my($i)=0; $i<$#_; $i+=2) {

		next unless numberp($_[$i+1]);

		$sumSq += $_[$i] * $_[$i+1]**2;

		$n += $_[$i];

	}

	return ($n>0) ? sqrt($sumSq/$n) : nan;

}

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

# Handle Metadata

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

if (%UL_P) {

	croak("$0: inconsistent vertical resolution\n")

		unless ($P{LADCPproc::vertical_resolution} == $UL_P{LADCPproc::vertical_resolution});

	unless ($P{LADCPproc::bin_length} == $UL_P{LADCPproc::bin_length}) {

		&antsInfo("Warnining: different DL/UL bin lengths; derived spectra cannot be corrected");

		&antsAddParams('LADCPproc::bin_length','',

					   'LADCPproc::DL_bin_length',$P{LADCPproc::bin_length},

					   'LADCPproc::UL_bin_length',$UL_P{LADCPproc::bin_length});

	}

}

$depthF = fnrNoErr('depth');									# layout of [LADCPproc] output

unless (defined($depthF)) {

	if (defined($opt_u)) {

		croak("No 'depth' field in primary shear file (extraneous -u?)\n");

	} else {

		croak("No 'depth' field in primary shear file\n");

	}

}

$dc_nshF = fnrNoErr('dc_nshear');

$dc_nshF = fnr('dc_nsamp') unless defined($dc_nshF);

$dc_uzF  = fnr('dc_u_z');

$dc_uzsF = fnrNoErr('dc_u_z.sig');

$dc_uzsF = fnr('dc_u_z_sig') unless defined($dc_uzsF);

$dc_vzF  = fnr('dc_v_z');

$dc_vzsF = fnrNoErr('dc_v_z.sig');

$dc_vzsF = fnr('dc_v_z_sig') unless defined($dc_vzsF);

$dc_wzF  = fnr('dc_w_z');

$dc_wzsF = fnrNoErr('dc_w_z.sig');

$dc_wzsF = fnr('dc_w_z_sig') unless defined($dc_wzsF);

$dc_elapsedF = fnr('dc_elapsed');

$uc_nshF = fnrNoErr('uc_nshear');

$uc_nshF = fnr('uc_nsamp') unless defined($uc_nshF);

$uc_uzF  = fnr('uc_u_z');

$uc_uzsF = fnrNoErr('uc_u_z.sig');

$uc_uzsF = fnr('uc_u_z_sig') unless defined($uc_uzsF);

$uc_vzF  = fnr('uc_v_z');

$uc_vzsF = fnrNoErr('uc_v_z.sig');

$uc_vzsF = fnr('uc_v_z_sig') unless defined($uc_vzsF);

$uc_wzF  = fnr('uc_w_z');

$uc_wzsF = fnrNoErr('uc_w_z.sig');

$uc_wzsF = fnr('uc_w_z_sig') unless defined($uc_wzsF);

$uc_elapsedF = fnr('uc_elapsed');

my(@gaps); my($curGap) = 0;

for (my($r)=0; &antsIn(); $r++) {

	my(@UL_);

	if (defined($opt_u)) {

		@UL_ = &antsFileIn(ULF);							# read UL shear data

		undef($opt_u) unless (@UL_);						# cheap trick

	}

	$depth[$r] = $ants_[0][$depthF];						## depth grid values

	croak("$opt_u: inconsistent depth record $r (DL: $depth[$r]; UL: $UL_[$depthF])\n")

		if defined($opt_u) && ($UL_[$depthF] != $depth[$r]);

	$dc_nsamp = $ants_[0][$dc_nshF];						# number of shear samples

	$uc_nsamp = $ants_[0][$uc_nshF];

	if (defined($opt_u)) {

		$dl_nsamp = $dc_nsamp + $uc_nsamp;

		$dc_nsamp += $UL_[$dc_nshF];

		$uc_nsamp += $UL_[$uc_nshF];

		$ul_nsamp = $dc_nsamp + $uc_nsamp;

	}

	$dc_nsamp[$r] = $dc_nsamp;								# save for each record

	$uc_nsamp[$r] = $uc_nsamp;

	$nsamp[$r] = $dc_nsamp + $uc_nsamp;

	if (defined($opt_u)) {

		$ul_nsamp[$r] = $ul_nsamp;

		$dl_nsamp[$r] = $dl_nsamp;

	}

	if (defined($opt_u)) {									# dual-head instrument

		if ($dc_nsamp > 0) {								# downcast shear

			my($DLf) = $ants_[0][$dc_nshF] / $dc_nsamp;

			my($ULf) =      $UL_[$dc_nshF] / $dc_nsamp;

			if ($DLf>0 && $ULf>0) {

				$dc_uz[$r] = $DLf*$ants_[0][$dc_uzF] + $ULf*$UL_[$dc_uzF];

				$dc_vz[$r] = $DLf*$ants_[0][$dc_vzF] + $ULf*$UL_[$dc_vzF];

	            $dc_wz[$r] = $DLf*$ants_[0][$dc_wzF] + $ULf*$UL_[$dc_wzF];

	            $dc_elapsed[$r] = $DLf*$ants_[0][$dc_elapsedF] + $ULf*$UL_[$dc_elapsedF];

	        } elsif ($DLf > 0) {

				$dc_uz[$r] = $ants_[0][$dc_uzF];

				$dc_vz[$r] = $ants_[0][$dc_vzF];

				$dc_wz[$r] = $ants_[0][$dc_wzF];

				$dc_elapsed[$r] = $ants_[0][$dc_elapsedF];

	        } else {

				$dc_uz[$r] = $UL_[$dc_uzF];

				$dc_vz[$r] = $UL_[$dc_vzF];

				$dc_wz[$r] = $UL_[$dc_wzF];

				$dc_elapsed[$r] = $UL_[$dc_elapsedF];

	        }

		} else {

			$dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = $dc_elapsed[$r] = nan;

	    }

		if ($uc_nsamp > 0) {								# upcast shear

			my($DLf) = $ants_[0][$uc_nshF] / $uc_nsamp;

			my($ULf) =      $UL_[$uc_nshF] / $uc_nsamp;

			if ($DLf>0 && $Ulf>0) {

				$uc_uz[$r] = $DLf*$ants_[0][$uc_uzF] + $ULf*$UL_[$uc_uzF];

				$uc_vz[$r] = $DLf*$ants_[0][$uc_vzF] + $ULf*$UL_[$uc_vzF];

				$uc_wz[$r] = $DLf*$ants_[0][$uc_wzF] + $ULf*$UL_[$uc_wzF];

				$uc_elapsed[$r] = $DLf*$ants_[0][$uc_elapsedF] + $ULf*$UL_[$uc_elapsedF];

	        } elsif ($DLf > 0) {

				$uc_uz[$r] = $ants_[0][$uc_uzF];

				$uc_vz[$r] = $ants_[0][$uc_vzF];

				$uc_wz[$r] = $ants_[0][$uc_wzF];

				$uc_elapsed[$r] = $ants_[0][$uc_elapsedF];

	        } else {

				$uc_uz[$r] = $UL_[$uc_uzF];

				$uc_vz[$r] = $UL_[$uc_vzF];

				$uc_wz[$r] = $UL_[$uc_wzF];

				$uc_elapsed[$r] = $UL_[$uc_elapsedF];

	        }

		} else {

			$uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = $uc_elapsed[$r] = nan;

	    }

	} else {	# downlooker only

		if ($dc_nsamp > 0) {							# downcast shear

			$dc_uz[$r] = $ants_[0][$dc_uzF];

			$dc_vz[$r] = $ants_[0][$dc_vzF];

			$dc_wz[$r] = $ants_[0][$dc_wzF];

			$dc_elapsed[$r] = $ants_[0][$dc_elapsedF];

		} else {

			$dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = $dc_elapsed[$r] = nan;

	    }

		if ($uc_nsamp > 0) {							# upcast shear

			$uc_uz[$r] = $ants_[0][$uc_uzF];

			$uc_vz[$r] = $ants_[0][$uc_vzF];

			$uc_wz[$r] = $ants_[0][$uc_wzF];

			$uc_elapsed[$r] = $ants_[0][$uc_elapsedF];

		} else {

			$uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = $uc_elapsed[$r] = nan;

	    }

    }

	if ($nsamp[$r] > 0) {

		my($dcf) = $dc_nsamp / $nsamp[$r];

		my($ucf) = $uc_nsamp / $nsamp[$r];

		if ($dcf>0 && $ucf>0) {

			$uz[$r] = $dcf*$dc_uz[$r] + $ucf*$uc_uz[$r];

			$vz[$r] = $dcf*$dc_vz[$r] + $ucf*$uc_vz[$r];

            $wz[$r] = $dcf*$dc_wz[$r] + $ucf*$uc_wz[$r];

            $elapsed[$r] = $dcf*$dc_elapsed[$r] + $ucf*$uc_elapsed[$r];

        } elsif ($dcf > 0) {

			$uz[$r] = $dc_uz[$r];

			$vz[$r] = $dc_vz[$r];

			$wz[$r] = $dc_wz[$r];

			$elapsed[$r] = $dc_elapsed[$r];

        } else {

			$uz[$r] = $uc_uz[$r];

			$vz[$r] = $uc_vz[$r];

			$wz[$r] = $uc_wz[$r];

			$elapsed[$r] = $uc_elapsed[$r];

        }

	} else {

		$uz[$r] = $vz[$r] = $wz[$r] = $elapsed[$r] = nan;

	}

	if (numberp($uz[$r]) && $curGap>0) {						# end of gap

		push(@gaps,$curGap)	unless ($r == $curGap);				# do not report "gap" at beginning of profile

#		print(STDERR "$curGap-gap at $depth[$r]m\n");

		$curGap = 0;

    } elsif (!numberp($uz[$r])) {								# currently in gap

#    	print(STDERR "in gap at $depth[$r]m (nsamp = $nsamp[$r], $dc_nsamp,$uc_nsamp)\n");

		$curGap++;

    }

	if ($nsamp[$r] > 0) {

		if (defined($opt_u)) {

			$uzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_uzsF],

							           $UL_[$dc_nshF],     $UL_[$dc_uzsF],

							      $ants_[0][$uc_nshF],$ants_[0][$uc_uzsF],

							           $UL_[$uc_nshF],     $UL_[$uc_uzsF]);

			$vzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_vzsF],

									   $UL_[$dc_nshF],     $UL_[$dc_vzsF],

							      $ants_[0][$uc_nshF],$ants_[0][$uc_vzsF],

									   $UL_[$uc_nshF],     $UL_[$uc_vzsF]);

			$wzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_wzsF],

									   $UL_[$dc_nshF],     $UL_[$dc_wzsF],

							      $ants_[0][$uc_nshF],$ants_[0][$uc_wzsF],

	                                   $UL_[$uc_nshF],     $UL_[$uc_wzsF]);

		} else { # DL only

			$uzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_uzsF],

							      $ants_[0][$uc_nshF],$ants_[0][$uc_uzsF]);

			$vzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_vzsF],

							      $ants_[0][$uc_nshF],$ants_[0][$uc_vzsF]);

			$wzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_wzsF],

	                              $ants_[0][$uc_nshF],$ants_[0][$uc_wzsF]);

	    }

	} else {

		$uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan;

	}

	if ($dc_nsamp > 0) {									# same calc for downcast only

		if (defined($opt_u)) {

			$dc_uzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_uzsF],

							              $UL_[$dc_nshF],     $UL_[$dc_uzsF]);

			$dc_vzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_vzsF],

							              $UL_[$dc_nshF],     $UL_[$dc_vzsF]);

			$dc_wzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_wzsF],

							              $UL_[$dc_nshF],     $UL_[$dc_wzsF]);

		} else {

			$dc_uzsig[$r] = $ants_[0][$dc_uzsF];

			$dc_vzsig[$r] = $ants_[0][$dc_vzsF];

			$dc_wzsig[$r] = $ants_[0][$dc_wzsF];

	    }

	} else {

		$uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan;

	}

	if ($uc_nsamp > 0) {									# same calc for upcast only

		if (defined($opt_u)) {

			$uc_uzsig[$r] = wavg_sig($ants_[0][$uc_nshF],$ants_[0][$uc_uzsF],

							              $UL_[$uc_nshF],     $UL_[$uc_uzsF]);

			$uc_vzsig[$r] = wavg_sig($ants_[0][$uc_nshF],$ants_[0][$uc_vzsF],

							              $UL_[$uc_nshF],     $UL_[$uc_vzsF]);

			$uc_wzsig[$r] = wavg_sig($ants_[0][$uc_nshF],$ants_[0][$uc_wzsF],

							              $UL_[$uc_nshF],     $UL_[$uc_wzsF]);

		} else {

			$uc_uzsig[$r] = $ants_[0][$uc_uzsF];

			$uc_vzsig[$r] = $ants_[0][$uc_vzsF];

			$uc_wzsig[$r] = $ants_[0][$uc_wzsF];

	    }

	} else {

		$uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan;

	}

}

if (@gaps) {

	&antsAddParams('shear_gaps',"@gaps");

	print(STDERR "shear gaps: @gaps\n");

} else {

	&antsAddParams('shear_gaps',0);

}

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

# Step 2: Low-Pass filter high-quality shear data; not yet implemented

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

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

# Step 3: Integrate Shear

#	- z(vel) = z(sh) + DZ/2

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

my($DZ) = $depth[1] - $depth[0];

for (my($r)=my($u)=my($v)=my($w)=my($dc_u)=my($dc_v)=my($dc_w)=my($uc_u)=my($uc_v)=my($uc_w)=0;

	 $r<@depth; $r++) {

	if ($nsamp[$r] >= $opt_n) {

		$u = $u[$r] = $u + $DZ*$uz[$r];

		$v = $v[$r] = $v + $DZ*$vz[$r];

		$w = $w[$r] = $w + $DZ*$wz[$r];

	}

	if ($dc_nsamp[$r] >= $opt_n) {

		$dc_u = $dc_u[$r] = $dc_u + $DZ*$dc_uz[$r];

		$dc_v = $dc_v[$r] = $dc_v + $DZ*$dc_vz[$r];

		$dc_w = $dc_w[$r] = $dc_w + $DZ*$dc_wz[$r];

	}

	if ($uc_nsamp[$r] >= $opt_n) {

		$uc_u = $uc_u[$r] = $uc_u + $DZ*$uc_uz[$r];

		$uc_v = $uc_v[$r] = $uc_v + $DZ*$uc_vz[$r];

		$uc_w = $uc_w[$r] = $uc_w + $DZ*$uc_wz[$r];

	}

}

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

# Step 4: Reference Velocities

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

my($refU,$refV,$refW,$dc_refU,$dc_refV,$dc_refW,$uc_refU,$uc_refV,$uc_refW);

if (defined($opt_r)) {											# reference using velocity profile

	print(STDERR "Loading reference-velocity data from $opt_r...\n")

		if ($opt_d);

	open(BTF,$opt_r) || croak("$opt_r: $!\n");

	my(@BTL) = &antsFileLayout(BTF);

	if (@BTL) {													# valid ANTS file

		my($BTdF,$BTndF,$BTuF,$BTvF,$BTwF,$BTusF,$BTvsF,$BTwsF);

		for (my($f)=0; $f<@BTL; $f++) {

			$BTdF = $f if ($BTL[$f] eq 'depth');

			$BTuF = $f if ($BTL[$f] eq 'u');

			$BTvF = $f if ($BTL[$f] eq 'v');

			$BTwF = $f if ($BTL[$f] eq 'w');

			$BTusF = $f if ($BTL[$f] eq 'u.sig');

			$BTvsF = $f if ($BTL[$f] eq 'v.sig');

			$BTwsF = $f if ($BTL[$f] eq 'w.sig');

			$BTndF = $f if ($BTL[$f] eq 'nsamp');

			$BTerrF= $f if ($BTL[$f] eq 'err');

		}

		if (defined($BTdF) && defined($BTuF) &&					# from LADCPproc

		    defined($BTvF) && defined($BTwF) &&

		    defined($BTusF) && defined($BTvsF) &&

			defined($BTwsF) && defined($BTndF)) {	

			while (my(@BTr) = &antsFileIn(BTF)) {

				my($gi) = int($BTr[$BTdF] / $DZ);

				next unless ($BTr[$BTndF] >= $minBTsamp);

				$BT_nsamp[$gi] = $BTr[$BTndF];

				$BT_u[$gi] = $BTr[$BTuF];

				$BT_v[$gi] = $BTr[$BTvF];

				$BT_w[$gi] = $BTr[$BTwF];

				$BT_u_var[$gi] = $BTr[$BTusF]**2;

				$BT_v_var[$gi] = $BTr[$BTvsF]**2;

				$BT_w_var[$gi] = $BTr[$BTwsF]**2;

			}

			&fixLowSampStat(\@BT_u_var,@BT_nsamp);				# remove spurious small variances

			&fixLowSampStat(\@BT_v_var,@BT_nsamp);

	        &fixLowSampStat(\@BT_w_var,@BT_nsamp);

        } elsif (defined($BTdF) && defined($BTuF) &&			# LDEO_IX ANTS format

				 defined($BTvF) && defined($BTerrF)) {

			croak("$0: -m not supported for LDEO_IX output\n")

				if defined($opt_m);

			while (my(@BTr) = &antsFileIn(BTF)) {

				my($gi) = int($BTr[$BTdF] / $DZ);

				$BT_u[$gi] = $BTr[$BTuF];

				$BT_v[$gi] = $BTr[$BTvF];

				$BT_u_var[$gi] = $BTr[$BTerrF]**2;

				$BT_v_var[$gi] = $BTr[$BTerrF]**2;

#				print(STDERR "$gi $BT_u[$gi] $BT_v[$gi] $BT_u_var[$gi] $BT_v_var[$gi]\n");

			}

        } else {

			croak("$opt_r: not a valid reference-velocity file (ANTS format)\n");

		}

	} else {													# non-ANTS file (LDEO_IX assumed)

		croak("$0: -m not supported for LDEO_IX output\n")

			if defined($opt_m);

		while (<BTF>) {

			last if /^Columns\s+=\s+z:u:v:err/;

		}

		croak("$opt_r: not a valid reference-velocity file (non-ANTS format)\n")

			unless /^Columns\s+=\s+z:u:v:err/;

		while (<BTF>) {

			my($depth,$u,$v,$err) = split;

			croak("$0: cannot handle non-numeric BT data (implementation restriction)\n")

				unless numberp($depth) && numberp($u) && numberp($v) && numberp($err);

			my($gi) = int($depth / $DZ);

			$BT_u[$gi] = $u;

			$BT_v[$gi] = $v;

			$BT_u_var[$gi] = $err**2;

			$BT_v_var[$gi] = $err**2;

#			print(STDERR "$gi $BT_u[$gi] $BT_v[$gi] $BT_u_var[$gi] $BT_v_var[$gi]\n");

		}

    }

    close(BTF);

	my($sumU,$sumV,$sumW,$dc_sumU,$dc_sumV,$dc_sumW,		# average integrated-shear velocities

	   $uc_sumU,$uc_sumV,$uc_sumW);

	my($nSumVel,$dc_nSumVel,$uc_nSumVel);

	my($wSumBTu,$wSumBTv,$wSumBTw);							# weighted sums of BT-ref'd velocities

	my($dc_wSumBTu,$dc_wSumBTv,$dc_wSumBTw);

	my($uc_wSumBTu,$uc_wSumBTv,$uc_wSumBTw);

	my($sumVarBTu,$sumVarBTv,$sumVarBTw);					# sum of variances of BT-ref'd vels

	my($dc_sumVarBTu,$dc_sumVarBTv,$dc_sumVarBTw);

	my($uc_sumVarBTu,$uc_sumVarBTv,$uc_sumVarBTw);

	for (my($r)=0; $r<@depth; $r++) {

		if (numberp($BT_u[$r]) && numberp($u[$r])) {

			$nSumVel++;

			$sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r];

			$wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $sumVarBTu += 1/$BT_u_var[$r];

			$wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $sumVarBTv += 1/$BT_v_var[$r];

			if (@BT_w) {

				$wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $sumVarBTw += 1/$BT_w_var[$r];

			}

		}

		if (numberp($BT_u[$r]) && numberp($dc_u[$r])) {

			$dc_nSumVel++;

			$dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r];

			$dc_wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $dc_sumVarBTu += 1/$BT_u_var[$r];

			$dc_wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $dc_sumVarBTv += 1/$BT_v_var[$r];

			if (@BT_w) {

				$dc_wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $dc_sumVarBTw += 1/$BT_w_var[$r];

			}

		}

		if (numberp($BT_u[$r]) && numberp($uc_u[$r])) {

			$uc_nSumVel++;

			$uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r];

			$uc_wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $uc_sumVarBTu += 1/$BT_u_var[$r];

			$uc_wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $uc_sumVarBTv += 1/$BT_v_var[$r];

			if (@BT_w) {