--- a/LADCP_w Fri Jun 17 13:33:47 2011 -0400
+++ b/LADCP_w Sat Jul 09 15:14:33 2011 -0400
@@ -2,19 +2,22 @@
#======================================================================
# L A D C P _ W
# doc: Fri Dec 17 18:11:13 2010
-# dlm: Wed Feb 16 14:03:53 2011
+# dlm: Fri Jul 8 05:11:59 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 45 29 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 67 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
-# calculate CTD acceleration from CTD velocity
-# remove water-depth from BT code, which is not really used and a bit of an outlier
-# because mean and stddev are used instead of median/mad
+# - make timelagging work for short casts (make sure 10% is not more than half window size)
+# - own seabed detection (P403)
+# - use instrument tilt in sidelobe editing?
+# - remove profile-binning code but add option to output BT-referenced data
+# - make upcast-flag valid for yoyo casts
+# - make diagnostic output 3-beam field work for Earth coordinates
+# - read assumed ADCP soundspeed from data file, instead of assuming 1500m/s
+# - calculate CTD acceleration from CTD velocity; probably not useful
+# - remove water-depth from BT code, which is not really used and a bit of an outlier
+# because mean and stddev are used instead of median/mad
$antsSummary = 'calculate vertical velocities from LADCP & CTD time series';
@@ -43,6 +46,33 @@
# - BUG: division by zero if no valid data
# Jan 5, 2010: - adapted to allow for gaps in CTD time series
# Feb 16, 2011: - cosmetics
+# Jun 22, 2011: - cosmetics
+# Jun 23, 2011: - disabled error on large rms reflr w
+# - added -l
+# - removed CTD headers from output
+# Jun 26, 2011: - added -u
+# - changed package correction from acceleration to velocity, because of
+# Stan's Antarctic data set where accelerations are zero but package effects are
+# there
+# Jul 2, 2011: - increased tilt default to 15 degrees
+# Jul 3, 2011: - replaced old package-velocity correction -x code by new beamvel correction
+# - removed -p
+# - replaced -d by residual (diagnostics) output
+# Jul 4, 2011: - saved a snapshot in [Jul_04_2011]
+# - removed output of binned profile (but not calculation code, which is required for residual)
+# - BUG: firstGoodEns or lastGoodEns could end up in a reflr_w gap when valid LADCP data begin
+# or end with CTD in water
+# - moved rarely used option -s to -k
+# - added -s)kip <ens> option
+# - had to very very slightly relax an assertion (by 1e-10 seconds...)
+
+# PROCESSING PARAMETER FILE
+# - # is comment character
+# - invalid entries ignored
+# - valid entries begin with <ADCP-file>: (NO INITIAL SPACE ALLOWED)
+# - remainder of line is added to usage before ADCP file and LADCP_w is restarted
+# - no argument with spaces allowed!
+# - -0 suppresses acting on -l
# CTD REQUIREMENTS
# - elapsed elapsed seconds; see note below
@@ -73,6 +103,15 @@
# - as a result, a profile only starts with elapsed==0 if the CTD
# is turned on when the LADCP is already in the water
+# OUTPUT
+# - default is "long" output with all the data; required to generate Visbeck-type plots
+# - residuals are calculated with respect to down-/upcast medians in binned profiles
+# - -p requests old profile output
+# - BT-referenced profiles are only found in -p output
+# - equivalence, assuming -o 10:
+# 1) 004DL.prof dc_w depth
+# 2) bindata -Sdowncast:1 -Fw.median,depth -n 20 depth 0 1000 10 004DL.w
+
# TIME LAGGING
# - occasionally, the time lagging algorithm fails, in particular if the CTD is turned on
# some time after the package enters the water
@@ -110,8 +149,11 @@
require "$PERL_TOOLS/RDI_BB_Read.pl";
require "$PERL_TOOLS/RDI_Coords.pl";
+@ARGS = @ARGV; # save opts & arguments
+
$antsParseHeader = 0;
-&antsUsage('3:4a:b:c:d:e:f:g:h:i:m:n:o:p:r:s:t:v:w:x:',1,
+&antsUsage('03:4a:b:c:e:f:g:h:i:k:l:m:n:o:p:r:s:t:uv:w:x:',1,
+ '[-l)oad processing parameters from <file>]',
'[-v)erbosity <level[1]>]',
'[require -4)-beam solutions]',
'[-r)ef-layer <bin[2],bin[6]>]',
@@ -120,13 +162,14 @@
'[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>]',
+ '[-i)nitial CTD time offset <guestimate> [-u)se as final]]',
'[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_wave_corr_coeffs/uc_wave_corr_coeffs>]',
+# '[-x <pkg_vel_corr_coeffs>]',
+ '[-x <beamvel-scale-fac>]',
'[valid LADCP -b)ins <bin[2],bin[*]>',
- '[-o)utput bin <resolution[50m]>] [require -s)amples <min[20]>]',
- '[-f) write time-series <file>] [-d)ump depth-bins to <basename>] [-p)ackage effect <file>]',
+ '[-o)utput bin <resolution[50m]>] [-k) require <min[20]> samples]',
+ '[-f) write time-series <file>] [output binned -p)rofile to <file>]',
'<LADCP-file> [CTD-file]');
&antsCardOpt(\$opt_v,1); # suppress regular info
@@ -134,19 +177,24 @@
$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_t,15); # 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
+&antsFloatOpt($opt_i); # externally supplied lag
+&antsUsageError() if ($opt_u && !defined($opt_i));
+
+&antsCardOpt(\$opt_s,0); # skip # initial ensembles
+
$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
+&antsCardOpt(\$opt_k,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
@@ -168,25 +216,37 @@
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_x)) { # decode corrections
+# (@pkg_vel_corr_poly) = split(',',$opt_x);
+# croak("$0: cannot decode -x $opt_x\n")
+# unless (@pkg_vel_corr_poly);
+# &antsAddParams('pkg_vel_corr_intercept',$pkg_vel_corr_poly[0],'pkg_vel_corr_slope',$pkg_vel_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"));
-}
+&antsFloatOpt(\$opt_x,1);
$LADCP_file = &antsFileArg();
+if (defined($opt_l) && !defined($opt_0)) { # load per-cast processing parameters
+ my(@cast_params);
+ open(PF,$opt_l) || croak("$opt_l: $!\n");
+ while (<PF>) {
+ s/#.*//;
+ @cast_params=split(/\s+/),last if /^$LADCP_file:/;
+ }
+ close(PF);
+
+ if (@cast_params) { # found valid entry
+ if ($ARGS[$#ARGS] eq $LADCP_file) { # CTD data on <stdin>
+ exec($0,@ARGS[0..$#ARGS-1],@cast_params[1..$#cast_params],'-0',$ARGS[$#ARGS]);
+ die("exec($0,@ARGS[0..$#ARGS-1],@cast_params[1..$#cast_params],'-0',$ARGS[$#ARGS]) failed!\n");
+ } else { # CTD file specified on cmdline
+ exec($0,@ARGS[0..$#ARGS-2],@cast_params[1..$#cast_params],'-0',$LADCP_file,$ARGS[$#ARGS]);
+ die("exec($0,@ARGS[0..$#ARGS-2],@cast_params[1..$#cast_params],'-0',$LADCP_file,$ARGS[$#ARGS]) failed!\n");
+ }
+ }
+}
+
#----------------------------------------------------------------------
# Screen Logging
# - warning levels:
@@ -268,7 +328,7 @@
}
croak("$LADCP_file: no valid data\n") unless ($nvv > 0);
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));
+ progress("\tattitude threshold (-t %d): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));
} else {
progress("Editing velocity data...\n");
$nvv = $cte = 0;
@@ -279,7 +339,7 @@
}
croak("$LADCP_file: no valid data\n") unless ($nvv > 0);
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));
+ progress("\tattitude threshold (-t %d): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));
}
#-------------------------------------------------------------------
@@ -297,6 +357,10 @@
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ for (my($beam)=0; $beam<4; $beam++) {
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam] *= $opt_x # HACK
+ if defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
+ }
($LADCP{ENSEMBLE}[$ens]->{U}[$bin],
$LADCP{ENSEMBLE}[$ens]->{V}[$bin],
$LADCP{ENSEMBLE}[$ens]->{W}[$bin],
@@ -331,6 +395,9 @@
$LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];
$LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];
$LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];
+ $LADCP{ENSEMBLE}[$ens]->{W}[$bin] *= $opt_x # HACK
+ if defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+
}
}
progress("\t$nvw valid velocities\n");
@@ -357,7 +424,7 @@
progress("Calculating LADCP time-series...\n");
($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) =
- calcLADCPts(\%LADCP,$refLr_firstBin,$refLr_lastBin,$opt_g);
+ calcLADCPts(\%LADCP,$opt_s,$refLr_firstBin,$refLr_lastBin,$opt_g);
croak("$LADCP_file: no good ensembles\n")
unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0);
@@ -390,6 +457,7 @@
progress("Reading CTD data...\n");
croak("$0: no CTD data\n") unless (&antsIn());
+undef($antsOldHeaders);
($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w,$CTD_w_t) =
&fnr('elapsed','depth','ss','w','w_t');
$CTD_temp = &fnrNoErr('temp');
@@ -456,11 +524,16 @@
$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});
+if ($opt_u) {
+ progress("\tskipping time lagging\n");
+} else {
+ $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});
+}
-$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});
#------------------------------------------------
@@ -478,10 +551,13 @@
next;
}
if ($skipped > 0) {
+ $firstGoodEns++,$skipped++,next # in gap
+ unless defined($LADCP{ENSEMBLE}[$firstGoodEns]->{REFLR_W});
info("$skipped initial LADCP ensembles skipped because CTD data begin with LADCP in water\n");
$skipped = 0;
}
if ($scan > $#{$CTD{ELAPSED}}) {
+ while (!defined($LADCP{ENSEMBLE}[$ens-1]->{REFLR_W})) { $ens--; } # in gap
info(sprintf("%d final LADCP ensembles skipped because CTD data end with LADCP in water\n",
$lastGoodEns-$ens+1));
$lastGoodEns = $ens-1;
@@ -494,7 +570,7 @@
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);
+ ) unless (abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2 + 1e-10);
$LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED} = $CTD{ELAPSED}[$scan]; # elapsed field for output
@@ -530,8 +606,8 @@
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);
+# 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));
@@ -608,14 +684,18 @@
progress("Binning velocities...\n");
+my($min_depth) = 9e99;
+my($max_depth) = 0;
+
progress("\tdowncast...\n");
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
-
- $LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin] =
+ $min_depth = $bindepth[$bin] if ($bindepth[$bin] < $min_depth);
+ $max_depth = $bindepth[$bin] if ($bindepth[$bin] > $max_depth);
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
@@ -626,27 +706,29 @@
push(@{$DNCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$DNCAST{BIN}[$bi]},$bin);
push(@{$DNCAST{DEPTH}[$bi]},$bindepth[$bin]);
- push(@{$DNCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin]);
+ push(@{$DNCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
- push(@{$DNCAST{PKGCORR_W}[100*round($CTD{W_T}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],0.01)+500]},
- $LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin])
- if defined($opt_p) && abs($CTD{W_T}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]) < 5;
+ # bin apparent ocean velocities as
+ # a function of package velocity
+# push(@{$DNCAST{PKGCORR_W}[10*round($CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],0.1)+50]},
+# $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin])
+# if defined($opt_p);
}
}
-if (defined($opt_x)) {
- progress("\t\tapplying surface-wave correction...\n");
- for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # first apply polynomial correction
- for (my($i)=0; $i<@{$DNCAST{W}[$bi]}; $i++) {
- for (my($e)=0; $e<@dc_corr_poly; $e++) {
- $DNCAST{W}[$bi][$i] -= $dc_corr_poly[$e] *
- $CTD{W_T}[$LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{CTD_SCAN}]**$e;
- }
- }
- }
-}
+#if (defined($opt_x)) { # apply polynomial package-velocity correction
+# progress("\t\tapplying package_velocity correction...\n");
+# for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) {
+# for (my($i)=0; $i<@{$DNCAST{W}[$bi]}; $i++) {
+# for (my($e)=0; $e<@pkg_vel_corr_poly; $e++) {
+# $DNCAST{W}[$bi][$i] -= $pkg_vel_corr_poly[$e] *
+# $CTD{W}[$LADCP{ENSEMBLE}[$DNCAST{ENSEMBLE}[$bi][$i]]->{CTD_SCAN}]**$e;
+# }
+# }
+# }
+#}
-for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) {
+for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) { # bin data
$DNCAST{MEAN_DEPTH}[$bi] = avg(@{$DNCAST{DEPTH}[$bi]});
$DNCAST{MEAN_ELAPSED}[$bi] = avg(@{$DNCAST{ELAPSED}[$bi]});
$DNCAST{MEDIAN_W}[$bi] = median(@{$DNCAST{W}[$bi]});
@@ -655,14 +737,16 @@
}
-progress("\tupcast...\n"); # upcast
+progress("\tupcast...\n"); # upcast
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
- $LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin] =
+ $min_depth = $bindepth[$bin] if ($bindepth[$bin] < $min_depth);
+ $max_depth = $bindepth[$bin] if ($bindepth[$bin] > $max_depth);
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
@@ -673,25 +757,25 @@
push(@{$UPCAST{CTD_W}[$bi]},$CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
push(@{$UPCAST{BIN}[$bi]},$bin);
push(@{$UPCAST{DEPTH}[$bi]},$bindepth[$bin]);
- push(@{$UPCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin]);
+ push(@{$UPCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin]);
- push(@{$UPCAST{PKGCORR_W}[100*round($CTD{W_T}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],0.01)+500]},
- $LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin])
- if defined($opt_p) && abs($CTD{W_T}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]) < 5;
+# push(@{$UPCAST{PKGCORR_W}[10*round($CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],0.1)+50]},
+# $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin])
+# if defined($opt_p);
}
}
-if (defined($opt_x)) {
- progress("\t\tapplying surface-wave correction...\n");
- for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) { # first apply polynomial correction
- for (my($i)=0; $i<@{$UPCAST{W}[$bi]}; $i++) {
- for (my($e)=0; $e<@uc_corr_poly; $e++) {
- $UPCAST{W}[$bi][$i] -= $uc_corr_poly[$e] *
- $CTD{W_T}[$LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{CTD_SCAN}]**$e;
- }
- }
- }
-}
+#if (defined($opt_x)) {
+# progress("\t\tapplying package-velocity correction...\n");
+# for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
+# for (my($i)=0; $i<@{$UPCAST{W}[$bi]}; $i++) {
+# for (my($e)=0; $e<@pkg_vel_corr_poly; $e++) {
+# $UPCAST{W}[$bi][$i] -= $pkg_vel_corr_poly[$e] *
+# $CTD{W}[$LADCP{ENSEMBLE}[$UPCAST{ENSEMBLE}[$bi][$i]]->{CTD_SCAN}]**$e;
+# }
+# }
+# }
+#}
for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
$UPCAST{MEAN_DEPTH}[$bi] = avg(@{$UPCAST{DEPTH}[$bi]});
@@ -701,6 +785,10 @@
$UPCAST{N_SAMP}[$bi] = @{$UPCAST{W}[$bi]};
}
+&antsAddParams('min_depth',$min_depth,'max_depth',$max_depth,
+ 'min_elapsed',$LADCP{ENSEMBLE}[$firstGoodEns]->{CTD_ELAPSED},
+ 'max_elapsed',$LADCP{ENSEMBLE}[$lastGoodEns]->{CTD_ELAPSED});
+
#--------------------------------------------------
# Calculate BT-referenced vertical-velocity profile
#--------------------------------------------------
@@ -712,7 +800,7 @@
my($sumSq) = my($n) = 0;
for (my($bi)=0; $bi<=$#{$BT{MEDIAN_W}}; $bi++) {
next unless defined($BT{MEDIAN_W}[$bi]);
- next unless ($BT{N_SAMP}[$bi]>=$opt_s && $DNCAST{N_SAMP}[$bi]>=$opt_s && $UPCAST{N_SAMP}[$bi]>=$opt_s);
+ next unless ($BT{N_SAMP}[$bi]>=$opt_k && $DNCAST{N_SAMP}[$bi]>=$opt_k && $UPCAST{N_SAMP}[$bi]>=$opt_k);
$sumSq += ($BT{MEDIAN_W}[$bi] - $DNCAST{MEDIAN_W}[$bi]/2 - $UPCAST{MEDIAN_W}[$bi]/2)**2;
$n++;
}
@@ -723,76 +811,149 @@
}
}
+#------------
+# full output
+#------------
+
+# NB: residual field is calculated with respect to down-/upcast medians in -o-size bins
+
+progress("Writing data...\n");
+
+@antsNewLayout = ('ensemble','bin','elapsed','depth','CTD_depth','downcast',
+ 'w','residual','CTD_w','LADCP_w','errvel',
+ 'pitch','roll','tilt','heading','3_beam','svel');
+
+for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ my(@bindepth) = calc_binDepths($ens);
+ for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ &antsOut(
+ $ens,$bin,$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},1,
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $DNCAST{MEDIAN_W}[$bi],
+ $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $LADCP{ENSEMBLE}[$ens]->{W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{PITCH},
+ $LADCP{ENSEMBLE}[$ens]->{ROLL},
+ $LADCP{ENSEMBLE}[$ens]->{TILT},
+ $LADCP{ENSEMBLE}[$ens]->{HEADING},
+ (defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + # only works for beam coords
+ defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
+ defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
+ defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
+ $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ );
+ }
+}
+
+for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) { # upcast
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ my(@bindepth) = calc_binDepths($ens);
+ for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ &antsOut(
+ $ens,$bin,$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $bindepth[$bin],$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},0,
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{SSCORRECTED_OCEAN_W}[$bin] - $UPCAST{MEDIAN_W}[$bi],
+ $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $LADCP{ENSEMBLE}[$ens]->{W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{PITCH},
+ $LADCP{ENSEMBLE}[$ens]->{ROLL},
+ $LADCP{ENSEMBLE}[$ens]->{TILT},
+ $LADCP{ENSEMBLE}[$ens]->{HEADING},
+ (defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]) + # only works for beam coords
+ defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]) +
+ defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]) +
+ defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3])) < 4 ? 1 : 0,
+ $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ );
+ }
+}
+
#---------------
# Output profile
#---------------
-progress("Writing vertical-velocity profile...\n");
-
-@antsNewLayout = ('depth','dc_depth','dc_elapsed','dc_w','dc_w.mad','dc_w.N',
- 'uc_depth','uc_elapsed','uc_w','uc_w.mad','uc_w.N',
- 'elapsed','w','w.mad','w.N',
- 'BT_w','BT_w.mad','BT_w.N');
+if (defined($opt_p)) {
+ progress("Writing vertical-velocity profiles to <$opt_p>...\n");
-for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) {
- &antsOut(($bi+0.5)*$opt_o, # nominal depth
- $DNCAST{MEAN_DEPTH}[$bi],$DNCAST{MEAN_ELAPSED}[$bi],
- $DNCAST{N_SAMP}[$bi]>=$opt_s?$DNCAST{MEDIAN_W}[$bi]:nan,
- $DNCAST{MAD_W}[$bi],$DNCAST{N_SAMP}[$bi],
- $UPCAST{MEAN_DEPTH}[$bi],$UPCAST{MEAN_ELAPSED}[$bi],
- $UPCAST{N_SAMP}[$bi]>=$opt_s?$UPCAST{MEDIAN_W}[$bi]:nan,
- $UPCAST{MAD_W}[$bi],$UPCAST{N_SAMP}[$bi],
- $DNCAST{MEAN_ELAPSED}[$bi]/2+$UPCAST{MEAN_ELAPSED}[$bi]/2,
- $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$opt_s ?
- ($DNCAST{MEDIAN_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MEDIAN_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
- nan,
- $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$opt_s ?
- ($DNCAST{MAD_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MAD_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
- nan,
- $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi],
- $BT{N_SAMP}[$bi]>=$opt_s?$BT{MEDIAN_W}[$bi]:nan,
- $BT{MAD_W}[$bi],$BT{N_SAMP}[$bi]
- );
-}
+ @antsNewLayout = ('depth','dc_depth','dc_elapsed','dc_w','dc_w.mad','dc_w.N',
+ 'uc_depth','uc_elapsed','uc_w','uc_w.mad','uc_w.N',
+ 'elapsed','w','w.mad','w.N',
+ 'BT_w','BT_w.mad','BT_w.N');
-#-------------------------
-# surface-wave effect file
-#-------------------------
-
-if (defined($opt_p)) {
- progress("Writing surface-wave-correction data to <$opt_p>...\n");
-
- @antsNewLayout = ('CTD_w_t','dc_w','dc_w.mad','dc_w.N','uc_w','uc_w.mad','uc_w.N','dc_w_corr','uc_w_corr');
&antsOut('EOF');
close(STDOUT);
open(STDOUT,">$opt_p") || croak("$opt_p: $!\n");
-
- for (my($bi)=0; $bi<=max($#{$DNCAST{PKGCORR_W}},$#{$UPCAST{PKGCORR_W}}); $bi++) {
- my($dc_N) = scalar(@{$DNCAST{PKGCORR_W}[$bi]});
- my($uc_N) = scalar(@{$UPCAST{PKGCORR_W}[$bi]});
- next unless ($dc_N>0 || $uc_N>0);
- my($dc_w) = median(@{$DNCAST{PKGCORR_W}[$bi]});
- my($uc_w) = median(@{$UPCAST{PKGCORR_W}[$bi]});
- my($w_t) = ($bi-500) / 100;
- if (defined($opt_x)) {
- my($dc_corr) = my($uc_corr) = 0;
- for (my($e)=0; $e<@dc_corr_poly; $e++) {
- $dc_corr += $dc_corr_poly[$e]*$w_t**$e;
- }
- for (my($e)=0; $e<@uc_corr_poly; $e++) {
- $uc_corr += $uc_corr_poly[$e]*$w_t**$e;
- }
- &antsOut($w_t,$dc_w,mad2($dc_w,@{$DNCAST{PKGCORR_W}[$bi]}),$dc_N,
- $uc_w,mad2($uc_w,@{$UPCAST{PKGCORR_W}[$bi]}),$uc_N,
- $dc_corr,$uc_corr);
- } else {
- &antsOut($w_t,$dc_w,mad2($dc_w,@{$DNCAST{PKGCORR_W}[$bi]}),$dc_N,
- $uc_w,mad2($uc_w,@{$UPCAST{PKGCORR_W}[$bi]}),$uc_N);
- }
+
+ for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) {
+ &antsOut(($bi+0.5)*$opt_o, # nominal depth
+ $DNCAST{MEAN_DEPTH}[$bi],$DNCAST{MEAN_ELAPSED}[$bi],
+ $DNCAST{N_SAMP}[$bi]>=$opt_k?$DNCAST{MEDIAN_W}[$bi]:nan,
+ $DNCAST{MAD_W}[$bi],$DNCAST{N_SAMP}[$bi],
+ $UPCAST{MEAN_DEPTH}[$bi],$UPCAST{MEAN_ELAPSED}[$bi],
+ $UPCAST{N_SAMP}[$bi]>=$opt_k?$UPCAST{MEDIAN_W}[$bi]:nan,
+ $UPCAST{MAD_W}[$bi],$UPCAST{N_SAMP}[$bi],
+ $DNCAST{MEAN_ELAPSED}[$bi]/2+$UPCAST{MEAN_ELAPSED}[$bi]/2,
+ $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$opt_k ?
+ ($DNCAST{MEDIAN_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MEDIAN_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
+ nan,
+ $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]>=$opt_k ?
+ ($DNCAST{MAD_W}[$bi]*$DNCAST{N_SAMP}[$bi]+$UPCAST{MAD_W}[$bi]*$UPCAST{N_SAMP}[$bi]) / ($DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi]) :
+ nan,
+ $DNCAST{N_SAMP}[$bi]+$UPCAST{N_SAMP}[$bi],
+ $BT{N_SAMP}[$bi]>=$opt_k?$BT{MEDIAN_W}[$bi]:nan,
+ $BT{MAD_W}[$bi],$BT{N_SAMP}[$bi]
+ );
}
}
+#-----------------------------
+# package-velocity effect file
+#-----------------------------
+
+#if (defined($opt_p)) {
+# progress("Writing package-velocity effect data to <$opt_p>...\n");
+#
+# @antsNewLayout = ('CTD_w','dc_w','dc_w.mad','dc_w.N','uc_w','uc_w.mad','uc_w.N','dc_w_corr','uc_w_corr');
+# &antsOut('EOF');
+#
+# close(STDOUT);
+# open(STDOUT,">$opt_p") || croak("$opt_p: $!\n");
+#
+# for (my($bi)=0; $bi<=max($#{$DNCAST{PKGCORR_W}},$#{$UPCAST{PKGCORR_W}}); $bi++) {
+# my($dc_N) = scalar(@{$DNCAST{PKGCORR_W}[$bi]});
+# my($uc_N) = scalar(@{$UPCAST{PKGCORR_W}[$bi]});
+# next unless ($dc_N>0 || $uc_N>0);
+# my($dc_w) = median(@{$DNCAST{PKGCORR_W}[$bi]});
+# my($uc_w) = median(@{$UPCAST{PKGCORR_W}[$bi]});
+# my($w) = ($bi-50) / 10;
+## if (defined($opt_x)) {
+## my($dc_corr) = my($uc_corr) = 0;
+## for (my($e)=0; $e<@pkg_vel_corr_poly; $e++) {
+## $dc_corr += $pkg_vel_corr_poly[$e]*$w**$e;
+## }
+## for (my($e)=0; $e<@pkg_vel_corr_poly; $e++) {
+## $uc_corr += $pkg_vel_corr_poly[$e]*$w**$e;
+## }
+## &antsOut($w,$dc_w,mad2($dc_w,@{$DNCAST{PKGCORR_W}[$bi]}),$dc_N,
+## $uc_w,mad2($uc_w,@{$UPCAST{PKGCORR_W}[$bi]}),$uc_N,
+## $dc_corr,$uc_corr);
+## } else {
+# &antsOut($w,$dc_w,mad2($dc_w,@{$DNCAST{PKGCORR_W}[$bi]}),$dc_N,
+# $uc_w,mad2($uc_w,@{$UPCAST{PKGCORR_W}[$bi]}),$uc_N);
+## }
+# }
+#}
+
#--------------------------------------
# write time-series output if requested
#--------------------------------------
@@ -833,71 +994,4 @@
undef($antsHeadersPrinted);
}
-#--------------------------------------------------------------------------------------------
-# Output all bins as separate files if requested
-# NB: - vertical LADCP velocities are corrected inaccurately for sound-speed variations!!!!
-# - full correction is used, on the other hand, for ocean velocities (w)
-#--------------------------------------------------------------------------------------------
-
-if (defined($opt_d)) {
-
- sub outProfBinRec($$$)
- {
- my($ens,$bin,$depth) = @_;
- my($sscorr) = $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]/1500;
-
- &antsOut($ens,
- $bin,
- $CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
- $depth,
- $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
- $LADCP{ENSEMBLE}[$ens]->{PITCH},
- $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH},
- $LADCP{ENSEMBLE}[$ens]->{ROLL},
- $LADCP{ENSEMBLE}[$ens]->{TILT},
- $LADCP{ENSEMBLE}[$ens]->{HEADING},
- $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
- $CTD{W_T}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
- $LADCP{ENSEMBLE}[$ens]->{W}[$bin]*$sscorr,
- $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin],
- $LADCP{ENSEMBLE}[$ens]->{REFLR_W},
- $LADCP{ENSEMBLE}[$ens]->{REFLR_W_ERR},
- $LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin]);
- }
-
- progress("Writing profile-bin data of downcast...\n");
-
- $commonParams = $antsCurParams;
- @antsNewLayout = ('ens','bin','elapsed','depth','sound_speed','pitch','gimbal_pitch',
- 'roll','tilt','heading','CTD_w','CTD_w_t','LADCP_w','LADCP_errvel',
- 'LADCP_reflr_w','LADCP_reflr_w_err','w');
-
- for (my($bi)=0; $bi<=$#{$DNCAST{ENSEMBLE}}; $bi++) {
- my($fn) = sprintf("$opt_d%03d.dncast",$bi);
- &antsOut('EOF');
- close(STDOUT);
- open(STDOUT,">$fn") || croak("$fn: $!\n");
- $antsCurParams = $commonParams;
- &antsAddParams('CTD_w',avg(@{$DNCAST{CTD_W}[$bi]}));
- for (my($eii)=0; $eii<=$#{$DNCAST{ENSEMBLE}[$bi]}; $eii++) {
- &outProfBinRec($DNCAST{ENSEMBLE}[$bi][$eii],$DNCAST{BIN}[$bi][$eii],$DNCAST{DEPTH}[$bi][$eii]);
- }
- }
-
- progress("Writing profile-bin data of upcast...\n");
-
- for (my($bi)=0; $bi<=$#{$UPCAST{ENSEMBLE}}; $bi++) {
- my($fn) = sprintf("$opt_d%03d.upcast",$bi);
- &antsOut('EOF');
- close(STDOUT);
- open(STDOUT,">$fn") || croak("$fn: $!\n");
- $antsCurParams = $commonParams;
- &antsAddParams('CTD_w',avg(@{$UPCAST{CTD_W}[$bi]}));
- for (my($eii)=0; $eii<=$#{$UPCAST{ENSEMBLE}[$bi]}; $eii++) {
- &outProfBinRec($UPCAST{ENSEMBLE}[$bi][$eii],$UPCAST{BIN}[$bi][$eii],$UPCAST{DEPTH}[$bi][$eii]);
- }
- close(STDOUT);
- }
-}
-
&antsExit();