.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/LADCP_VKE Sat Apr 04 07:22:55 2015 +0000
@@ -0,0 +1,214 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P _ V K E
+# doc: Tue Oct 14 11:05:16 2014
+# dlm: Fri Nov 7 13:54:06 2014
+# (c) 2012 A.M. Thurnherr
+# uE-Info: 33 15 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
+$antsMinLibVersion = 6.0;
+
+# HISTORY:
+# Oct 14, 2014: - created from [LADCPfs]
+# Oct 15, 2014: - added parameterization output
+# Oct 17, 2014: - changed parameterization constant $c to 0.021
+# Nov 6, 2014: - restored from backup and adapted to ANTS V6.0
+# Nov 7, 2014: - changed parameterization constant $c to 0.0215
+
+# NOTES:
+# - requires power densities
+# - output spectra (-s) have ADCP-related corrections applied unless -c is set
+
+($ANTSBIN) = ($0 =~ m{^(.*)/[^/]*$});
+($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
+require "$ANTSLIB/ants.pl";
+require "$ANTSLIB/libLADCP.pl";
+
+#----------------------------------------------------------------------
+# Usage
+#----------------------------------------------------------------------
+
+my($c) = 0.0215;
+
+&antsUsage('cne:r:s:',2,
+ '[suppress spectral -c)orrection] [tilt-correction -r)ange <max[0m]>]',
+ '[output -s)pectra <basename>]',
+ "[-e)ps-parameterization <scale[$c]>",
+ '<lambda.min> <lambda.max> [file...]');
+
+$lmin = &antsFloatArg(); # wavelength limits
+$lmax = &antsFloatArg();
+
+&antsFloatOpt(\$opt_e,$c); # default parameterization
+
+$imin = 0; # find frequency bin limits
+for ($nfreq=1; defined($P{"binpgrams::lambda.$nfreq"}); $nfreq++) {
+ $imin = $nfreq if ($imin==0 && $P{"binpgrams::lambda.$nfreq"}<=$lmax);
+ $imax = $nfreq if ($P{"binpgrams::lambda.$nfreq"} >= $lmin);
+}
+croak("$0: <lambda.min=$lmin> < min(lambda)")
+ unless defined($imax);
+
+$pg_fmin = fnr('pwrdens.0'); # first power field in spectra
+$fs_fmin = $pg_fmin + $imin; # first power field in finescale range
+$fs_fmax = $pg_fmin + $imax; # last power field in finescale range
+
+$mindf = fnr('depth.min');
+$maxdf = fnr('depth.max');
+
+#----------------------------------------------------------------------
+# Calculate Finescale Quantities (depending on input quantity)
+#----------------------------------------------------------------------
+
+sub integrate_fs_power($) # integrate fs spectrum (always done)
+{
+ my($r) = @_;
+
+ $ants_[$r][$fspwrf] = 0;
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
+ $ants_[$r][$fspwrf] += $ants_[0][$f];
+ }
+ $ants_[$r][$fspwrf] *= $P{'binpgrams::resolution_bandwidth'};
+}
+
+
+sub fit_universal_w_spec($) # vertical velocity => p0
+{
+ my($r) = @_;
+ my($nsamp) = $fs_fmax - $fs_fmin + 1;
+
+ #---------------------------------------------------
+ # fit slope-2 line in log-log space (main estimator)
+ #---------------------------------------------------
+
+ if ($nsamp >= 2) { # require min 2 wavenumber samples
+
+ my($sumd,$sumx,$sumy) = (0,0,0);
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { # loop over wavenumbers
+ my($i) = $f - $pg_fmin;
+ $sumd += log10($ants_[$r][$f]) + 2*log10($P{"binpgrams::k.$i"});
+ $sumx += log10($P{"binpgrams::k.$i"});
+ $sumy += log10($ants_[$r][$f]);
+ }
+ my($p0) = $sumd/$nsamp;
+ $ants_[$r][$p0f] = 10**$p0;
+
+ my($avgx) = $sumx/$nsamp; # avg for r calc
+ my($avgy) = $sumy/$nsamp;
+
+ my($sumsq,$sxx,$syy,$sxy) = (0,0,0,0); # calc rms misfit & correlation coeff
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
+ my($i) = $f - $pg_fmin;
+ my($xt) = log10($P{"binpgrams::k.$i"}) - $avgx;
+ my($yt) = log10($ants_[$r][$f]) - $avgy;
+ $sumsq += ($p0 - 2*log10($P{"binpgrams::k.$i"}) - log10($ants_[$r][$f]))**2;
+ $sxx += $xt * $xt; $syy += $yt * $yt; $sxy += $xt * $yt;
+ }
+ $ants_[$r][$rmsf] = sqrt($sumsq/$nsamp);
+ $ants_[$r][$rf] = $sxy/(sqrt($sxx * $syy) + $SMALL_AMOUNT);
+
+ } else {
+ &antsInfo("WARNING: no fit --- need min 2 samples");
+ }
+ $ants_[$r][$sampf] = $nsamp;
+}
+
+#----------------------------------------------------------------------
+# Process File
+#----------------------------------------------------------------------
+
+$fspwrf = &antsNewField('fs_pwr'); # derived fields
+
+if ($P{binpgrams::yfname} eq 'dc_w' || # w-finestructure parameterization
+ $P{binpgrams::yfname} eq 'uc_w') {
+ $p0f = &antsNewField('p0');
+ $rmsf = &antsNewField('p0fit.rms');
+ $sampf = &antsNewField('p0fit.nsamp');
+ $rf = &antsNewField('p0fit.r');
+ $wepsf = &antsNewField('eps.w');
+}
+
+my(@outLayout) = @antsNewLayout; # save for later
+for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input
+ $totf = $f if ($outLayout[$f] eq 'pwr.tot');
+ $tnsf = $f if ($outLayout[$f] eq 'pwr.tot.nsamp');
+}
+croak("$0: cannot find fields 'pwr.tot' or 'pwr.tot.nsamp' in input\n")
+ unless defined($totf) || defined($tnsf);
+$lsf = defined($tnsf) ? $tnsf : $totf;
+
+&antsInstallBufFull(0); # load entire file
+&antsIn();
+my($Hbuf) = $antsOldHeaders; # save for later
+
+for (my($r)=0; $r<@ants_; $r++) { # loop over all windows
+
+ #--------------------------
+ # apply spectral correction
+ #--------------------------
+
+ unless ($opt_c) {
+ for (my($i)=0; $i<$nfreq; $i++) { # loop over wavenumbers
+ if ($P{binpgrams::yfname} eq 'dc_w' || # vertical velocity
+ $P{binpgrams::yfname} eq 'uc_w') {
+ $ants_[$r][$i+$pg_fmin] *= T_w($P{"binpgrams::k.$i"},$P{ADCP_bin_length},$P{binpgrams::depth_resolution},$opt_r);
+ } elsif ($P{binpgrams::yfname} eq 'dc_w_depth'|| # dw/dz
+ $P{binpgrams::yfname} eq 'uc_w_depth') {
+ $ants_[$r][$i+$pg_fmin] *= T_w_z($P{"binpgrams::k.$i"},$P{ADCP_bin_length},$P{binpgrams::depth_resolution},$opt_r);
+ } else {
+ croak("$0: cannot calculate finescale $P{'binpgrams::yfname'} spectra (implementation restriction)\n");
+ }
+ }
+ }
+
+ #------------------------
+ # calculate fs quantities
+ #------------------------
+
+ integrate_fs_power($r); # always integrate spectra
+ fit_universal_w_spec($r)
+ if (($P{binpgrams::yfname} eq 'dc_w') || ($P{binpgrams::yfname} eq 'uc_w'));
+ $ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2;
+
+ #---------------
+ # produce output
+ #---------------
+
+ if (defined($opt_s)) { # output individual spectra
+ open(STDOUT_DUP,">&",STDOUT) || croak("$0: cannot dup STDOUT\n");
+ @antsNewLayout = ('k','lambda','pwrdens','finescale','pwrdens_fit');
+ $ofname = sprintf("${opt_s}_%02d.wspec",$r+1);
+ open(STDOUT,">$ofname") || croak("$ofname: $!\n");
+ undef($antsOldHeaders);
+ &antsActivateOut();
+
+ my($saveParams) = $antsCurParams; # add all extra input fields as %PARAMs
+ &antsAddParams('binpgrams::depth.min',$ants_[$r][$mindf],
+ 'binpgrams::depth.max',$ants_[$r][$maxdf]);
+ for (my($f)=$lsf+1; $f<@outLayout; $f++) {
+ &antsAddParams($outLayout[$f],$ants_[$r][$f]);
+ }
+
+ for (my($f)=$pg_fmin; $f<$pg_fmin+$nfreq; $f++) {
+ my($k) = eval(sprintf('$P{"binpgrams::k.%d"}',$f-$pg_fmin));
+ my($l) = eval(sprintf('$P{"binpgrams::lambda.%d"}',$f-$pg_fmin));
+ &antsOut($k,$l,$ants_[$r][$f],
+ ($f>=$fs_fmin && $f<=$fs_fmax),
+ numberp($ants_[$r][$p0f]) ? $ants_[$r][$p0f] * $k**(-2) : nan);
+ }
+
+ &antsOut('EOF');
+ open(STDOUT,">&",STDOUT_DUP) || croak("$0: cannot restore STDOUT\n");
+ close(STDOUT_DUP);
+ $antsCurParams = $saveParams;
+ }
+}
+
+@antsNewLayout = @outLayout;
+$antsOldHeaders = $Hbuf;
+$antsHeadersPrinted = 0;
+
+&antsFlush(); # output record with results
+&antsExit();
--- a/LADCP_w Mon Oct 27 13:36:21 2014 +0000
+++ b/LADCP_w Sat Apr 04 07:22:55 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W
# doc: Fri Dec 17 18:11:13 2010
-# dlm: Wed Oct 15 23:22:34 2014
+# dlm: Thu Mar 26 15:57:53 2015
# (c) 2010 A.M. Thurnherr
-# uE-Info: 232 46 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 988 33 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# TODO:
@@ -18,6 +18,7 @@
# ? use instrument tilt in sidelobe editing?
$antsSummary = 'calculate vertical velocities from LADCP & CTD time series';
+$antsMinLibVersion = 6.0;
# HISTORY:
# Dec 17, 2010: - created from [mergeCTD+LADCP]
@@ -145,14 +146,28 @@
# - BUG: soundspeed profile extrapolation code could bomb when seabed
# is shallower than profile depth
# - BUG: usage message had wrong -t default value
+# Oct 27, 2014: - removed dc/uc averaged w from profile output
+# Oct 30, 2014: - adapted to V6.0
+# Oct 31, 2014: - BUG: wrong -v default in usage message
+# - cosmetics
+# Nov 3, 2014: - removed unnecessary var @ARGS
+# - renamed CTD sound speed var to sspd
+# Nov 4, 2014: - BUG: PPI filtering was enabled by default for 300kHz instruments
+# - BUG: non-existing bottom was used to do sidelobe filtering
+# Nov 6, 2014: - added number of bins to diagnostic output
+# Nov 7, 2014: - added consistency check to make sure sound-speed profiles is complete
+# - BUG: soundspeed extrapolation to surface was not executed
+# Mar 19, 2015: - added %start_time from CTD params
+# Mar 26, 2015: - BUG: debug die() statement left in place
# CTD REQUIREMENTS
# - elapsed elapsed seconds; see note below
# - depth
-# - ss sound speed
+# - sspd sound speed
# - w ddepth/dt
# - temp OPTIONAL; used for backscatter calculation (i.e. not very important)
# - %lat/%lon OPTIONAL
+# - %start_time OPTIONAL
# 2-BEAM SOLUTIONS
# - for beam-coordinate data, two separate two-beam solutions (w12 & w34) are calculated
@@ -161,7 +176,7 @@
# NUMERICAL OPTIONS
# - the first option in the list cannot be numerical!
-# - if need be, use -v 1 as a dummy option
+# - if need be, use -v 2 as a dummy option
# ELAPSED TIMES
# - there are 2 different elapsed times used in this program:
@@ -221,11 +236,9 @@
# Usage
#------
-@ARGS = @ARGV; # save opts & arguments
-
$antsParseHeader = 0;
&antsUsage('3:4a:b:c:e:g:h:i:k:m:n:o:p:qs:t:uv:w:x:',1,
- '[-v)erbosity <level[1]>]',
+ '[-v)erbosity <level[2]>]',
'[-q)uick (no single-ping denoising)]',
'[require -4)-beam solutions] [apply beamvel-m)ask <file> if it exists]',
'[valid LADCP -b)ins <bin[2],bin[*]>',
@@ -239,7 +252,7 @@
'[pressure-sensor -a)cceleration correction <residual/CTD_w_tt>]',
'[-o)utput bin <resolution[10m]>] [-k) require <min[20]> samples]',
'[e-x)ecute <perl-expr>]',
- '<station-id> [run-label]');
+ '<profile id> [run-label]');
&antsUsageError() if ($opt_u && !defined($opt_i));
@@ -341,7 +354,7 @@
progress("Reading LADCP data from <$LADCP_file>...\n");
readData($LADCP_file,\%LADCP);
-progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
+progress("\t%d ensembles with %d bins each\n",scalar(@{$LADCP{ENSEMBLE}}),$LADCP{N_BINS});
croak("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -r\n")
unless ($LADCP{N_BINS} >= $refLr_lastBin);
@@ -616,10 +629,15 @@
open(STDIN,$CTD_file) || croak("$CTD_file: $!\n");
croak("$CTD_file: no data\n") unless (&antsIn());
undef($antsOldHeaders);
-($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w');
+
+($CTD_elapsed,$CTD_depth,$CTD_w) = &fnr('elapsed','depth','w'); # required fields
+$CTD_svel = &fnrNoErr('sspd'); # optional fields
+$CTD_svel = &fnr('ss') unless defined($CTD_svel); # backward compatibility
$CTD_temp = &fnrNoErr('temp');
-&antsAddParams('lat',$P{lat}) if defined($P{lat});
+
+&antsAddParams('lat',$P{lat}) if defined($P{lat}); # optional %PARAMs
&antsAddParams('lon',$P{lon}) if defined($P{lon});
+&antsAddParams('start_time',$P{start_time}) if defined($P{start_time});
$CTD_maxdepth = -1;
@@ -674,7 +692,7 @@
my($min_depth) = 9e99;
for (my($s)=0; $s<=$CTD_atbottom; $s+=$scans_per_sec) {
next unless ($CTD{DEPTH}[$s] >= 0 && numberp($CTD{SVEL}[$s]));
- $min_depth = $s if ($s < $min_depth);
+ $min_depth = $CTD{DEPTH}[$s] if ($CTD{DEPTH}[$s] < $min_depth);
$sVelProf[int($CTD{DEPTH}[$s])] = $CTD{SVEL}[$s];
}
while ($min_depth > 0) { # fill surface gap
@@ -685,8 +703,15 @@
$sVelProf[$d] = $sVelProf[$d-1]
unless defined($sVelProf[$d]);
}
+
+if (abs($#sVelProf - $CTD_maxdepth) > 10) {
+ croak(sprintf("$0: max CTD depth = %d m, bottom of sound-speed profile = %d m\n",
+ $CTD_maxdepth,$#sVelProf));
+} elsif (abs($#sVelProf - $CTD_maxdepth) > 1) {
+ warning(2,sprintf("max CTD depth = %d m, bottom of sound-speed profile = %d m",
+ $CTD_maxdepth,$#sVelProf));
+}
-
#-------------------
# Determine time lag
#-------------------
@@ -896,6 +921,8 @@
my($dd) = abs($water_depth_BT - $water_depth);
warning(2,sprintf("Large RDI vs. own water-depth difference (%.1fm)\n",$dd))
if ($dd > 5);
+ } else {
+ $water_depth = $sig_water_depth = undef;
}
}
@@ -915,7 +942,7 @@
($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,$water_depth);
progress("\t$nvrm velocities from $nerm ensembles removed\n");
- if ($PPI_editing) {
+ if (eval($PPI_editing_required)) {
progress("Editing data to remove PPI from seabed...\n");
progress("\tConstructing depth-average soundspeed profile...\n");
my($dz) = $water_depth - $#sVelProf; # $#sVelProf = max_depth(profile) in meters
@@ -958,7 +985,7 @@
#----------------------------------------------------------------------
# Data Editing after LADCP and CTD data have been merged
# 1) surface layer editing
-# 2) user-supplied $edit_data_hook
+# 2) user-supplied $post_merge_hook
#----------------------------------------------------------------------
progress("Removing data from instrument at surface...\n");
@@ -1279,6 +1306,7 @@
if (defined($out_w)) {
progress("Writing vertical-velocity data to $out_w...\n");
+ @k = keys(%P);
@antsNewLayout = ('ensemble','bin','elapsed','depth','CTD_depth','downcast',
'w','w12','w34','residual','CTD_w','CTD_w_tt','LADCP_w','errvel',
@@ -1367,36 +1395,27 @@
@antsNewLayout = ('depth','dc_depth','dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp','dc_w12','dc_w34',
'uc_depth','uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp','uc_w12','uc_w34',
- 'elapsed','w','w.mad','w.nsamp',
- 'BT_w','BT_w.mad','BT_w.nsamp');
+ 'BT_w','BT_w.mad','BT_w.nsamp');
open(STDOUT,"$out_profile") || croak("$out_profile: $!\n");
for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) {
- &antsOut(($bi+0.5)*$opt_o, # nominal depth
-# DOWNCAST
- $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],
- $DNCAST{N_SAMP}[$bi]>=$opt_k?$DNCAST{MEDIAN_W12}[$bi]:nan,
- $DNCAST{N_SAMP}[$bi]>=$opt_k?$DNCAST{MEDIAN_W34}[$bi]:nan,
-# UPCAST
- $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],
- $UPCAST{N_SAMP}[$bi]>=$opt_k?$UPCAST{MEDIAN_W12}[$bi]:nan,
- $UPCAST{N_SAMP}[$bi]>=$opt_k?$UPCAST{MEDIAN_W34}[$bi]:nan,
-# COMBINED
- $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]
+ &antsOut(($bi+0.5)*$opt_o, # nominal depth
+# DOWNCAST
+ $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],
+ $DNCAST{N_SAMP}[$bi]>=$opt_k?$DNCAST{MEDIAN_W12}[$bi]:nan,
+ $DNCAST{N_SAMP}[$bi]>=$opt_k?$DNCAST{MEDIAN_W34}[$bi]:nan,
+# UPCAST
+ $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],
+ $UPCAST{N_SAMP}[$bi]>=$opt_k?$UPCAST{MEDIAN_W12}[$bi]:nan,
+ $UPCAST{N_SAMP}[$bi]>=$opt_k?$UPCAST{MEDIAN_W34}[$bi]:nan,
+# BOTTOM TRACK
+ $BT{N_SAMP}[$bi]>=$opt_k?$BT{MEDIAN_W}[$bi]:nan,
+ $BT{MAD_W}[$bi],$BT{N_SAMP}[$bi]
);
}
&antsOut('EOF'); open(STDOUT,">&2");
--- a/LWplot_prof_2beam Mon Oct 27 13:36:21 2014 +0000
+++ b/LWplot_prof_2beam Sat Apr 04 07:22:55 2015 +0000
@@ -2,14 +2,15 @@
#======================================================================
# L W P L O T _ P R O F _ 2 B E A M
# doc: Fri Oct 14 09:42:36 2011
-# dlm: Thu Nov 21 10:22:15 2013
+# dlm: Mon Nov 3 22:01:45 2014
# (c) 2011 A.M. Thurnherr
-# uE-Info: 82 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 13 59 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# May 15, 2013: - created from [LWplot_prof]
# Oct 30, 2103: - got rid of non-portable echo -e
+# Nov 3, 2014: - adapted to updated layout of .prof file
# NOTES:
# - this version plots 2-beam solutions instead of final w
@@ -41,7 +42,7 @@
set -- $fields
[ "$1" = depth -a "$7" = dc_w12 -a "$8" = dc_w34 -a "$5" = dc_w.mad -a "$6" = dc_w.nsamp -a \
"${14}" = uc_w12 -a "${15}" = uc_w34 -a "${12}" = uc_w.mad -a "${13}" = uc_w.nsamp -a \
- "${20}" = BT_w -a "${21}" = BT_w.mad -a "${22}" = BT_w.nsamp ] || {
+ "${16}" = BT_w -a "${17}" = BT_w.mad -a "${18}" = BT_w.nsamp ] || {
echo "$0: file layout error" >&2
exit 1
}
@@ -85,17 +86,17 @@
awk '{print $8, $1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,coral,4_6:0 >> "$eps_file"
awk '{print $14,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,SeaGreen,6_2:0 >> "$eps_file"
awk '{print $15,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,SeaGreen,4_6:0 >> "$eps_file"
-awk '{print $20,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,black >> "$eps_file"
+awk '{print $16,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,black >> "$eps_file"
# MEAN ABSOLUTE DEVIATIONS (COMBINED SOLUTION)
awk '{print $5,$1, $4}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gcoral >> "$eps_file"
awk '{print $12,$1,$11}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -GSeaGreen >> "$eps_file"
-awk '{print $21,$1,$20}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gblack >> "$eps_file"
+awk '{print $17,$1,$20}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gblack >> "$eps_file"
# NUMBER OF SAMPLES (COMBINED SOLUTION)
awk '{print $6,$1, $4}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/coral >> "$eps_file"
awk '{print $13,$1,$11}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/SeaGreen >> "$eps_file"
-awk '{print $22,$1,$20}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/black >> "$eps_file"
+awk '{print $18,$1,$20}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/black >> "$eps_file"
# LABELS
echo 0.02 0.02 12 0 0 TL $out_basename $run_label | pstext -O -K $U $X >> "$eps_file"
--- a/PostProcess.sh Mon Oct 27 13:36:21 2014 +0000
+++ b/PostProcess.sh Sat Apr 04 07:22:55 2015 +0000
@@ -2,15 +2,16 @@
#======================================================================
# P O S T P R O C E S S . S H
# doc: Mon Oct 17 16:42:08 2011
-# dlm: Thu Oct 11 14:23:01 2012
+# dlm: Wed Nov 5 10:39:37 2014
# (c) 2011 A.M. Thurnherr
-# uE-Info: 13 47 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 14 62 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Oct 17, 2011: - created
# Oct 21, 2011: - made user friendly
# Oct 11, 2012: - added support for TL output
+# Nov 5, 2014: - changed extension of w samples output file
[ $# -gt 0 ] || {
echo "Usage: $0 <out-basename> [<run-label> [<data-subdir> [<plot-subdir> [<log-subdir>]]]]" >&2
@@ -23,7 +24,7 @@
[ -n "$4" ] && P_SUBDIR=$4 || P_SUBDIR=$RUN
[ -n "$5" ] && L_SUBDIR=$5 || L_SUBDIR=$RUN
-chmod +x $D_SUBDIR/$OBN.prof $D_SUBDIR/$OBN.w $D_SUBDIR/$OBN.TL
+chmod +x $D_SUBDIR/$OBN.prof $D_SUBDIR/$OBN.samp $D_SUBDIR/$OBN.TL
tile -s 1.3 -y -90 $P_SUBDIR/${OBN}_prof.eps $P_SUBDIR/${OBN}_w.eps \
$P_SUBDIR/${OBN}_BR.eps $P_SUBDIR/${OBN}_residuals.eps \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/SBE_w Sat Apr 04 07:22:55 2015 +0000
@@ -0,0 +1,364 @@
+#!/usr/bin/perl
+#======================================================================
+# S B E _ W
+# doc: Mon Nov 3 17:34:19 2014
+# dlm: Fri Nov 7 15:52:27 2014
+# (c) 2014 A.M. Thurnherr
+# uE-Info: 22 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
+$antsMinLibVersion = 6.0;
+
+# HISTORY:
+# Nov 3, 2014: - created
+# Nov 4, 2014: - improved
+# Nov 6, 2014: - BUG: sound speed was not calculated correctly
+# - added -a
+# - added conductivity & temperature editing
+# Nov 7, 2014: - loosened outlier editing
+# - added no-valid-data error message
+# - modified binning criterion to allow any sampling
+# frequency (not just divisors of 24)
+
+($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
+
+require "$ANTS/ants.pl";
+require "$ANTS/fft.pl";
+require "$ANTS/libSBE.pl";
+require "$ANTS/libEOS83.pl";
+
+$antsParseHeader = 0; # usage
+&antsUsage('al:rs:v:',1,
+ '[-v)erbosity <level[0]>]',
+ '[use -a)lternate sensor pair]',
+ '[-r)etain all data (no editing)]',
+ '[-s)ampling <rate[6Hz]>]',
+ '[-l)owpass w <cutoff[2s]>]',
+ '<SBE CNV file>');
+
+&antsFloatOpt(\$opt_l,2); # defaults
+&antsCardOpt(\$opt_s,6);
+&antsCardOpt(\$opt_v,0);
+
+$CNVfile = $ARGV[0]; # open CNV file
+open(F,&antsFileArg());
+&antsActivateOut(); # activate ANTS file
+
+#----------------------------------------------------------------------
+# Read Data
+#----------------------------------------------------------------------
+
+print(STDERR "Reading $CNVfile...") if ($opt_v);
+
+($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) = # decode SBE header
+ SBE_parseHeader(F,0,0); # SBE field names, no time check
+
+croak("$CNVfile: unexpected sampling interval $sampint\n")
+ unless (abs($sampint-1/24) < 1e-5);
+croak("$CNVfile: unknown latitude\n")
+ unless numberp($lat);
+
+$pressF = fnr('prDM');
+
+if ($opt_a) { # temp/cond alternate sensor pair
+ $tempF = fnr('t190C');
+ $condF = fnrNoErr('c1S/m');
+ if (defined($condF)) {
+ $condHistRes = 10; # 0.1 bins
+ } else {
+ $condF = fnr('c1mS/cm');
+ $condHistRes = 1; # 1.0 bins
+ }
+} else { # primary sensor pair
+ $tempF = fnr('t090C');
+ $condF = fnrNoErr('c0S/m');
+ if (defined($condF)) {
+ $condHistRes = 10;
+ } else {
+ $condF = fnr('c0mS/cm');
+ $condHistRes = 1;
+ }
+}
+
+&antsInstallBufFull(0); # read entire CNV file
+&SBEin(F,$ftype,$nfields,$nscans,$bad);
+
+printf(STDERR "\n\t%d scans",$nscans) if ($opt_v > 1);
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Edit Data
+# - pressure outliers & spikes
+# - conductivity outliers & spikes
+#----------------------------------------------------------------------
+
+unless ($opt_r) {
+ print(STDERR "Editing Data...") if ($opt_v);
+
+ #----------------------------------------
+ # trim initial records with nan pressures
+ #----------------------------------------
+ my($trimmed) = 0; # trim leading nan pressures
+ shift(@ants_),$trimmed++
+ until numberp($ants_[0][$pressF]);
+ printf(STDERR "\n\t%d initial records with nan pressure and/or conductivity trimmed",$trimmed) if ($opt_v > 1);
+ my($lvp) = $ants_[0][$pressF];
+
+ #------------------------------------------------
+ # edit pressure outliers outside contiguous range
+ #------------------------------------------------
+ my($outliers) = 0; my($min,$max);
+ for (my($s)=0; $s<$nscans; $s++) {
+ $phist[$ants_[$s][$pressF]+1000]++
+ if ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);
+ }
+ for ($max=25; $phist[$max+1000]||$phist[$max+1001]; $max++) {} # outliers after 2 gaps
+ for ($min=$max; $phist[$min+1000]||$phist[$min+ 999]; $min--) {}
+ for (my($s)=0; $s<$nscans; $s++) {
+ if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }
+ if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }
+ }
+ &antsAddParams("pressure_outliers",sprintf("%d",$outliers));
+ printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
+
+ #----------------------------------------------------
+ # edit conductivity outliers outside contiguous range
+ #----------------------------------------------------
+ $outliers = 0;
+ my($modeSamp)=0;
+ undef(@phist);
+ for (my($s)=0; $s<$nscans; $s++) {
+ next unless ($ants_[$s][$condF] > 0);
+ my($b) = $ants_[$s][$condF]*$condHistRes; # 1/10 S/m histogram resolution (1 mS/cm)
+ $phist[$b]++;
+ next unless ($phist[$b] > $modeSamp);
+ $modeSamp = $phist[$b]; $modeBin = $b;
+ }
+ for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {} # outliers after 2 gaps
+ for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}
+ $max /= $condHistRes; $min /= $condHistRes;
+ for (my($s)=0; $s<$nscans; $s++) {
+ if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; }
+ if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; }
+ }
+ &antsAddParams("conductivity_outliers",sprintf("%d",$outliers));
+ printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
+
+ #----------------------------------------------------
+ # edit temperature outliers outside contiguous range
+ #----------------------------------------------------
+ $outliers = 0;
+ my($modeSamp)=0;
+ undef(@phist);
+ for (my($s)=0; $s<$nscans; $s++) {
+ next unless ($ants_[$s][$tempF] > 0);
+ my($b) = $ants_[$s][$tempF]*10; # 10th of a degree histogram resolution
+ $phist[$b]++;
+ next unless ($phist[$b] > $modeSamp);
+ $modeSamp = $phist[$b]; $modeBin = $b;
+ }
+ for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {} # outliers after 2 gaps
+ for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}
+ $max /= 10; $min /= 10;
+ for (my($s)=0; $s<$nscans; $s++) {
+ if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; }
+ if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; }
+ }
+ &antsAddParams("temperature_outliers",sprintf("%d",$outliers));
+ printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
+
+ #----------------------------------------
+ # edit pressure spikes based on gradients
+ #----------------------------------------
+
+ for (my($s)=1; $s<$nscans; $s++) { # calculated pressure gradients (across gaps)
+ if (numberp($ants_[$s][$pressF])) {
+ $dp[$s-1] = $ants_[$s][$pressF] - $lvp;
+ $lvp = $ants_[$s][$pressF];
+ } else {
+ $dp[$s-1] = nan;
+ }
+ }
+
+ my($ns1,$ns2) = (0,0);
+ for (my($s)=0; $s<$nscans-2; $s++) { # consecutive large pressure gradients of opposite sign
+ if (($dp[$s]*$dp[$s+1] < 0) && # tests return false if either of the dps is not defined
+ (abs($dp[$s]) > 10) &&
+ (abs($dp[$s+1]) > 10)) {
+ $ants_[$s+1][$pressF] = nan;
+ $dp[$s] = $dp[$s+1] = undef;
+ $ns1++;
+ }
+ }
+ for (my($s)=0; $s<$nscans-3; $s++) { # 3 consecutive large pressure gradients of opposite sign
+ if (($dp[$s]>2 && $dp[$s+1]<-4 && $dp[$s+2]>2) ||
+ ($dp[$s]<-2 && $dp[$s+1]>4 && $dp[$s+2]<-2)) {
+ $ants_[$s+1][$pressF] = $ants_[$s+2][$pressF] = nan;
+ $dp[$s] = $dp[$s+1] = $dp[$s+2] = undef;
+ $ns2+=2;
+ }
+ }
+ &antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2));
+ printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v > 1);
+
+ printf(STDERR "\n") if ($opt_v);
+
+} # if $opt_r
+
+#----------------------------------------------------------------------
+# Correcting for pressure bias
+#----------------------------------------------------------------------
+
+print(STDERR "Correcting for pressure bias...") if ($opt_v);
+
+my($minP) = 9e99;
+for (my($s)=0; $s<$nscans; $s++) {
+ $minP = $ants_[$s][$pressF]
+ if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP);
+}
+croak("$CNVfile: no valid CTD pressure data below 25dbar\n")
+ unless ($minP < 9e99);
+&antsAddParams('pressure_bias',$minP);
+printf(STDERR "\n\tsubtracting %.1f dbar",$minP) if ($opt_v > 1);
+for (my($s)=0; $s<$nscans; $s++) {
+ $ants_[$s][$pressF] -= $minP
+ if numberp($ants_[$s][$pressF]);
+}
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Binning data
+#----------------------------------------------------------------------
+
+my($sps) = round(1 / $sampint / $opt_s);
+print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v);
+&antsAddParams('sampling_frequency',1/$opt_s);
+&antsAddParams('sampling_rate',$opt_s);
+
+my(@press,@temp,@cond);
+my($sp,$np,$st,$nt,$sc,$nc);
+
+$sp = $st = $sc = $np = $nt = $nc = 0;
+for (my($rec)=1,my($s)=0; $s<$nscans; $s++) {
+ if ($s*$sampint > $rec/$opt_s) {
+ $rec++;
+ push(@press,$np>0?$sp/$np:nan);
+ push(@temp, $nt>0?$st/$nt:nan);
+ push(@cond, $nc>0?$sc/$nc:nan);
+ $sp = $st = $sc = $np = $nt = $nc = 0;
+ }
+ $sp+=$ants_[$s][$pressF],$np++ if numberp($ants_[$s][$pressF]);
+ $st+=$ants_[$s][$tempF],$nt++ if numberp($ants_[$s][$tempF]);
+ $sc+=$ants_[$s][$condF],$nc++ if numberp($ants_[$s][$condF]);
+}
+
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Calculating derived quantities
+#----------------------------------------------------------------------
+
+print(STDERR "Calculating vertical package velocity & sound speed...") if ($opt_v);
+
+for (my($r)=0; $r<@press; $r++) {
+ $elapsed[$r] = $r/$opt_s;
+ $depth[$r] = &depth($press[$r],$lat);
+ croak(sprintf("$CNVfile: unrealistic depth %d m at elapsed=%.1f s (r=$r)\n",
+ $depth[$r],$elapsed[$r]))
+ if numberp($depth[$r]) && ($depth[$r]<0 || $depth[$r]>6100);
+ $sspd[$r] = &sVel(&salin($cond[$r],$temp[$r],$press[$r]),$temp[$r],$press[$r]);
+ croak(sprintf("$CNVfile: unrealistic soundspeed %dm/s at elapsed=%.1fs & depth=%.1fm ($cond[$r],$temp[$r],$press[$r])\n",
+ $sspd[$r],$elapsed[$r],$depth[$r]))
+ if numberp($sspd[$r]) && ($sspd[$r]<1400 || $sspd[$r]>1600);
+}
+
+$w[0] = nan;
+for (my($r)=1; $r<@depth-1; $r++) {
+ $w[$r] = numbersp($depth[$r-1],$depth[$r+1])
+ ? ($depth[$r+1] - $depth[$r-1]) * $opt_s
+ : nan;
+}
+push(@w,nan);
+
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Low-pass filter velocity data
+# - interpolate missing vertical velocities first
+#----------------------------------------------------------------------
+
+if ($opt_l > 0) {
+ print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v);
+ &antsAddParams('w_lowpass_cutoff',$opt_l);
+
+ my($trimmed) = 0;
+ shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++
+ until numberp($w[0]);
+ my($interpolated) = 0;
+ for ($r=1; $r<@w; $r++) {
+ next if numberp($w[$r]);
+ my($lv) = $r-1;
+ for ($nv=$r+1; $nv<@depth && !numberp($w[$nv]); $nv++) {}
+ if ($nv < @depth) {
+ while ($r < $nv) {
+ $w[$r] = $w[$lv] + ($r-$lv)/($nv-$lv) * ($w[$nv]-$w[$lv]);
+ $interpolated++;
+ $r++;
+ }
+
+ } else {
+ $trimmed += @w-$r;
+ splice(@w,$r); splice(@depth,$r);
+ splice(@elapsed,$r); splice(@sspd,$r);
+ }
+ }
+ &antsAddParams('w_interpolated',$interpolated);
+ printf(STDERR "\n\t%d/%d vertical velocities trimmed/interpolated",$trimmed,$interpolated) if ($opt_v > 1);
+
+
+#--------------------
+# Zero Pad Data
+#--------------------
+
+ for ($pot=1; $pot<@w; $pot<<=1) {} # determine power of two
+
+ for ($r=0; $r<@w; $r++) { # copy data
+ $fftbuf[2*$r] = $w[$r];
+ $fftbuf[2*$r+1] = 0;
+ }
+ printf(STDERR "\n\t%d zero records added",$pot-$r) if ($opt_v > 1);
+ while ($r < $pot) { # pad with zeroes
+ $fftbuf[2*$r] = $fftbuf[2*$r+1] = 0;
+ $r++;
+ }
+
+#--------------------
+# Low-Pass Filter
+#--------------------
+
+ @fco = &FOUR1(-1,@fftbuf); # forward FFT
+ $n = @fco/2;
+ for (my($ip)=2; $ip<=$n; $ip+=2) { # +ve freq fco
+ my($in) = 2*$n-$ip; # -ve freq fco
+ my($f) = $ip/2/$n*$opt_s; # frequency
+ $fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0
+ if ($f > 1/$opt_l); # low-pass filter
+ }
+ @w_lp = &FOUR1(1,@fco); # inverse FFT
+
+ printf(STDERR "\n") if ($opt_v);
+} else {
+ @w_lp = @w;
+}
+
+#----------------------------------------------------------------------
+
+print(STDERR "Writing output...\n") if ($opt_v);
+
+@antsNewLayout = ('elapsed','depth','sspd','w.raw','w');
+for ($r=0; $r<@w; $r++) {
+ &antsOut($elapsed[$r],$depth[$r],$sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp);
+}
+
+exit(0); # don't flush @ants_
--- a/Utilities/README Mon Oct 27 13:36:21 2014 +0000
+++ b/Utilities/README Sat Apr 04 07:22:55 2015 +0000
@@ -1,13 +1,21 @@
======================================================================
R E A D M E
doc: Wed Oct 17 11:57:40 2012
- dlm: Wed Oct 17 11:57:40 2012
+ dlm: Thu Mar 26 16:01:54 2015
(c) 2012 A.M. Thurnherr
- uE-Info: 5 38 NIL 0 0 72 0 2 4 NIL ofnI
+ uE-Info: 21 63 NIL 0 0 72 0 2 4 NIL ofnI
======================================================================
+=Usage=
+
In oder to use the utilities in this directory, include something like this:
require "$WCALC/Utilities/post_merge_TL_check.pl";
in the relevant ProcessingParams file.
+
+
+=List of Utilities=
+
+[post_merge_TL_check.pl] check time lagging (creates .TLcheck files)
+[post_merge_dwdz_filt.pl] filter ensembles with large |dw/dz|
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utilities/post_merge_dwdz_filt.pl Sat Apr 04 07:22:55 2015 +0000
@@ -0,0 +1,58 @@
+#======================================================================
+# P O S T _ M E R G E _ D W D Z _ F I L T . P L
+# doc: Thu Mar 26 16:02:09 2015
+# dlm: Thu Mar 26 17:11:24 2015
+# (c) 2015 A.M. Thurnherr
+# uE-Info: 10 25 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# Try to improve 2010 GoM data set by filtering ensembles with large dwdz
+# - no apparent benefit
+
+# HISTORY:
+# Mar 26, 2015: - created
+
+$post_merge_hook = sub { # arguments: firstGoodEns, lastGoodEns
+ my($fe,$le) = @_;
+
+ progress("\tediting ensembles with large dw/dz...\n");
+ my($nedt) = 0;
+
+ @antsNewLayout = ('ensemble','dwdz');
+ open(STDOUT,">$data_subdir/$out_basename.dwdz")
+ || croak("$data_subdir/$out_basename.dwdz: $!\n");
+
+ for (my($ens)=$fe; $ens<=$le; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); # skip non-valid
+
+ my($w1,$w2,$w3,$w4);
+ my($b1,$b2,$b3,$b4);
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ if (defined($w1) && defined($w2)) { # find last two
+ $w3 = $w4; $b3 = $b4;
+ $w4 = $LADCP{ENSEMBLE}[$ens]->{W}[$bin]; $b4 = $bin;
+ } else { # find first two
+ if (defined($w1)) { $w2 = $LADCP{ENSEMBLE}[$ens]->{W}[$bin]; $b2 = $bin; }
+ else { $w1 = $LADCP{ENSEMBLE}[$ens]->{W}[$bin]; $b1 = $bin; }
+ }
+ }
+
+ $nedt++,undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}),next # require at least 4 samples
+ unless defined($w1) && defined($w2) && defined($w3) && defined($w4);
+
+ my($dwdz) = (($w1+$w2)/2 - ($w3+$w4)/2) / (($b1+$b2)/2 - ($b3+$b4)/2)*$LADCP{BIN_LENGTH};
+ &antsOut($ens,$dwdz);
+ if (abs($dwdz) > 0.05) { # TWEAKABLE PARAMETER
+ undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
+ $nedt++,
+ next;
+ }
+
+ }
+
+ &antsOut('EOF'); open(STDOUT,'>&2');
+ progress("\t\t$nedt ensembles removed (%d%% of total)...\n",100*$nedt/($le-$fe+1));
+};
+
+1;
--- a/acoustic_backscatter.pl Mon Oct 27 13:36:21 2014 +0000
+++ b/acoustic_backscatter.pl Sat Apr 04 07:22:55 2015 +0000
@@ -1,9 +1,9 @@
#======================================================================
# A C O U S T I C _ B A C K S C A T T E R . P L
# doc: Wed Oct 20 13:02:27 2010
-# dlm: Thu Apr 17 08:49:21 2014
+# dlm: Fri Nov 7 14:45:44 2014
# (c) 2010 A.M. Thurnherr
-# uE-Info: 19 34 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 80 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -16,7 +16,8 @@
# - SV now saved in ensemble
# Oct 21, 2011: - BUG: made code work for uplooker again
# Mar 4, 2014: - added support for missing PITCH/ROLL (TILT)
-# Apr 17, 2017: - BUG: missing ;
+# Apr 17, 2014: - BUG: missing ;
+# Nov 7, 2014: - BUG: calc_binDepths() was called without valid CTD depth
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Deines (IEEE 1999)
@@ -75,6 +76,7 @@
my($cosBeamAngle) = cos(rad($LADCP{BEAM_ANGLE}));
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
my(@bd) = calc_binDepths($ens);
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
my($depth) = int($bd[$bin]);
--- a/defaults.pl Mon Oct 27 13:36:21 2014 +0000
+++ b/defaults.pl Sat Apr 04 07:22:55 2015 +0000
@@ -1,9 +1,9 @@
#======================================================================
# D E F A U L T S . P L
# doc: Tue Oct 11 17:11:21 2011
-# dlm: Wed Oct 15 23:21:48 2014
+# dlm: Tue Nov 4 10:35:07 2014
# (c) 2011 A.M. Thurnherr
-# uE-Info: 39 68 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 44 63 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -37,27 +37,16 @@
# May 20, 2014: - added support for $PPI_editing
# May 21, 2014: - added $PPI_extend_upper_limit
# Oct 15, 2014: - investigated, modified and documented -t default
-
-# Variable Names:
-# - variables that are only used in a particular library are
-# prefixed with a 2-caps code
+# Oct 27, 2014: - removed CTD-acceleration-effects and w spectral plots
+# - removed time-lagging stats output file by default
+# Oct 31, 2014: - re-arranged order of things
+# - .w => .samp output
+# Nov 4, 2014: - BUG: PPI_editing did not work as advertised
#======================================================================
# Data Input
#======================================================================
-# File to load cruise/cast-specific processing parameters from
-
-if (-r "ProcessingParams.$RUN") {
- $processing_param_file = "ProcessingParams.$RUN";
-} elsif (-r "ProcessingParams.default") {
- $processing_param_file = "ProcessingParams.default";
-} elsif (-r "ProcessingParams") {
- $processing_param_file = "ProcessingParams";
-} else {
- croak("$0: cannot find either <ProcessingParams.$RUN> or <ProcessingParams[.default]>\n");
-}
-
# CTD depth adjustment
# - set with -d (-a up to 2013/05/16)
# - value is added to CTD pressure
@@ -116,46 +105,6 @@
$data_subdir = $plot_subdir = $log_subdir = $RUN;
-# main w output and all its plots:
-# _w.eps vertical velocities
-# _residuals.eps residual vertical velocities
-# _Sv.eps volume scattering coefficient after Deimes (1999)
-# _corr.eps correlation [DISABLED 2013/05/16]
-
-$out_w = "| LWplot_residuals $plot_subdir/${out_basename}_residuals.eps" .
- "| LWplot_Sv $plot_subdir/${out_basename}_Sv.eps" .
-# "| LWplot_corr $plot_subdir/${out_basename}_corr.eps" .
- "| LWplot_w $plot_subdir/${out_basename}_w.eps" .
- "> $data_subdir/$out_basename.w";
-
-
-# w profile output
-
-$out_profile = "| LWplot_prof_2beam $plot_subdir/${out_basename}_prof.eps" .
- "| LWplot_spec $plot_subdir/${out_basename}_spec.eps" .
- "> $data_subdir/$out_basename.prof";
-
-# log output
-
-$out_log = "$log_subdir/$out_basename.log";
-
-
-# time-series output (CTD acceleration effect)
-
-$out_timeseries = "| LWplot_CAE $plot_subdir/${out_basename}_CAE.eps" .
- "> $data_subdir/$out_basename.tis";
-
-
-# per-bin residual output (plot only)
-
-$out_BR = "| LWplot_BR $plot_subdir/${out_basename}_BR.eps";
-
-
-# time-lagging output
-
-$out_TL = "| LWplot_TL $plot_subdir/${out_basename}_TL.eps" .
- "> $data_subdir/$out_basename.TL";
-
#======================================================================
# Data Editing
#======================================================================
@@ -233,6 +182,8 @@
# PPI editing as described in [edit_data.pl]
# - enabled by default for WH150 data
+# - variable sets an expression, to be evaluated once the data
+# have been loaded
# - 2014 CLIVAR P16 #47 has a slight discontinuity at 4000m; this
# discontinuity is there without PPI filtering but gets slightly
# worse with PPI filtering. Setting $PPI_extend_upper_limit to
@@ -243,7 +194,7 @@
# set by the shortest acoustic path between the ADCP and the
# seabed.
-$PPI_editing = ($LADCP{BEAM_FREQUENCY} < 300);
+$PPI_editing_required = '($LADCP{BEAM_FREQUENCY} < 300)';
#$PPI_extend_upper_limit = 1.03; # arbitrarily increase calculated max dist from seabed by 3%
@@ -337,5 +288,84 @@
$BT_max_w_error = 0.03;
#======================================================================
+# ProcessingParams file selection
+#======================================================================
+
+if (-r "ProcessingParams.$RUN") {
+ $processing_param_file = "ProcessingParams.$RUN";
+} elsif (-r "ProcessingParams.default") {
+ $processing_param_file = "ProcessingParams.default";
+} elsif (-r "ProcessingParams") {
+ $processing_param_file = "ProcessingParams";
+} else {
+ croak("$0: cannot find either <ProcessingParams.$RUN> or <ProcessingParams[.default]>\n");
+}
+
+#----------------------------------------------------------------------
+# Processing log (diagnostic messages) output
+#----------------------------------------------------------------------
+
+$out_log = "$log_subdir/$out_basename.log";
+
+
+#----------------------------------------------------------------------
+# Vertical-velocity profile output and plots:
+# Data:
+# *.prof vertical velocity profiles
+# Standard Plots:
+# *_prof.eps vertical velocity profiles (main output plot)
+# Optional Plots:
+# *_wspec.eps vertical-velocity wavenumber spectra
+#----------------------------------------------------------------------
+
+$out_profile = "| LWplot_prof_2beam $plot_subdir/${out_basename}_prof.eps" .
+# "| LWplot_spec $plot_subdir/${out_basename}_spec.eps" .
+ "> $data_subdir/$out_basename.prof";
+
+
+#----------------------------------------------------------------------
+# Vertical-velocity sample data output and plots:
+# Data:
+# *.samp w sample data
+# Standard Plots:
+# *_w.eps vertical velocity time-depth plot
+# *_residuals.eps residual vertical velocity time-depth plot
+# *_Sv.eps volume scattering coefficient time-depth plot
+# Optional Plots:
+# *_corr.eps correlation time-depth plot [REMOVED FROM DEFAULTS 2013/05/16]
+#----------------------------------------------------------------------
+
+$out_w = "| LWplot_residuals $plot_subdir/${out_basename}_residuals.eps" .
+ "| LWplot_Sv $plot_subdir/${out_basename}_Sv.eps" .
+# "| LWplot_corr $plot_subdir/${out_basename}_corr.eps" .
+ "| LWplot_w $plot_subdir/${out_basename}_w.eps" .
+ "> $data_subdir/$out_basename.samp";
+
+
+#----------------------------------------------------------------------
+# Time-series output
+# Data:
+# *.tis combined CTD/LADCP time-series data, including
+# package- and LADCP reference layer w
+# Optional Plots:
+# *_CAE.eps plot of CTD acceleration effects on reference-layer w
+#----------------------------------------------------------------------
+
+$out_timeseries =
+# "| LWplot_CAE $plot_subdir/${out_basename}_CAE.eps" .
+ "> $data_subdir/$out_basename.tis";
+
+#----------------------------------------------------------------------
+# Per-bin vertical-velocity residuals (plot only)
+#----------------------------------------------------------------------
+
+$out_BR = "| LWplot_BR $plot_subdir/${out_basename}_BR.eps";
+
+
+#----------------------------------------------------------------------
+# Time-lagging correlation statistics (plot only)
+#----------------------------------------------------------------------
+
+$out_TL = "| LWplot_TL $plot_subdir/${out_basename}_TL.eps";
1; # return true