2
|
1 |
#======================================================================
|
|
2 |
# B O T T O M _ T R A C K I N G . P L
|
|
3 |
# doc: Wed Oct 20 21:05:37 2010
|
|
4 |
# dlm: Fri Dec 31 14:21:14 2010
|
|
5 |
# (c) 2010 A.M. Thurnherr
|
|
6 |
# uE-Info: 32 30 NIL 0 0 72 2 2 4 NIL ofnI
|
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# HISTORY:
|
|
10 |
# Oct 20, 2010: - created
|
|
11 |
# Dec 30, 2010: - adapted for use with LADCP_w
|
|
12 |
|
|
13 |
# This code is essentially identical to the one used in LADCPproc. Differences:
|
|
14 |
# 1) velocity editing is simpler: no wake editing, no PPI editing, no shear
|
|
15 |
# editing, no w outlier
|
|
16 |
# 2) median/mad calculated instead of mean/stddev
|
|
17 |
# 3) u,v not calculated
|
|
18 |
|
|
19 |
# Don't look for BT-referenced velocities if package is more than $BT_max_range
|
|
20 |
# above seabed. This parameter is frequency dependent and the current value is
|
|
21 |
# appropriate (if rather high) for 300kHz Workhorse intruments.
|
|
22 |
|
|
23 |
my($BT_max_range) = 300;
|
|
24 |
|
|
25 |
|
|
26 |
# The code only tries to bin BT-referenced velocities if a consistent bottom
|
|
27 |
# is available in all 4 beams. Ensembles where the range of bin numbers where
|
|
28 |
# the maximum echo is found is greater than $max_BIT_bin_range_diff are rejected.
|
|
29 |
# In addition to flukes this also rejects ensembles collected with large
|
|
30 |
# instrument tilts. The value of 3 is a first guess that has not been explored.
|
|
31 |
|
|
32 |
my($BT_max_bin_range_diff) = 3;
|
|
33 |
|
|
34 |
|
|
35 |
# If the difference between measured vertical velocity of the seabed (i.e.
|
|
36 |
# the package vertical velocity referenced by the seabed) and the vertical
|
|
37 |
# velocity of the CTD (from dp/dt) si greater than $BT_max_w_error the current
|
|
38 |
# ensemble is ignored and $nBTwFlag is increased. The value of
|
|
39 |
# 3cm/s is taken from listBT developed on A0304 cruise.
|
|
40 |
|
|
41 |
my($BT_max_w_error) = 0.03;
|
|
42 |
|
|
43 |
#----------------------------------------------------------------------
|
|
44 |
# bin valid BT-referenced velocities
|
|
45 |
# input: ensemble number, water depth (with uncertainty)
|
|
46 |
# output: @BTu,@BTv,@BTw main result
|
|
47 |
# editflags for information
|
|
48 |
# @BTbtmu, @BTbtmv, @BTtilt stats to find reasons for quality differences
|
|
49 |
#----------------------------------------------------------------------
|
|
50 |
|
|
51 |
my($nBTfound,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0);
|
|
52 |
my(@BTu,@BTv,@BTw);
|
|
53 |
|
|
54 |
sub binBTprof($$$)
|
|
55 |
{
|
|
56 |
my($ens,$wd,$sig_wd) = @_;
|
|
57 |
|
|
58 |
my(@ea_max) = (0,0,0,0); my(@ea_max_bin) = (nan,nan,nan,nan);
|
|
59 |
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
|
|
60 |
$ea_max[0] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0],
|
|
61 |
$ea_max_bin[0] = $bin
|
|
62 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0] > $ea_max[0]);
|
|
63 |
$ea_max[1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1],
|
|
64 |
$ea_max_bin[1] = $bin
|
|
65 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1] > $ea_max[1]);
|
|
66 |
$ea_max[2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2],
|
|
67 |
$ea_max_bin[2] = $bin
|
|
68 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2] > $ea_max[2]);
|
|
69 |
$ea_max[3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3],
|
|
70 |
$ea_max_bin[3] = $bin
|
|
71 |
if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3] > $ea_max[3]);
|
|
72 |
}
|
|
73 |
|
|
74 |
return # disregard boundary maxima
|
|
75 |
unless (min(@ea_max_bin) > $LADCP_firstBin-1) &&
|
|
76 |
(max(@ea_max_bin) < $LADCP_lastBin-1);
|
|
77 |
return # inconsistent range to seabed
|
|
78 |
unless (max(@ea_max_bin)-min(@ea_max_bin) <= $BT_max_bin_range_diff);
|
|
79 |
|
|
80 |
$nBTfound++;
|
|
81 |
my($seafloor_bin) = round(avg(@ea_max_bin));
|
|
82 |
|
|
83 |
my(@bd) = calc_binDepths($ens);
|
|
84 |
$nBTdepthFlag++,return # BT range inconsistent with water depth
|
|
85 |
unless (abs($wd-$bd[$seafloor_bin]) < max($sig_wd,$LADCP{BIN_LENGTH}));
|
|
86 |
|
|
87 |
# try vertical velocities at seabed bin plus one above and below
|
|
88 |
# this does not really work because, often, only one of the bins has valid velocities
|
|
89 |
my($w1) = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin-1];
|
|
90 |
my($w2) = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin ];
|
|
91 |
my($w3) = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin+1];
|
|
92 |
|
|
93 |
$w1 = 9e99 unless numberp($w1); # invalid velocity sentinels
|
|
94 |
$w2 = 9e99 unless numberp($w1);
|
|
95 |
$w3 = 9e99 unless numberp($w1);
|
|
96 |
|
|
97 |
my($seafloor_u,$seafloor_v,$seafloor_w);
|
|
98 |
|
|
99 |
# determine which of the three trial bins is most consistent with reflr vertical velocities
|
|
100 |
die("assertion failed") unless numberp($LADCP{ENSEMBLE}[$ens]->{REFLR_W});
|
|
101 |
if (abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w2) &&
|
|
102 |
abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w3)) {
|
|
103 |
$seafloor_u = $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$seafloor_bin-1];
|
|
104 |
$seafloor_v = $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$seafloor_bin-1];
|
|
105 |
$seafloor_w = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin-1];
|
|
106 |
} elsif (abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w2)) {
|
|
107 |
$seafloor_u = $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$seafloor_bin+1];
|
|
108 |
$seafloor_v = $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$seafloor_bin+1];
|
|
109 |
$seafloor_w = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin+1];
|
|
110 |
} else {
|
|
111 |
$nBTvalidVelFlag++,return # none of 3 trial bins has valid velocity
|
|
112 |
if ($w2 == 9e99);
|
|
113 |
$seafloor_u = $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$seafloor_bin];
|
|
114 |
$seafloor_v = $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$seafloor_bin];
|
|
115 |
$seafloor_w = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin];
|
|
116 |
}
|
|
117 |
|
|
118 |
$nBTwFlag++,return # velocity error is too great
|
|
119 |
if (abs($seafloor_w-$LADCP{ENSEMBLE}[$ens]->{REFLR_W}) > $BT_max_w_error);
|
|
120 |
|
|
121 |
push(@BTbtmu,$seafloor_u); # calc stats to try to find reasons for quality
|
|
122 |
push(@BTbtmv,$seafloor_v);
|
|
123 |
push(@BTbtmw,$seafloor_w);
|
|
124 |
push(@BTtilt,$LADCP{ENSEMBLE}[$ens]->{TILT});
|
|
125 |
|
|
126 |
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
|
|
127 |
next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
|
|
128 |
my($gi) = int($bd[$bin]) / $opt_o;
|
|
129 |
push(@{$BTw[$gi]},$LADCP{ENSEMBLE}[$ens]->{W}[$bin]-$seafloor_w);
|
|
130 |
}
|
|
131 |
}
|
|
132 |
|
|
133 |
#----------------------------------------------------------------------
|
|
134 |
# calculate BT-referenced velocity profile
|
|
135 |
# input: start,end LADCP ensembles, water depth with uncertainty
|
|
136 |
# output: %BT{MEDIAN_W,MAD_W,N_SAMP}
|
|
137 |
#----------------------------------------------------------------------
|
|
138 |
|
|
139 |
sub calc_BTprof($$$$)
|
|
140 |
{
|
|
141 |
my($LADCP_start,$LADCP_end,$wd,$sig_wd) = @_;
|
|
142 |
|
|
143 |
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
|
|
144 |
next unless ($wd-$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} < $BT_max_range);
|
|
145 |
binBTprof($ens,$wd,$sig_wd);
|
|
146 |
}
|
|
147 |
|
|
148 |
progress("\t$nBTfound BT ensembles found\n");
|
|
149 |
progress("\t$nBTdepthFlag flagged bad because of wrong bottom depth\n");
|
|
150 |
progress("\t$nBTvalidVelFlag flagged bad because of no valid velocities\n");
|
|
151 |
progress("\t$nBTwFlag flagged bad because of incorrect vertical velocity\n");
|
|
152 |
|
|
153 |
for (my($gi)=0; $gi<@BTw; $gi++) { # calc grid medians & mads
|
|
154 |
$BT{N_SAMP}[$gi] = @{$BTw[$gi]};
|
|
155 |
$BT{MEDIAN_W}[$gi] = median(@{$BTw[$gi]});
|
|
156 |
$BT{MAD_W}[$gi] = mad2($BT{W}[$gi],@{$BTw[$gi]});
|
|
157 |
}
|
|
158 |
|
|
159 |
&antsAddParams('BT_rms_seafloor_u',round(rms(@BTbtmu),0.01),
|
|
160 |
'BT_rms_seafloor_v',round(rms(@BTbtmv),0.01),
|
|
161 |
'BT_rms_seafloor_w',round(rms(@BTbtmw),0.01),
|
|
162 |
'BT_avg_seafloor_u',round(avg(@BTbtmu),0.01),
|
|
163 |
'BT_avg_seafloor_v',round(avg(@BTbtmv),0.01),
|
|
164 |
'BT_avg_seafloor_w',round(avg(@BTbtmw),0.01),
|
|
165 |
'BT_rms_tilt',round(rms(@BTtilt),0.1));
|
|
166 |
}
|
|
167 |
|
|
168 |
1;
|