--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/LADCP_w Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,1037 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P _ W
+# doc: Fri Dec 17 18:11:13 2010
+# dlm: Sun Jul 3 12:58:19 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 844 45 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?
+# 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';
+
+# HISTORY:
+# Dec 17, 2010: - created from [mergeCTD+LADCP]
+# Dec 18, 2010: - made to work
+# Dec 19, 2010: - improved considerably
+# Dec 20, 2010: - onward
+# - BUG: depth-binning was off by 1 bin?!
+# - added binning correction for instrument tilt
+# Dec 21, 2010: - added -h (seafloor depth)
+# Dec 22, 2010: - BUG: had not applied soundspeed-correction to w
+# - debugged opt_d
+# Dec 23, 2010: - continued implementation of soundspeed corrections
+# Dec 24, 2010: - added winch_w, wave_w
+# - removed beampair velocities from code
+# Dec 25, 2010: - adapted for surface-wave correction in terms of acceleration (CTD_w_t)
+# - removed elapsed_mismatch
+# - removed winch_w, wave_w
+# Dec 26, 2010: - made -p output layout independent of -x to avoid Makefile problems
+# Dec 30, 2010: - cleaned up some
+# - folded-in backscatter calculation from shear method
+# - folded-in BT calculation from shear method
+# Dec 31, 2010: - added weighted mean w profile to output
+# Jan 2, 2010: - BUG: BT_w.mad could bomb with division by 0
+# - 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
+
+# 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
+# - depth
+# - ss sound speed
+# - w ddepth/dt
+# - w_t dw/dt
+# - temp OPTIONAL; used for backscatter calculation (i.e. not very important)
+
+# NUMERICAL OPTIONS
+# - the first option in the list cannot be numerical!
+# - if need be, use -v 1 as a dummy option
+
+# ELAPSED TIMES
+# - there are 2 different elapsed times used in this program:
+# 1) elapsed based on firstgoodens in the LADCP time series
+# 2) CTD elapsed time
+# - CTD elapsed time does not have to start with zero!
+# - do not use the Seabird elapsed field, which is only reported to
+# 3 significant digits, causing significant jitter in dt; however,
+# at least up to 2010 Seabird simply calculates elapsed time by
+# assuming a 24Hz sampling rate and no record drop; therefore,
+# it is best to calculate elapsed time as %RECNO/24
+# - the elapsed field of the output is the elapsed time from the CTD
+# file; this is required in order to be able to compare the times
+# from the uplooker and downlooker-derived vertical velocity
+# profiles
+# - as a result, a profile only starts with elapsed==0 if the CTD
+# is turned on when the LADCP is already in the water
+
+# TIME LAGGING
+# - occasionally, the time lagging algorithm fails, in particular if the CTD is turned on
+# some time after the package enters the water
+# - in this case, an initial guess can be provided with -i
+# - e.g. plot 'CTD/24Hz/054.1Hz elapsed w','LADCP/raw/054UL.prof =$elapsed-850 w' => -i -850
+
+# VELOCITY AMBIGUITY ERRORS
+# - quite extensive tests with DIMES US2 station 146, which has a lot of
+# ambiguity velocity errors, reveal that -m 1 catches those errors
+# quite nicely
+# - even when the errors are not filtered with -m 1, they do not
+# affect the w profiles, as long as the median bin values are used
+
+# SCREEN LOGGING
+# - there are 4 verbosity levels, selected by -v
+# 0: only print errors
+# 1: default, UNIX-like (warnings and info messages that are not produced for every cast)
+# 2: progress messages and useful information
+# >2: debug messges
+# - the most useful ones of these are 1 & 2
+
+($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
+($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
+($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
+
+require "$ANTS/ants.pl";
+require "$ANTS/libstats.pl";
+require "$WCALC/edit_data.pl";
+require "$WCALC/time_series.pl";
+require "$WCALC/time_lag.pl";
+require "$WCALC/find_seabed.pl";
+require "$WCALC/svel_corrections.pl";
+require "$WCALC/acoustic_backscatter.pl";
+require "$WCALC/bottom_tracking.pl";
+require "$PERL_TOOLS/RDI_BB_Read.pl";
+require "$PERL_TOOLS/RDI_Coords.pl";
+
+@ARGS = @ARGV; # save opts & arguments
+
+$antsParseHeader = 0;
+&antsUsage('03:4a:b:c:d:e:f:g:h:i:l:m:n:o: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]>]',
+ '[-c)orrelation <min[70]>] [-t)ilt <max[10deg]> [-e)rr-vel <max[0.1m/s]>]',
+ '[-h water <depth>]',
+ '[max LADCP time-series -g)ap <length[60s]>]',
+ '[-m)ax vertical <velocity[1m/s]>',
+ '[-a)djust CTD depth <by[0m]>]',
+ '[-i)nitial CTD time offset <guestimate> [-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 <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-velocity effect <file>]',
+ '[-f) write time-series <file>] [output residual -d)iagnostics to <file>]',
+ '<LADCP-file> [CTD-file]');
+
+&antsCardOpt(\$opt_v,1); # suppress regular info
+
+$RDI_Coords::minValidVels = 4 if ($opt_4); # suppress 3-beam solutions
+
+&antsFloatOpt(\$opt_c,70); # min correlation
+&antsFloatOpt(\$opt_t,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));
+
+$opt_n = '10,100' unless defined($opt_n); # number of time-lags to carry out
+$opt_w = '240,20' unless defined($opt_w); # time-lag search window (full width)
+&antsFloatOpt(\$opt_3,0.6);
+
+&antsFloatOpt(\$opt_o,50); # output bin size
+&antsCardOpt(\$opt_s,20); # min samples required
+
+$opt_r = '2,6' unless defined($opt_r); # reference layer bins for w for time matching
+$opt_b = '2,*' unless defined($opt_b); # bins to use in w calculations
+
+@n_lags = split(',',$opt_n); # decode -n
+croak("$0: cannot decode -n $opt_n\n")
+ unless numberp($n_lags[0]) && numberp($n_lags[1]);
+
+@w_size = split(',',$opt_w); # decode -w
+croak("$0: cannot decode -w $opt_w\n")
+ unless numberp($w_size[0]) && numberp($w_size[1]);
+
+($refLr_firstBin,$refLr_lastBin) = split(',',$opt_r); # decode -r
+croak("$0: cannot decode -r $opt_r\n")
+ unless numberp($refLr_firstBin) && numberp($refLr_lastBin);
+
+($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b); # decode -b
+croak("$0: cannot decode -b $opt_b\n")
+ unless (numberp($LADCP_firstBin) &&
+ ($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));
+
+#if (defined($opt_x)) { # decode corrections
+# (@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]);
+#}
+
+&antsFloatOpt(\$opt_x,1);
+
+#if (defined($opt_d)) { # make sure output directory is clean
+# croak("$0: old depth-bin files <${opt_d}[0-9][0-9][0-9].dncast> found --- remove before creating new ones!\n")
+# if (glob("${opt_d}[0-9][0-9][0-9].dncast"));
+# croak("$0: old depth-bin files <${opt_d}[0-9][0-9][0-9].upcast> found --- remove before creating new ones!\n")
+# if (glob("${opt_d}[0-9][0-9][0-9].upcast"));
+#}
+
+$LADCP_file = &antsFileArg();
+
+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:
+# 0 probably unimportant, e.g. nonsensical parameters that probably won't affect solution
+# 1 may be somewhat important
+# 2 important
+#----------------------------------------------------------------------
+
+sub progress(@)
+{ printf(STDERR @_) if ($opt_v > 1); }
+
+sub info(@)
+{
+ print(STDERR ($opt_v > 1) ? "\t" : "$LADCP_file: ");
+ printf(STDERR @_) if ($opt_v > 0);
+}
+
+sub warning(@)
+{
+ my($lvl,@msg) = @_;
+ return if ($opt_v == 0);
+ print(STDERR "\n") if ($opt_v > 1);
+ print(STDERR "WARNING (level $lvl): ");
+ printf(STDERR @msg);
+ print(STDERR "\n") if ($opt_v > 1);
+}
+
+sub debugmsg(@)
+{ printf(STDERR @_) if ($opt_v > 2); }
+
+#----------------
+# Read LADCP data
+#----------------
+
+progress("Reading LADCP data ($LADCP_file)...\n");
+readData($LADCP_file,\%LADCP);
+progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
+
+croak("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -r\n")
+ unless ($LADCP{N_BINS} >= $refLr_lastBin);
+
+croak("$0: first reference-layer bin outside valid range\n")
+ unless ($refLr_firstBin>=1 && $refLr_firstBin<=$LADCP{N_BINS});
+croak("$0: last reference-layer bin outside valid range\n")
+ unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS});
+croak("$0: first reference-layer bin > last reference-layer bin\n")
+ unless ($refLr_firstBin <= $refLr_lastBin);
+
+$LADCP_lastBin = $LADCP{N_BINS}-1
+ if ($LADCP_lastBin eq '*');
+croak("$0: first valid LADCP bin outside valid range\n")
+ unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS});
+croak("$0: last valid LADCP bin outside valid range\n")
+ unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS});
+croak("$0: first valid LADCP bin > last valid LADCP bin\n")
+ unless ($LADCP_firstBin <= $LADCP_lastBin);
+
+warning(0,"first reference-layer bin < first valid LADCP bin\n")
+ unless ($refLr_firstBin >= $LADCP_firstBin);
+warning(0,"last reference-layer bin > last valid LADCP bin\n")
+ unless ($refLr_lastBin <= $LADCP_lastBin);
+
+warning(1,"if at all, bin 1 should not be used for short blank-after-transmit values\n")
+ if ($LADCP{BLANKING_DISTANCE}<$LADCP{BIN_LENGTH} && $refLr_firstBin==1);
+
+#------------------------------------------------------------
+# Edit beam-velocity data
+# 1) correlation threshold
+# 2) tilt threshold (also sets TILT field in all ensembles)
+#------------------------------------------------------------
+
+if ($LADCP{BEAM_COORDINATES}) {
+ progress("Editing beam-velocity data...\n");
+ $nvv = $cte = 0;
+ for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $nvv += countValidBeamVels($ens);
+ $cte += editCorr($ens,$opt_c);
+ $pte += editTilt($ens,$opt_t);
+ }
+ 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 (-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;
+ for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $nvv += countValidBeamVels($ens);
+ $cte += editCorr_Earthcoords($ens,$opt_c);
+ $pte += editTilt($ens,$opt_t);
+ }
+ 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 (-t %d): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));
+}
+
+#-------------------------------------------------------------------
+# Calculate earth velocities
+# - this is done for all bins (not just valid ones), to allow
+# useless possibility that invalid bins are used for reflr calcs
+# - also calculate separate beam-pair velocities
+# - the UNEDITED velocities are saved for the BT calculations
+# (W is required, U & V are only used for stats that have not
+# been very useful so far)
+#-------------------------------------------------------------------
+
+if ($LADCP{BEAM_COORDINATES}) {
+ progress("Calculating earth-coordinate velocities...\n");
+ $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],
+ $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) =
+ velInstrumentToEarth(\%LADCP,$ens,
+ velBeamToInstrument(\%LADCP,
+ @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));
+ $nvw += defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+ gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+ $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];
+ }
+ }
+ progress("\t$nvw valid velocities\n");
+ progress("\t3-beam solutions : $RDI_Coords::threeBeam_1 " .
+ "$RDI_Coords::threeBeam_2 " .
+ "$RDI_Coords::threeBeam_3 " .
+ "$RDI_Coords::threeBeam_4\n")
+ unless ($opt_4);
+} else {
+ progress("Counting valid vertical velocities...\n");
+ $nvw = 0;
+ for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ ($LADCP{ENSEMBLE}[$ens]->{U}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{V}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{W}[$bin],
+ $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) = @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]};
+ $nvw += defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ $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");
+}
+
+#----------------------------------------------
+# S1 STEP: Edit earth-coordinate -velocity data
+# 1) error-velocity threshold
+#----------------------------------------------
+
+progress("Editing earth-coordinate velocity data...\n");
+
+$ete = 0;
+for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $ete += editErrVel($ens,$opt_e);
+}
+progress("\t error-velocity threshold (-e %.1f): %d velocites removed (%d%% of total)\n",
+ $opt_e,$ete,round(100*$ete/$nvw));
+
+#----------------------------
+# Calculate LADCP time series
+#----------------------------
+
+progress("Calculating LADCP time-series...\n");
+
+($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) =
+ calcLADCPts(\%LADCP,$refLr_firstBin,$refLr_lastBin,$opt_g);
+croak("$LADCP_file: no good ensembles\n")
+ unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0);
+
+my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} -
+ $LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED};
+croak("$0: implausibly short cast ($cast_duration seconds)\n")
+ unless ($cast_duration > 600);
+
+$LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1);
+
+progress("\tStart of cast : %s (#%5d)\n",
+ $LADCP{ENSEMBLE}[$firstGoodEns]->{TIME},
+ $LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER});
+progress("\tBottom of cast : %s (#%5d) @ dz~%.1fm\n",
+ $LADCP{ENSEMBLE}[$LADCP_atbottom]->{TIME},
+ $LADCP{ENSEMBLE}[$LADCP_atbottom]->{NUMBER},
+ $LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH});
+progress("\tEnd of cast : %s (#%5d)\n",
+ $LADCP{ENSEMBLE}[$lastGoodEns]->{TIME},
+ $LADCP{ENSEMBLE}[$lastGoodEns]->{NUMBER});
+progress("\tCast duration : %.1f hours (pinging for %.1f hours)\n",
+ $cast_duration / 3600,
+ ($LADCP{ENSEMBLE}[$#{$LADCP{ENSEMBLE}}]->{UNIX_TIME} -
+ $LADCP{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);
+progress("\tMean ping interval: %.1f seconds\n",$LADCP{MEAN_DT});
+
+#--------------
+# Read CTD data
+#--------------
+
+progress("Reading CTD data...\n");
+croak("$0: no CTD data\n") unless (&antsIn());
+undef($antsOldHeaders);
+($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w,$CTD_w_t) =
+ &fnr('elapsed','depth','ss','w','w_t');
+$CTD_temp = &fnrNoErr('temp');
+
+$CTD_maxdepth = -1;
+
+do {
+ croak("$0: cannot deal with non-numeric CTD elapsed time\n")
+ unless &antsNumbers($CTD_elapsed);
+ push(@{$CTD{ELAPSED}},$ants_[0][$CTD_elapsed]);
+ push(@{$CTD{DEPTH}}, $ants_[0][$CTD_depth]+$opt_a);
+ push(@{$CTD{SVEL}}, $ants_[0][$CTD_svel]);
+ push(@{$CTD{W}}, $ants_[0][$CTD_w]);
+ push(@{$CTD{W_T}}, $ants_[0][$CTD_w_t]);
+ push(@{$CTD{TEMP}}, $ants_[0][$CTD_temp]) if defined($CTD_temp);
+ if ($ants_[0][$CTD_depth]+$opt_a > $CTD_maxdepth) {
+ $CTD_maxdepth = $ants_[0][$CTD_depth]+$opt_a;
+ $CTD_atbottom = $#{$CTD{DEPTH}};
+ }
+} while (&antsIn());
+
+$CTD{DT} = ($CTD{ELAPSED}[$#{$CTD{ELAPSED}}] - $CTD{ELAPSED}[0]) / $#{$CTD{ELAPSED}};
+
+progress("\t%d scans at %.1fHz\n",scalar(@{$CTD{DEPTH}}),1/$CTD{DT});
+progress("\tstart depth = %.1fm\n",$CTD{DEPTH}[0]);
+croak("$0: CTD start depth must be numeric\n")
+ unless numberp($CTD{DEPTH}[0]);
+progress("\tmax depth = %.1fm (# $CTD_atbottom)\n",$CTD_maxdepth);
+
+#--------------------------------------------------------------------
+# Construct sound-speed correction profile from CTD 1Hz downcast data
+# very simple algorithm that stores the last value found
+# in each 1m bin
+#--------------------------------------------------------------------
+
+progress("Constructing sound-speed correction profile\n");
+
+my($scans_per_sec) = int(1/$CTD{DT}+0.5);
+for (my($s)=0; $s<=$CTD_atbottom; $s+=$scans_per_sec) {
+ next unless ($CTD{DEPTH}[$s] >= 0 && numberp($CTD{SVEL}[$s]));
+ $sVelProf[int($CTD{DEPTH}[$s])] = $CTD{SVEL}[$s];
+}
+
+#-------------------
+# Determine time lag
+#-------------------
+
+if (defined($opt_i)) {
+ progress("Setting initial time lag...\n");
+ $CTD{TIME_LAG} = $opt_i;
+ progress("\t-i => elapsed(CTD) ~ elapsed(LADCP) + %.1fs\n",$CTD{TIME_LAG});
+} else {
+ progress("Guestimating time lag...\n");
+
+ my($CTD_10pct_down) = my($LADCP_10pct_down) = 0;
+ $CTD_10pct_down++
+ until ($CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0] >= 0.1*($CTD_maxdepth-$CTD{DEPTH}[0]));
+ $LADCP_10pct_down++
+ until ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} >= 0.1*$LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH});
+
+ $CTD{TIME_LAG} = $CTD{ELAPSED}[$CTD_10pct_down] - $LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED};
+
+ progress("\telapsed(dz(CTD)=%.1fm) ~ elapsed(dz(LADCP)=%.1fm) + %.1fs\n",
+ $CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH},$CTD{TIME_LAG});
+}
+
+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});
+}
+
+&antsAddParams('CTD_time_lag',$CTD{TIME_LAG});
+
+#------------------------------------------------
+# Merge CTD with LADCP data
+# - after this step, reflr w is sound-speed corrected!!!
+#------------------------------------------------
+
+progress("Merging CTD with LADCP data...\n");
+
+for (my($skipped)=0,my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+ my($scan) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
+ if ($scan <= 0) { # NB: must be <=, rather than <, or assertion below sometimes fails
+ $skipped++;
+ $firstGoodEns = $ens+1;
+ next;
+ }
+ if ($skipped > 0) {
+ info("$skipped initial LADCP ensembles skipped because CTD data begin with LADCP in water\n");
+ $skipped = 0;
+ }
+ if ($scan > $#{$CTD{ELAPSED}}) {
+ info(sprintf("%d final LADCP ensembles skipped because CTD data end with LADCP in water\n",
+ $lastGoodEns-$ens+1));
+ $lastGoodEns = $ens-1;
+ last;
+ }
+
+ die("assertion failed!\n" .
+ "\ttest: abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2\n" .
+ "\tens = $ens, scan = $scan\n" .
+ sprintf("\tadjusted LADCP time = %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG}) .
+ sprintf("\tCTD($scan) time = %f\n",$CTD{ELAPSED}[$scan]) .
+ "=> Did you use SeaBird elapsed time? Don't!"
+ ) unless (abs($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$scan]) <= $CTD{DT}/2);
+
+ $LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED} = $CTD{ELAPSED}[$scan]; # elapsed field for output
+
+ if (defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W}) # not a gap
+ && numberp($CTD{DEPTH}[$scan])) {
+ $LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $CTD{SVEL}[$scan]/1500; # correct (inadequately) for sound-speed variations
+ croak(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs (use -a?)\n",
+ $CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan]))
+ unless ($CTD{DEPTH}[$scan] >= 0);
+ $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan];
+ $LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan;
+ my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan];
+ if (abs($reflr_ocean_w) <= $opt_m) {
+ $sumWsq += &SQR($reflr_ocean_w);
+ $nWsq++;
+ if ($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} > 100 &&
+ $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} < $LADCP{ENSEMBLE}[$LADCP_atbottom]->{CTD_DEPTH}-100) {
+ $sumWsqI += &SQR($reflr_ocean_w);
+ $nWsqI++;
+ }
+ } else {
+ undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); # DON'T USE THIS ENSEMBLE LATER
+ }
+ } else{
+ undef($LADCP{ENSEMBLE}[$ens]->{REFLR_W}); # don't output in time-series file
+ undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); # old DEPTH from calcLADCPts()
+ }
+}
+
+if ($nWsq > 0 && $nWsqI > 0) {
+ &antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI));
+ progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n",
+ 100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI));
+ warning(0,"%.2f cm/s reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI))
+ if (sqrt($sumWsqI/$nWsqI) > 0.05);
+# croak("$0: rms reference-layer w_ocean is too large\n")
+# unless (sqrt($sumWsqI/$nWsqI) < 0.07);
+} elsif ($nWsq > 0) {
+ &antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',nan);
+ progress("\t%.2f cm/s rms reference-layer w_ocean\n",100*sqrt($sumWsq/$nWsq));
+} else {
+ croak("$0: no valid vertical velocities\n");
+}
+
+#----------------------------------------------------------------------------
+# Remove data contaminated by sidelobe reflection from seabed and sea surface
+#----------------------------------------------------------------------------
+
+if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {
+ &antsAddParams('ADCP_orientation','downlooker');
+
+ if (numberp($opt_h)) {
+ progress("Setting water depth (-h)\n");
+ $water_depth = $opt_h;
+ $sig_water_depth = 0;
+ } else {
+ progress("Finding seabed...\n");
+ calc_backscatter_profs($firstGoodEns,$lastGoodEns);
+ ($water_depth,$sig_water_depth) =
+ find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_atbottom]->{CTD_DEPTH});
+ ($water_depth_BT,$sig_water_depth_BT) =
+ find_seabed(\%LADCP,$LADCP_atbottom,$LADCP{BEAM_COORDINATES});
+ if (defined($water_depth_BT)) {
+ my($dd) = abs($water_depth_BT - $water_depth);
+ warning(2,sprintf("Large RDI vs. own water-depth difference (%.1fm)\n",$dd))
+ if ($dd > 5);
+ }
+ }
+
+ &antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth);
+
+ if (defined($water_depth)) {
+ if (defined($water_depth_BT)) {
+ progress("\t%.1f(%.1f) m water depth (%.1f(%.1f)m from BT)\n",
+ $water_depth,$sig_water_depth,$water_depth_BT,$sig_water_depth_BT);
+ } else {
+ progress("\t%.1f(%.1f) m water depth\n",$water_depth,$sig_water_depth);
+ }
+ warning(1,sprintf("large uncertainty in water-depth estimation (%.1fm)\n",$sig_water_depth))
+ if ($sig_water_depth > $LADCP{BIN_LENGTH});
+
+ progress("Editing data to remove sidelobe interference from seabed...\n");
+ ($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,$water_depth);
+ progress("\t$nvrm velocities from $nerm ensembles removed\n");
+
+ } else {
+ info("no seabed found in backscatter profiles --- no sidelobe editing done\n");
+ }
+
+} else {
+ &antsAddParams('ADCP_orientation','uplooker');
+ progress("Editing data to remove sidelobe interference from sea surface...\n");
+ ($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,undef);
+ progress("\t$nvrm velocities from $nerm ensembles removed\n");
+}
+
+#---------------------------------------------------------------------------
+# Depth-bin LADCP velocity data
+#
+# NOTES:
+# 1) ensemble and bin numbers are saved for maximum flexibility
+# 2) only ensemble/bins with valid vertical velocities are saved
+# 3) applying the full soundspeed correction to w is most likely pointless in
+# practice, but hey!, CPU cycles are cheap; [in a cast in the Gulf of Mexico
+# which has fairly pronounce soundspeed gradients, the max value of Kn
+# is 1.00004160558372, which gives rise to a correction of less than 0.2mm/s
+# at a winch+wave speed of 3m/s....]
+# 4) as far as I can tell, the soundspeed correction for bin length also
+# has only a minute effect
+#---------------------------------------------------------------------------
+
+progress("Binning velocities...\n");
+
+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]->{SSCORRECTED_OCEAN_W}[$bin] =
+ sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
+ $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
+ $bindepth[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ push(@{$DNCAST{ENSEMBLE}[$bi]},$ens);
+ push(@{$DNCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
+ 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]->{SSCORRECTED_OCEAN_W}[$bin]);
+
+ # 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)) { # 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++) { # 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]});
+ $DNCAST{MAD_W}[$bi] = mad2($DNCAST{MEDIAN_W}[$bi],@{$DNCAST{W}[$bi]});
+ $DNCAST{N_SAMP}[$bi] = @{$DNCAST{W}[$bi]};
+}
+
+
+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]->{SSCORRECTED_OCEAN_W}[$bin] =
+ sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
+ $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH},
+ $bindepth[$bin]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ push(@{$UPCAST{ENSEMBLE}[$bi]},$ens);
+ push(@{$UPCAST{ELAPSED}[$bi]},$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]);
+ 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]->{SSCORRECTED_OCEAN_W}[$bin]);
+
+# 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 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]});
+ $UPCAST{MEAN_ELAPSED}[$bi] = avg(@{$UPCAST{ELAPSED}[$bi]});
+ $UPCAST{MEDIAN_W}[$bi] = median(@{$UPCAST{W}[$bi]});
+ $UPCAST{MAD_W}[$bi] = mad2($UPCAST{MEDIAN_W}[$bi],@{$UPCAST{W}[$bi]});
+ $UPCAST{N_SAMP}[$bi] = @{$UPCAST{W}[$bi]};
+}
+
+#--------------------------------------------------
+# Calculate BT-referenced vertical-velocity profile
+#--------------------------------------------------
+
+if (defined($water_depth)) {
+ progress("Calculating BT-referenced vertical velocities\n");
+ calc_BTprof($firstGoodEns,$lastGoodEns,$water_depth,$sig_water_depth);
+
+ 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);
+ $sumSq += ($BT{MEDIAN_W}[$bi] - $DNCAST{MEDIAN_W}[$bi]/2 - $UPCAST{MEDIAN_W}[$bi]/2)**2;
+ $n++;
+ }
+ if ($n > 0) {
+ my($rms) = round(sqrt($sumSq/$n),0.001);
+ &antsAddParams('BT_rms_w_discrepancy',$rms);
+ progress("\t$rms m/s rms vertical-velocity discrepancy\n");
+ }
+}
+
+#---------------
+# 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');
+
+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]
+ );
+}
+
+#-------------------------------------
+# diagnostic output with all residuals
+#-------------------------------------
+
+if (defined($opt_d)) {
+ progress("Writing residuals-diagnostic data to <$opt_d>...\n");
+
+ @antsNewLayout = ('ensemble','bin','elapsed','depth','downcast',
+ 'w','residual','CTD_w','LADCP_w','errvel',
+ 'pitch','roll','tilt','heading','3_beam','svel');
+ &antsOut('EOF');
+
+ close(STDOUT);
+ open(STDOUT,">$opt_d") || croak("$opt_d: $!\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]);
+ my($bi) = $bindepth[$bin]/$opt_o;
+ &antsOut(
+ $ens,$bin,$CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $bindepth[$bin],$ens<=$LADCP_atbottom,
+ $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,
+ $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]->{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,
+ $CTD{SVEL}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ );
+ }
+ }
+
+}
+
+
+#-----------------------------
+# 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
+#--------------------------------------
+
+if (defined($opt_f)) {
+ progress("Writing time-series data to <$opt_f>...\n");
+
+ @antsNewLayout = ('ens','elapsed',
+ 'depth','sound_speed','pitch','gimbal_pitch','roll','tilt','heading',
+ 'CTD_w','CTD_w_t','LADCP_reflr_w','LADCP_reflr_w_err',
+ 'ocean_reflr_w');
+ &antsOut('EOF');
+
+ close(STDOUT);
+ open(STDOUT,">$opt_f") || croak("$opt_f: $!\n");
+
+ for ($ens=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
+ my($reflr_w) = defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W})
+ ? $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}]
+ : undef;
+ &antsOut($ens,
+ $CTD{ELAPSED}[$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN}],
+ $LADCP{ENSEMBLE}[$ens]->{CTD_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]->{REFLR_W},
+ $LADCP{ENSEMBLE}[$ens]->{REFLR_W_ERR},
+ $reflr_w);
+ }
+
+ close(STDOUT);
+ 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]->{SSCORRECTED_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();
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/README Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,24 @@
+======================================================================
+ R E A D M E
+ doc: Sun Jul 3 13:28:47 2011
+ dlm: Sun Jul 3 13:28:47 2011
+ (c) 2011 A.M. Thurnherr
+ uE-Info: 16 70 NIL 0 0 72 3 2 4 NIL ofnI
+======================================================================
+
+During 2011_IWISE I began playing with Visbeck-type plots, which
+required different output. Additionally, based on processing a single
+in-ice station, I now think that there may be beam-velocity bias. I am
+not sure, however, because while multiplying the beam velocities with a
+factor of 1.03 appears to take care of the down-/up-cast bias, there is
+still a correlation between the residuals and CTD_w. In any case, I
+changed the semantics of the old -x option, that was supposed to
+correct for package-motion effects and removed the -p option. I also
+removed the old -d option that dumped the bins into individual files,
+and replaced it with a "residual diagnostics" output (for the
+Visbeck-type plots). All the old code is left in place, commented out.
+I soon realized that the main profile output can be easily and exactly
+reconstructed from the "residual diagnostic output" using bindata, i.e.
+it seems more useful to have the diagnostic output become the default.
+To make the code cleaner, I want to remove all the commented-out stuff.
+Therefore I saved the current version so I can return to it if needed.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/acoustic_backscatter.pl Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,149 @@
+#======================================================================
+# 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 Dec 30 22:22:02 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 131 36 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Oct 20, 2010: - created
+# Dec 10, 2010: - BUG: backscatter above sea surface made code bomb
+# when run with uplooker data
+# Dec 30, 2010: - adapted for use with [LADCP_w]
+
+#----------------------------------------------------------------------
+# Volume Scattering Coefficient, following Deines (IEEE 1999)
+# NOTES:
+# - instrument specific! (300kHz Workhorse)
+# - no sound-speed correction applied
+# - R in bin center, instead of last quarter
+# - transmit power assumes 33V batteries
+#----------------------------------------------------------------------
+
+# NB:
+# - correction seems to work for a subset of bins (~bins 3-9 for
+# 2010 P403 station 46)
+# - this may imply that noise level depends on bin
+# - far bins are important for seabed detection, i.e. cannot simply
+# be discarded at this stage
+
+sub Sv($$$$$)
+{
+ my($temp,$PL,$Er,$R,$EA) = @_;
+ my($C) = -143; # RDI Workhorse monitor
+ my($Ldbm) = 10 * log10($PL);
+ my($PdbW) = 14.0;
+ my($alpha) = 0.069;
+ my($Kc) = 0.45;
+
+ return $C + 10*log10(($temp+273)*$R**2) - $Ldbm - $PdbW
+ + 2*$alpha*$R + $Kc*($EA-$Er);
+}
+
+#----------------------------------------------------------------------
+# Calculate per-bin backscatter profiles
+# input: first and last valid LADCP ensemble
+# output: sSv[$depth][$bin] sum of volume scattering coefficients
+# nSv[$depth][$bin] number of samples in bin
+#----------------------------------------------------------------------
+
+my(@sSv,@nSv);
+
+sub calc_backscatter_profs($$)
+{
+ my($LADCP_start,$LADCP_end) = @_;
+
+ my(@Er) = (1e99,1e99,1e99,1e99); # echo intensity reference level
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ $Er[0] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][0] < $Er[0]);
+ $Er[1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][1] < $Er[1]);
+ $Er[2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][2] < $Er[2]);
+ $Er[3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][3] < $Er[3]);
+ }
+ debugmsg("per-beam noise levels = @Er\n");
+
+ my($cosBeamAngle) = cos(rad($LADCP{BEAM_ANGLE}));
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ my(@bd) = calc_binDepths($ens);
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ my($depth) = int($bd[$bin]);
+# next if ($depth < 0); # enable to use this code for uplookers
+ my($range_to_bin) = ($bd[$bin] - $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH})
+ / cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))
+ / $cosBeamAngle;
+ if (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP})) {
+ $sSv[$depth][$bin] += Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[0],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[1],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[2],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[3],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;
+ } else {
+ $sSv[$depth][$bin] += Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[0],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[1],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[2],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMPERATURE},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[3],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;
+ }
+ $nSv[$depth][$bin]++;
+ }
+ }
+}
+
+#----------------------------------------------------------------------
+# determine location of seabed from backscatter profiles
+# input: depth below seabed can possibly be (e.g. max CTD depth)
+# output: median/mad of estimated water depth
+#----------------------------------------------------------------------
+
+sub find_backscatter_seabed($)
+{
+ my($search_below) = int($_[0]); # grid index to begin search
+ my(@wdepth); # list of water_depth indices
+
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { # find backscatter min/max below $search_below in each bin
+ my($minSv,$maxSv,$depthmaxSv,$lastvalid) = (1e99,-1e99,-1,-1);
+ for (my($depth)=$search_below; $depth<@nSv; $depth++) {
+ next unless ($nSv[$depth][$bin] > 0);
+ my($Sv) = $sSv[$depth][$bin] / $nSv[$depth][$bin];
+ $lastvalid = $depth;
+ $minSv = $Sv if ($Sv < $minSv);
+ $maxSv = $Sv, $depthmaxSv = $depth if ($Sv > $maxSv);
+ }
+ if ($maxSv-$minSv>10 && $depthmax!=$lastvalid) { # ignore boundary maxima & scatter
+ push(@wdepth,$depthmaxSv);
+ }
+ }
+
+ my($wd) = median(@wdepth);
+ return ($wd,mad2($wd,@wdepth));
+
+}
+
+1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/bottom_tracking.pl Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,168 @@
+#======================================================================
+# B O T T O M _ T R A C K I N G . P L
+# doc: Wed Oct 20 21:05:37 2010
+# dlm: Fri Dec 31 14:21:14 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 32 30 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Oct 20, 2010: - created
+# Dec 30, 2010: - adapted for use with LADCP_w
+
+# This code is essentially identical to the one used in LADCPproc. Differences:
+# 1) velocity editing is simpler: no wake editing, no PPI editing, no shear
+# editing, no w outlier
+# 2) median/mad calculated instead of mean/stddev
+# 3) u,v not calculated
+
+# Don't look for BT-referenced velocities if package is more than $BT_max_range
+# above seabed. This parameter is frequency dependent and the current value is
+# appropriate (if rather high) for 300kHz Workhorse intruments.
+
+my($BT_max_range) = 300;
+
+
+# The code only tries to bin BT-referenced velocities if a consistent bottom
+# is available in all 4 beams. Ensembles where the range of bin numbers where
+# the maximum echo is found is greater than $max_BIT_bin_range_diff are rejected.
+# In addition to flukes this also rejects ensembles collected with large
+# instrument tilts. The value of 3 is a first guess that has not been explored.
+
+my($BT_max_bin_range_diff) = 3;
+
+
+# If the difference between measured vertical velocity of the seabed (i.e.
+# the package vertical velocity referenced by the seabed) and the vertical
+# velocity of the CTD (from dp/dt) si greater than $BT_max_w_error the current
+# ensemble is ignored and $nBTwFlag is increased. The value of
+# 3cm/s is taken from listBT developed on A0304 cruise.
+
+my($BT_max_w_error) = 0.03;
+
+#----------------------------------------------------------------------
+# bin valid BT-referenced velocities
+# input: ensemble number, water depth (with uncertainty)
+# output: @BTu,@BTv,@BTw main result
+# editflags for information
+# @BTbtmu, @BTbtmv, @BTtilt stats to find reasons for quality differences
+#----------------------------------------------------------------------
+
+my($nBTfound,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0);
+my(@BTu,@BTv,@BTw);
+
+sub binBTprof($$$)
+{
+ my($ens,$wd,$sig_wd) = @_;
+
+ my(@ea_max) = (0,0,0,0); my(@ea_max_bin) = (nan,nan,nan,nan);
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ $ea_max[0] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0],
+ $ea_max_bin[0] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0] > $ea_max[0]);
+ $ea_max[1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1],
+ $ea_max_bin[1] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1] > $ea_max[1]);
+ $ea_max[2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2],
+ $ea_max_bin[2] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2] > $ea_max[2]);
+ $ea_max[3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3],
+ $ea_max_bin[3] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3] > $ea_max[3]);
+ }
+
+ return # disregard boundary maxima
+ unless (min(@ea_max_bin) > $LADCP_firstBin-1) &&
+ (max(@ea_max_bin) < $LADCP_lastBin-1);
+ return # inconsistent range to seabed
+ unless (max(@ea_max_bin)-min(@ea_max_bin) <= $BT_max_bin_range_diff);
+
+ $nBTfound++;
+ my($seafloor_bin) = round(avg(@ea_max_bin));
+
+ my(@bd) = calc_binDepths($ens);
+ $nBTdepthFlag++,return # BT range inconsistent with water depth
+ unless (abs($wd-$bd[$seafloor_bin]) < max($sig_wd,$LADCP{BIN_LENGTH}));
+
+ # try vertical velocities at seabed bin plus one above and below
+ # this does not really work because, often, only one of the bins has valid velocities
+ my($w1) = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin-1];
+ my($w2) = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin ];
+ my($w3) = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin+1];
+
+ $w1 = 9e99 unless numberp($w1); # invalid velocity sentinels
+ $w2 = 9e99 unless numberp($w1);
+ $w3 = 9e99 unless numberp($w1);
+
+ my($seafloor_u,$seafloor_v,$seafloor_w);
+
+ # determine which of the three trial bins is most consistent with reflr vertical velocities
+ die("assertion failed") unless numberp($LADCP{ENSEMBLE}[$ens]->{REFLR_W});
+ if (abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w2) &&
+ abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w3)) {
+ $seafloor_u = $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$seafloor_bin-1];
+ $seafloor_v = $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$seafloor_bin-1];
+ $seafloor_w = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin-1];
+ } elsif (abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{REFLR_W}-$w2)) {
+ $seafloor_u = $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$seafloor_bin+1];
+ $seafloor_v = $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$seafloor_bin+1];
+ $seafloor_w = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin+1];
+ } else {
+ $nBTvalidVelFlag++,return # none of 3 trial bins has valid velocity
+ if ($w2 == 9e99);
+ $seafloor_u = $LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$seafloor_bin];
+ $seafloor_v = $LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$seafloor_bin];
+ $seafloor_w = $LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$seafloor_bin];
+ }
+
+ $nBTwFlag++,return # velocity error is too great
+ if (abs($seafloor_w-$LADCP{ENSEMBLE}[$ens]->{REFLR_W}) > $BT_max_w_error);
+
+ push(@BTbtmu,$seafloor_u); # calc stats to try to find reasons for quality
+ push(@BTbtmv,$seafloor_v);
+ push(@BTbtmw,$seafloor_w);
+ push(@BTtilt,$LADCP{ENSEMBLE}[$ens]->{TILT});
+
+ for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
+ next unless defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ my($gi) = int($bd[$bin]) / $opt_o;
+ push(@{$BTw[$gi]},$LADCP{ENSEMBLE}[$ens]->{W}[$bin]-$seafloor_w);
+ }
+}
+
+#----------------------------------------------------------------------
+# calculate BT-referenced velocity profile
+# input: start,end LADCP ensembles, water depth with uncertainty
+# output: %BT{MEDIAN_W,MAD_W,N_SAMP}
+#----------------------------------------------------------------------
+
+sub calc_BTprof($$$$)
+{
+ my($LADCP_start,$LADCP_end,$wd,$sig_wd) = @_;
+
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ next unless ($wd-$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} < $BT_max_range);
+ binBTprof($ens,$wd,$sig_wd);
+ }
+
+ progress("\t$nBTfound BT ensembles found\n");
+ progress("\t$nBTdepthFlag flagged bad because of wrong bottom depth\n");
+ progress("\t$nBTvalidVelFlag flagged bad because of no valid velocities\n");
+ progress("\t$nBTwFlag flagged bad because of incorrect vertical velocity\n");
+
+ for (my($gi)=0; $gi<@BTw; $gi++) { # calc grid medians & mads
+ $BT{N_SAMP}[$gi] = @{$BTw[$gi]};
+ $BT{MEDIAN_W}[$gi] = median(@{$BTw[$gi]});
+ $BT{MAD_W}[$gi] = mad2($BT{W}[$gi],@{$BTw[$gi]});
+ }
+
+ &antsAddParams('BT_rms_seafloor_u',round(rms(@BTbtmu),0.01),
+ 'BT_rms_seafloor_v',round(rms(@BTbtmv),0.01),
+ 'BT_rms_seafloor_w',round(rms(@BTbtmw),0.01),
+ 'BT_avg_seafloor_u',round(avg(@BTbtmu),0.01),
+ 'BT_avg_seafloor_v',round(avg(@BTbtmv),0.01),
+ 'BT_avg_seafloor_w',round(avg(@BTbtmw),0.01),
+ 'BT_rms_tilt',round(rms(@BTtilt),0.1));
+}
+
+1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/edit_data.pl Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,169 @@
+#======================================================================
+# E D I T _ D A T A . P L
+# doc: Sat May 22 21:35:55 2010
+# dlm: Thu Dec 30 20:24:08 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 127 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# May 22, 2010: - created
+# May 24, 2010: - added editSideLobesFromSeabed()
+# Oct 29, 2010: - added editCorr_Earthcoords
+# Dec 20, 2010: - BUG: DISTANCE_TO_BIN1_CENTER & BIN_LENGTH had been
+# interpreted as along-beam, rather than vertical
+# - replaced editPitchRoll by editTilt
+# Dec 25, 2010: - adapted to changes in [LADCP_w]
+
+# NOTES:
+# - all bins must be edited (not just valid ones), to allow
+# reflr calculations to use invalid bins
+
+#======================================================================
+# $vv = countValidVels($ens)
+#======================================================================
+
+sub countValidBeamVels($)
+{
+ my($ens) = @_;
+
+ my($vv) = 0;
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ $vv += defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][0]);
+ $vv += defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][1]);
+ $vv += defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][2]);
+ $vv += defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][3]);
+ }
+ return $vv;
+}
+
+#======================================================================
+# $edited = editCorr($ens,$threshold)
+#
+# NOTES:
+# - called before Earth vels have been calculated
+#======================================================================
+
+sub editCorr($$)
+{
+ my($ens,$lim) = @_;
+
+ my($nrm) = 0;
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ for (my($beam)=0; $beam<4; $beam++) {
+ next if ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] >= $lim ||
+ !defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]));
+ undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
+ $nrm++;
+ }
+ }
+ return $nrm;
+}
+
+sub editCorr_Earthcoords($$)
+{
+ my($ens,$lim) = @_;
+
+ my($nrm) = 0;
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ my($beam);
+ for ($beam=0; $beam<4; $beam++) {
+ last unless ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] >= $lim);
+ }
+ if ($beam < 4) {
+ for (my($c)=0; $c<4; $c++) {
+ next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$c]);
+ undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$c]);
+ $nrm++;
+ }
+ }
+ }
+ return $nrm;
+}
+
+#======================================================================
+# $edited = editTilt($ens,$threshold)
+#
+# NOTES:
+# - called before Earth vels have been calculated
+# - sets TILT field for each ensemble as a side-effect
+# - for consistency with editCorr() the individual velocities are counted
+#======================================================================
+
+sub editTilt($$)
+{
+ my($ens,$lim) = @_;
+
+ $LADCP{ENSEMBLE}[$ens]->{TILT} =
+ &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+
+ return 0 if ($LADCP{ENSEMBLE}[$ens]->{TILT} <= $lim);
+
+ my($nrm) = 0;
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ for (my($beam)=0; $beam<4; $beam++) {
+ next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
+ undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$beam]);
+ $nrm++;
+ }
+ }
+ return $nrm;
+}
+
+#======================================================================
+# $edited = editErrVel($ens,$threshold)
+#
+# NOTES:
+# - call after Earth vels have been calculated
+#======================================================================
+
+sub editErrVel($$)
+{
+ my($ens,$lim) = @_;
+
+ my($nrm) = 0;
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ next if (abs($LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) <= $lim);
+ undef($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
+ $nrm++
+ }
+ return $nrm;
+}
+
+#======================================================================
+# ($nvrm,$nerm) = editSideLobes($fromEns,$toEns,$range)
+#
+# NOTES:
+# 1) When this code is executed the sound speed is known. No attempt is made to correct for
+# along-beam soundspeed variation, but the soundspeed at the transducer is accounted for.
+#======================================================================
+
+sub editSideLobes($$$)
+{
+ my($fe,$te,$wd) = @_; # first & last ens to process, water depth for downlooker
+ my($nvrm) = 0; # of velocities removed
+ my($nerm) = 0; # of ensembles affected
+ for (my($e)=$fe; $e<=$te; $e++) {
+ next unless numberp($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH});
+ my($range) = $LADCP{ENSEMBLE}[$e]->{XDUCER_FACING_UP}
+ ? $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}
+ : $wd - $LADCP{ENSEMBLE}[$e]->{CTD_DEPTH};
+ my($sscorr) = $CTD{SVEL}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}] / 1500;
+ my($goodBins) = ($range - $sscorr*$LADCP{DISTANCE_TO_BIN1_CENTER})
+ / ($sscorr*$LADCP{BIN_LENGTH})
+ - 1.5;
+
+ my($dirty) = 0;
+ for (my($bin)=int($goodBins); $bin<$LADCP{N_BINS}; $bin++) { # NB: 2 good bins implies that bin 2 is bad
+ next unless ($bin>=0 && defined($LADCP{ENSEMBLE}[$e]->{W}[$bin]));
+ $dirty = 1;
+ $nvrm++;
+ undef($LADCP{ENSEMBLE}[$e]->{W}[$bin]);
+ }
+
+ $nerm += $dirty;
+ }
+ return ($nvrm,$nerm);
+}
+
+1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/find_seabed.pl Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,98 @@
+#======================================================================
+# F I N D _ S E A B E D . P L
+# doc: Sun May 23 20:26:11 2010
+# dlm: Thu Dec 30 21:45:02 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 12 0 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# May 23, 2010: - adapted from [perl-tools/RDI_Utils.pl]
+# Dec 25, 2010: - adapted to changes in [LADCP_w]
+
+# NOTES:
+# 1) BT range is corrected for sound speed at the transducer. This is not
+# accurate, but unlikely to be very wrong, at least for deep casts,
+# because the vertical sound speed variability near the seabed tends
+# to be small. The seabed depth is only used for sidelobe editing,
+# which is done with a generous safety margin (from the UH shear
+# method implementation).
+# 2) Acquisition sound speed of 1500m/s assumed.
+# 3) To be reasonably accurate, DEPTH must be from the CTD at this stage.
+
+#======================================================================
+# (seabed median depth, mad) = find_seabed(dta ptr, btm ensno, coord flag)
+#======================================================================
+
+my($search_width) = 200; # # of ensembles around bottom to search
+my($mode_width) = 10; # max range of bottom around mode
+my($min_dist) = 20; # min dist to seabed for good data
+my($z_offset) = 10000; # shift z to ensure +ve array indices
+
+sub find_seabed($$$)
+{
+ my($d,$be,$beamCoords) = @_;
+ my($i,$dd,$sd,$nd);
+ my(@guesses);
+
+ return undef unless ($be-$search_width >= 0 &&
+ $be+$search_width <= $#{$d->{ENSEMBLE}});
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless (defined($d->{ENSEMBLE}[$i]->{CTD_DEPTH}) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[0]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[1]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[2]) &&
+ defined($d->{ENSEMBLE}[$i]->{BT_RANGE}[3]));
+ my(@BT) = $beamCoords
+ ? velInstrumentToEarth($d,$i,
+ velBeamToInstrument($d,
+ @{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}}))
+ : velApplyHdgBias($d,$i,@{$d->{ENSEMBLE}[$i]->{BT_VELOCITY}});
+ next unless (abs($BT[3]) < 0.05);
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} =
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[0]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[1]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[2]/4 +
+ $d->{ENSEMBLE}[$i]->{BT_RANGE}[3]/4;
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} *= cos(rad($d->{BEAM_ANGLE}));
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} *= $CTD{SVEL}[$d->{ENSEMBLE}[$i]->{CTD_SCAN}]/1500;
+ next unless ($d->{ENSEMBLE}[$i]->{DEPTH_BT} >= $min_dist);
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} *= -1
+ if ($d->{ENSEMBLE}[$i]->{XDUCER_FACING_UP});
+ $d->{ENSEMBLE}[$i]->{DEPTH_BT} += $d->{ENSEMBLE}[$i]->{CTD_DEPTH};
+ if ($d->{ENSEMBLE}[$i]->{DEPTH_BT} > $d->{ENSEMBLE}[$be]->{CTD_DEPTH}) {
+ $guesses[int($d->{ENSEMBLE}[$i]->{DEPTH_BT})+$z_offset]++;
+ $nd++;
+ } else {
+ undef($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ }
+ }
+ return undef unless ($nd>5);
+
+ my($mode,$nmax);
+ for ($i=0; $i<=$#guesses; $i++) { # find mode
+ $nmax=$guesses[$i],$mode=$i-$z_offset
+ if ($guesses[$i] > $nmax);
+ }
+
+ $nd = 0;
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ if (abs($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$mode) <= $mode_width) {
+ $dd += $d->{ENSEMBLE}[$i]->{DEPTH_BT};
+ $nd++;
+ } else {
+ undef($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ }
+ }
+ return undef unless ($nd >= 2);
+
+ $dd /= $nd;
+ for ($i=$be-$search_width; $i<=$be+$search_width; $i++) {
+ next unless defined($d->{ENSEMBLE}[$i]->{DEPTH_BT});
+ $sd += ($d->{ENSEMBLE}[$i]->{DEPTH_BT}-$dd)**2;
+ }
+
+ return ($dd, sqrt($sd/($nd-1)));
+}
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/svel_corrections.pl Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,89 @@
+#======================================================================
+# S V E L _ C O R R E C T I O N S . P L
+# doc: Thu Dec 30 01:35:18 2010
+# dlm: Thu Dec 30 01:40:13 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 86 20 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+sub sscorr_w($$$$) # sound-speed correction for w
+{ # see RDI Coord. Trans. manual sec. 4.1, ...
+ my($wObs,$wCTD,$dADCP,$dBin) = @_; # but there is an error: the ^2 applies to the []
+ my($tanSqBeamAngle) = tan(rad($LADCP{BEAM_ANGLE}))**2;
+
+ $dADCP = int($dADCP); # @sVelProf is binned to 1m
+ $dBin = int($dBin);
+
+ while (!numberp($sVelProf[$dADCP])) { $dADCP--; } # skip gaps & bottom of profile
+ while (!numberp($sVelProf[$dBin ])) { $dBin--; }
+
+ my($Kn) = sqrt(1 + (1 - $sVelProf[$dBin]/$sVelProf[$dADCP])**2 * $tanSqBeamAngle);
+ return ($wObs*$sVelProf[$dBin]/1500 - $wCTD) / $Kn;
+}
+
+sub calc_binDepths($) # see RDI Coord Trans manual sec. 4.2
+{
+ my($ens) = @_;
+ my(@bindz);
+
+ my($tanSqBeamAngle) = tan(rad($LADCP{BEAM_ANGLE}))**2;
+ my($curdz) = 0; # calc avg sndspeed btw transducer & 1st bin
+ $curdz-- until numberp($sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)]);
+ my($avgss) = my($ADCPss) = $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)];
+
+ my($sumss) = my($nss) = 0;
+ if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) {
+ while ($curdz >= -$LADCP{DISTANCE_TO_BIN1_CENTER}*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))) {
+ if (numberp($sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)])) {
+ $sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)]; $nss++;
+ }
+ $curdz--;
+ }
+ } else {
+ while ($curdz <= $LADCP{DISTANCE_TO_BIN1_CENTER}*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))) {
+ if (numberp($sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)])) {
+ $sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)]; $nss++;
+ }
+ $curdz++;
+ }
+ }
+ $avgss = $sumss/$nss if ($nss>0);
+
+ my($Kn) = sqrt(1 + (1 - $avgss/$ADCPss)**2 * $tanSqBeamAngle);
+ $bindz[0] = $LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} ?
+ - $LADCP{DISTANCE_TO_BIN1_CENTER}*$Kn*$avgss/1500*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT})) :
+ + $LADCP{DISTANCE_TO_BIN1_CENTER}*$Kn*$avgss/1500*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}));
+
+ for (my($bin)=1; $bin<=$LADCP_lastBin-1; $bin++) {
+ $sumss = $nss = 0;
+ if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) {
+ while ($curdz >= $bindz[$bin-1]-$LADCP{BIN_LENGTH}*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))) {
+ if (numberp($sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)])) {
+ $sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)]; $nss++;
+ }
+ $curdz--;
+ }
+ } else {
+ while ($curdz <= $bindz[$bin-1]+$LADCP{BIN_LENGTH}*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}))) {
+ if (numberp($sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)])) {
+ $sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}+$curdz)]; $nss++;
+ }
+ $curdz++;
+ }
+ }
+ $avgss = $sumss/$nss if ($nss > 0); # otherwise, leave avgss as is
+
+ $Kn = sqrt(1 + (1 - $avgss/$ADCPss)**2 * $tanSqBeamAngle);
+ $bindz[$bin] = $LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} ?
+ $bindz[$bin-1] - $LADCP{BIN_LENGTH}*$Kn*$avgss/1500*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT})) :
+ $bindz[$bin-1] + $LADCP{BIN_LENGTH}*$Kn*$avgss/1500*cos(rad($LADCP{ENSEMBLE}[$ens]->{TILT}));
+ }
+
+ my(@bindepth);
+ for (my($i)=0; $i<@bindz; $i++) {
+ $bindepth[$i] = $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} + $bindz[$i];
+ }
+ return @bindepth;
+}
+
+1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/time_lag.pl Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,159 @@
+#======================================================================
+# T I M E _ L A G . P L
+# doc: Fri Dec 17 21:59:07 2010
+# dlm: Sun Jun 26 20:47:09 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 17 80 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Dec 17, 2010: - created
+# Dec 18, 2010: - adapted for multi-pass lagging
+# Dec 20: 2010: - added code to adjust start and end of profile ens
+# based on extent of CTD profile and guestimated time
+# ofset
+# Jun 26, 2010: - added heuristic to chose between weighted-mean and
+# unambiguously best offsets
+# - turned -3 criterion into warning when 3 lags are consecutive
+
+# TODO:
+# - better seabed code (from LADCPproc)
+# - intermediate-step timelagging guess
+# - flip aliased ensembles
+
+sub mad_w($$$) # mean absolute deviation
+{
+ my($fe,$le,$so) = @_; # first/last LADCP ens, CTD scan offset
+ my($sad) = my($n) = 0;
+
+ for (my($e)=$fe; $e<=$le; $e++) {
+ my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
+ die("assertion failed\n" .
+ "\ttest: abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2\n" .
+ "\te = $e, s = $s, ensemble = $LADCP{ENSEMBLE}[$e]->{NUMBER}"
+ ) unless (abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2);
+
+ my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W} - $CTD{W}[$s+$so];
+ next unless (abs($dw) <= $opt_m);
+
+ $sad += abs($dw);
+ $n++;
+ }
+ return ($n>0) ? $sad/$n : 9e99; # n == 0, e.g. in bottom gap
+}
+
+
+sub bestLag($$$$) # find best lag in window
+{
+ my($fe,$le,$ww,$soi) = @_; # first/last LADCP ens, window width, scan-offset increment
+
+ my($bestso) = 0; # error at first-guess offset
+ my($bestmad) = mad_w($fe,$le,0);
+
+ for (my($dso) = 1; $dso <= int($ww/2/$CTD{DT} + 0.5); $dso+=$soi) {
+ my($mad) = mad_w($fe,$le,-$dso);
+ $bestmad=$mad,$bestso=-$dso if ($mad < $bestmad);
+ $mad = mad_w($fe,$le,$dso);
+ $bestmad=$mad,$bestso=$dso if ($mad < $bestmad);
+ }
+ return ($bestso,$bestmad);
+}
+
+#----------------------------------------------------------------------
+# carry out lag correlations and keep tally of the results
+# - fist and last 10% of LADCP profile ignored
+#----------------------------------------------------------------------
+
+sub calc_lag($$$)
+{
+ my($n_windows,$w_size,$scan_increment) = @_;
+
+RETRY:
+ progress("Calculating $n_windows time-lags from ${w_size}s-long windows at %dHz resolution...\n",
+ int(1/$scan_increment/$CTD{DT}+0.5));
+
+ my($approx_CTD_profile_start_ens) =
+ $firstGoodEns + int(($CTD{ELAPSED}[0] - $CTD{TIME_LAG}) / $LADCP{MEAN_DT});
+ my($approx_CTD_profile_end_ens) =
+ $firstGoodEns + int(($CTD{ELAPSED}[$#{$CTD{ELAPSED}}] + $CTD{ELAPSED}[0] - $CTD{TIME_LAG}) / $LADCP{MEAN_DT});
+
+ my($approx_joint_profile_start_ens) = max($firstGoodEns,$approx_CTD_profile_start_ens);
+ my($approx_joint_profile_end_ens) = min($lastGoodEns,$approx_CTD_profile_end_ens);
+ debugmsg("profile start: $firstGoodEns -> $approx_joint_profile_start_ens\n");
+ debugmsg("profile end : $lastGoodEns -> $approx_joint_profile_end_ens\n");
+
+ my($skip_ens) = int(($approx_joint_profile_end_ens - $approx_joint_profile_start_ens) / 10 + 0.5);
+
+ my(%nBest);
+ for (my($wi)=0; $wi<$n_windows; $wi++) {
+ my($fe) = $approx_joint_profile_start_ens + $skip_ens + $wi*int(($approx_joint_profile_end_ens-$approx_joint_profile_start_ens-2*$skip_ens)/$n_windows+0.5);
+ my($so,$mad) = bestLag($fe,$fe+int($w_size/$LADCP{MEAN_DT}+0.5),$w_size,$scan_increment);
+ debugmsg("%.1f cm/s mad(w) at %3d scans offset\n",100*$mad,$so);
+ $nBest{$so}++;
+ }
+
+ my(@best_lag);
+ foreach my $i (keys(%nBest)) {
+ $best_lag[0] = $i if ($nBest{$i} > $nBest{$best_lag[0]});
+ }
+ foreach my $i (keys(%nBest)) {
+ next if ($i == $best_lag[0]);
+ $best_lag[1] = $i if ($nBest{$i} > $nBest{$best_lag[1]});
+ }
+ foreach my $i (keys(%nBest)) {
+ next if ($i == $best_lag[0] || $i == $best_lag[1]);
+ $best_lag[2] = $i if ($nBest{$i} > $nBest{$best_lag[2]});
+ }
+ progress("\t3 most popular offsets: %d (%d%%), %d (%d%%), %d (%d%%)\n",
+ $best_lag[0],int(($nBest{$best_lag[0]}/$n_windows)*100+0.5),
+ $best_lag[1],int(($nBest{$best_lag[1]}/$n_windows)*100+0.5),
+ $best_lag[2],int(($nBest{$best_lag[2]}/$n_windows)*100+0.5));
+
+ if ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} <= 6) {
+ warning(0,"cannot determine valid lag => trying again with doubled window size\n");
+ undef(%nBest);
+ $w_size *= 2;
+ goto RETRY;
+ }
+
+ unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $opt_3*$n_windows) {
+ if (max(@best_lag)-min(@best_lag) > 2) {
+ croak(sprintf("$0: cannot determine a valid lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
+ int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5)))
+ } else {
+ warning(1,"top 3 tags account for only %d%% of total\n",
+ int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5));
+ }
+ }
+
+ my($bmo);
+ if (max(@best_lag)-min(@best_lag) > 5 || $nBest{$best_lag[0]}/$n_windows >= 0.5) {
+ $bmo = $best_lag[0];
+ progress("\tunambiguously best offset = %d scans\n",$bmo);
+ } else {
+ $bmo = ($nBest{$best_lag[0]}*$best_lag[0] +
+ $nBest{$best_lag[1]}*$best_lag[1] +
+ $nBest{$best_lag[2]}*$best_lag[2]) / ($nBest{$best_lag[0]} +
+ $nBest{$best_lag[1]} +
+ $nBest{$best_lag[2]});
+ progress("\tweighted-mean offset = %.1f scans\n",$bmo);
+ }
+
+ if ($bmo > 0.9*$w_size/2/$CTD{DT}) {
+ warning(0,"lag too close to the edge of the window --- trying again after adjusting the guestimated offset\n");
+ $CTD{TIME_LAG} += $w_size/2;
+ undef(%nBest);
+ goto RETRY;
+ }
+ if (-$bmo > 0.9*$w_size/2/$CTD{DT}) {
+ warning(0,"lag too close to the edge of the window --- trying again after adjusting the guestimated offset\n");
+ $CTD{TIME_LAG} -= $w_size/2;
+ undef(%nBest);
+ goto RETRY;
+ }
+
+ return $CTD{TIME_LAG}+$bmo*$CTD{DT};
+}
+
+
+1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Jul_04_2011/time_series.pl Sat Jul 09 15:14:33 2011 -0400
@@ -0,0 +1,116 @@
+#======================================================================
+# T I M E _ S E R I E S . P L
+# doc: Sun May 23 16:40:53 2010
+# dlm: Sat Jul 2 21:12:35 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 15 48 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# May 23, 2010: - created from [perl-tools/RDI_Utils.pl]
+# Oct 20, 2010: - disabled max_gap profile-restarting code
+# Dec 17, 2010: - re-added {DEPTH} field to ensembles
+# Dec 18, 2010: - max gap re-enabled
+# Dec 20, 2010: - cosmetics
+# Jul 2, 2011: - tightened gap-detection code
+
+# NOTES:
+# - resulting DEPTH field based on integrated w without any sound speed correction
+# - single-ping ensembles assumed, i.e. no percent-good tests applied
+# - specified bin numbers are 1-relative
+
+sub ref_lr_w($$$$) # calc ref-layer vert vels
+{
+ my($dta,$ens,$rl_b0,$rl_b1) = @_;
+ my($i,@n,@bn,@v,@vel,@bv,@w);
+
+ for ($i=$rl_b0-1; $i<=$rl_b1-1; $i++) {
+ if (defined($dta->{ENSEMBLE}[$ens]->{W}[$i])) { # valid w
+ $vel[2] += $dta->{ENSEMBLE}[$ens]->{W}[$i]; $n[2]++;
+ $vel[3] += $dta->{ENSEMBLE}[$ens]->{ERRVEL}[$i], $n[3]++ if defined($dta->{ENSEMBLE}[$ens]->{ERRVEL}[$i]);
+ push(@w,$dta->{ENSEMBLE}[$ens]->{W}[$i]); # for stderr test
+ }
+ }
+
+ my($w) = $n[2] ? $vel[2]/$n[2] : undef; # w uncertainty
+ my($sumsq) = 0;
+ for ($i=0; $i<=$#w; $i++) {
+ $sumsq += ($w-$w[$i])**2;
+ }
+ my($stderr) = $n[2]>=2 ? sqrt($sumsq)/($n[2]-1) : undef;
+
+ if (defined($w)) { # valid w
+ $dta->{ENSEMBLE}[$ens]->{REFLR_W} = $w;
+ $dta->{ENSEMBLE}[$ens]->{REFLR_W_ERR} = $stderr;
+ }
+}
+
+#======================================================================
+# ($firstgood,$lastgood,$atbottom,$w_gap_time) =
+# calcLADCPts($dta,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap);
+#======================================================================
+
+sub calcLADCPts($$$$)
+{
+ my($dta,$rl_b0,$rl_b1,$max_gap) = @_;
+ my($firstgood,$lastgood,$atbottom,$w_gap_time,$max_depth);
+
+ for (my($depth)=0,my($e)=0; $e<=$#{$dta->{ENSEMBLE}}; $e++) {
+
+ ref_lr_w($dta,$e,$rl_b0,$rl_b1);
+
+ if (defined($firstgood)) {
+ $dta->{ENSEMBLE}[$e]->{ELAPSED} = # time since start
+ $dta->{ENSEMBLE}[$e]->{UNIX_TIME} -
+ $dta->{ENSEMBLE}[$firstgood]->{UNIX_TIME};
+ } else {
+ if (defined($dta->{ENSEMBLE}[$e]->{REFLR_W})) { # start of prof.
+ $firstgood = $lastgood = $e;
+ $dta->{ENSEMBLE}[$e]->{ELAPSED} = 0;
+ }
+ next;
+ }
+
+ unless (defined($dta->{ENSEMBLE}[$e]->{REFLR_W})) { # gap
+ $w_gap_time += $dta->{ENSEMBLE}[$e]->{UNIX_TIME} -
+ $dta->{ENSEMBLE}[$e-1]->{UNIX_TIME};
+ next;
+ }
+
+ my($dt) = $dta->{ENSEMBLE}[$e]->{UNIX_TIME} - # time step since
+ $dta->{ENSEMBLE}[$lastgood]->{UNIX_TIME}; # ... last good ens
+
+ if ($dt > $max_gap) {
+ if ($max_depth>50 && $depth<0.1*$max_depth) {
+ warning(1,"long gap (%ds) near end of profile --- terminated at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
+ last;
+ }
+ if ($depth < 10) {
+ warning(1,"long gap (%ds) near beginning of profile --- restarted at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
+ $firstgood = $lastgood = $e;
+ undef($atbottom); undef($max_depth);
+ $depth = 0;
+ $dta->{ENSEMBLE}[$e]->{ELAPSED} = 0;
+ $w_gap_time = 0;
+ next;
+ }
+ if ($dta->{ENSEMBLE}[$e]->{ELAPSED} < 200) {
+ warning(1,"long gap (%ds) at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}, %ds into the profile\n",
+ $dt,$dta->{ENSEMBLE}[$e]->{ELAPSED});
+ } else {
+ warning(1,"long gap (%ds) at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}, %.1fmin into the profile\n",
+ $dt,$dta->{ENSEMBLE}[$e]->{ELAPSED}/60);
+ }
+ }
+
+ $depth += $dta->{ENSEMBLE}[$lastgood]->{REFLR_W} * $dt; # integrate
+ $dta->{ENSEMBLE}[$e]->{DEPTH} = $depth;
+
+ $atbottom = $e, $max_depth = $depth if ($depth > $max_depth);
+ $lastgood = $e;
+ }
+
+ return ($firstgood,$lastgood,$atbottom,$w_gap_time);
+}
+
+1;
--- 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();
--- a/time_lag.pl Fri Jun 17 13:33:47 2011 -0400
+++ b/time_lag.pl Sat Jul 09 15:14:33 2011 -0400
@@ -1,9 +1,9 @@
#======================================================================
-# C A L C _ T I M E L A G S . P L
+# T I M E _ L A G . P L
# doc: Fri Dec 17 21:59:07 2010
-# dlm: Tue Dec 21 00:40:33 2010
+# dlm: Fri Jul 8 04:11:52 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 117 20 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 60 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -12,27 +12,52 @@
# Dec 20: 2010: - added code to adjust start and end of profile ens
# based on extent of CTD profile and guestimated time
# ofset
+# Jun 26, 2010: - added heuristic to chose between weighted-mean and
+# unambiguously best offsets
+# - turned -3 criterion into warning when 3 lags are consecutive
+# Jul 4, 2011: - increased MAX_ALLOWED_THREE_LAG_SPREAD from 2 to 3
+# Jul 7, 2011: - removed window-mean w before time lagging to allow lagging
+# of casts with large w
# TODO:
# - better seabed code (from LADCPproc)
# - intermediate-step timelagging guess
# - flip aliased ensembles
+my($MAX_ALLOWED_THREE_LAG_SPREAD) = 3; # this was initially set to 2 but found to be
+ # violated quite often during 2011_IWISE. A
+ # large spread may indicate dropped CTD scans.
+ # The optimum value may be cast-duration dependent.
+
sub mad_w($$$) # mean absolute deviation
{
my($fe,$le,$so) = @_; # first/last LADCP ens, CTD scan offset
my($sad) = my($n) = 0;
- for (my($e)=$fe; $e<=$le; $e++) {
+ my($LADCP_mean_w,$CTD_mean_w,$nsamp) = (0,0,0);
+ for (my($e)=$fe; $e<=$le; $e++) { # first, calculate mean w in window
my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
die("assertion failed\n" .
"\ttest: abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2\n" .
"\te = $e, s = $s, ensemble = $LADCP{ENSEMBLE}[$e]->{NUMBER}"
) unless (abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2);
+ next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W});
+ my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w);
+ next unless (abs($dw) <= $opt_m);
- my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W} - $CTD{W}[$s+$so];
+ $LADCP_mean_w += $LADCP{ENSEMBLE}[$e]->{REFLR_W};
+ $CTD_mean_w += $CTD{W}[$s+$so];
+ $nsamp++;
+ }
+ return 9e99 unless ($nsamp);
+ $LADCP_mean_w /= $nsamp;
+ $CTD_mean_w /= $nsamp;
+
+ for (my($e)=$fe; $e<=$le; $e++) { # now, calculate mad
+ my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
+ my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w);
+ next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W});
next unless (abs($dw) <= $opt_m);
-
$sad += abs($dw);
$n++;
}
@@ -113,31 +138,43 @@
goto RETRY;
}
- croak(sprintf("$0: cannot determine a valid lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
- int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5)))
- unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $opt_3*$n_windows);
+ unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $opt_3*$n_windows) {
+ if (max(@best_lag)-min(@best_lag) > $MAX_ALLOWED_THREE_LAG_SPREAD) {
+ croak(sprintf("$0: cannot determine a valid lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
+ int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5)))
+ } else {
+ warning(1,"top 3 tags account for only %d%% of total\n",
+ int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5));
+ }
+ }
- my($wmo) = ($nBest{$best_lag[0]}*$best_lag[0] +
+ my($bmo);
+ if (max(@best_lag)-min(@best_lag) > 5 || $nBest{$best_lag[0]}/$n_windows >= 0.5) {
+ $bmo = $best_lag[0];
+ progress("\tunambiguously best offset = %d scans\n",$bmo);
+ } else {
+ $bmo = ($nBest{$best_lag[0]}*$best_lag[0] +
$nBest{$best_lag[1]}*$best_lag[1] +
$nBest{$best_lag[2]}*$best_lag[2]) / ($nBest{$best_lag[0]} +
$nBest{$best_lag[1]} +
$nBest{$best_lag[2]});
- progress("\tweighted-mean offset = %.1f scans\n",$wmo);
+ progress("\tweighted-mean offset = %.1f scans\n",$bmo);
+ }
- if ($wmo > 0.9*$w_size/2/$CTD{DT}) {
+ if ($bmo > 0.9*$w_size/2/$CTD{DT}) {
warning(0,"lag too close to the edge of the window --- trying again after adjusting the guestimated offset\n");
$CTD{TIME_LAG} += $w_size/2;
undef(%nBest);
goto RETRY;
}
- if (-$wmo > 0.9*$w_size/2/$CTD{DT}) {
+ if (-$bmo > 0.9*$w_size/2/$CTD{DT}) {
warning(0,"lag too close to the edge of the window --- trying again after adjusting the guestimated offset\n");
$CTD{TIME_LAG} -= $w_size/2;
undef(%nBest);
goto RETRY;
}
- return $CTD{TIME_LAG}+$wmo*$CTD{DT};
+ return $CTD{TIME_LAG}+$bmo*$CTD{DT};
}
--- a/time_series.pl Fri Jun 17 13:33:47 2011 -0400
+++ b/time_series.pl Sat Jul 09 15:14:33 2011 -0400
@@ -1,9 +1,9 @@
#======================================================================
-# C A L C _ L A D C P _ T I S . P L
+# T I M E _ S E R I E S . P L
# doc: Sun May 23 16:40:53 2010
-# dlm: Tue Dec 21 00:34:15 2010
+# dlm: Mon Jul 4 01:53:04 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 87 22 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 16 47 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -12,6 +12,8 @@
# Dec 17, 2010: - re-added {DEPTH} field to ensembles
# Dec 18, 2010: - max gap re-enabled
# Dec 20, 2010: - cosmetics
+# Jul 2, 2011: - tightened gap-detection code
+# Jul 4, 2011: - added support for $skip_ens
# NOTES:
# - resulting DEPTH field based on integrated w without any sound speed correction
@@ -46,15 +48,15 @@
#======================================================================
# ($firstgood,$lastgood,$atbottom,$w_gap_time) =
-# calcLADCPts($dta,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap);
+# calcLADCPts($dta,$skip_ens,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap);
#======================================================================
sub calcLADCPts($$$$)
{
- my($dta,$rl_b0,$rl_b1,$max_gap) = @_;
+ my($dta,$skip_ens,$rl_b0,$rl_b1,$max_gap) = @_;
my($firstgood,$lastgood,$atbottom,$w_gap_time,$max_depth);
- for (my($depth)=0,my($e)=0; $e<=$#{$dta->{ENSEMBLE}}; $e++) {
+ for (my($depth)=0,my($e)=$skip_ens; $e<=$#{$dta->{ENSEMBLE}}; $e++) {
ref_lr_w($dta,$e,$rl_b0,$rl_b1);
@@ -80,17 +82,26 @@
$dta->{ENSEMBLE}[$lastgood]->{UNIX_TIME}; # ... last good ens
if ($dt > $max_gap) {
- if ($max_depth>50 && $depth<0.5*$max_depth) {
- warning(1,"long gap (%ds). Profile ended at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
+ if ($max_depth>50 && $depth<0.1*$max_depth) {
+ warning(1,"long gap (%ds) near end of profile --- terminated at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
last;
}
- warning(1,"long gap (%ds). Profile restarted at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
- $firstgood = $lastgood = $e;
- undef($atbottom); undef($max_depth);
- $depth = 0;
- $dta->{ENSEMBLE}[$e]->{ELAPSED} = 0;
- $w_gap_time = 0;
- next;
+ if ($depth < 10) {
+ warning(1,"long gap (%ds) near beginning of profile --- restarted at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}\n",$dt);
+ $firstgood = $lastgood = $e;
+ undef($atbottom); undef($max_depth);
+ $depth = 0;
+ $dta->{ENSEMBLE}[$e]->{ELAPSED} = 0;
+ $w_gap_time = 0;
+ next;
+ }
+ if ($dta->{ENSEMBLE}[$e]->{ELAPSED} < 200) {
+ warning(1,"long gap (%ds) at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}, %ds into the profile\n",
+ $dt,$dta->{ENSEMBLE}[$e]->{ELAPSED});
+ } else {
+ warning(1,"long gap (%ds) at ensemble #$dta->{ENSEMBLE}[$e]->{NUMBER}, %.1fmin into the profile\n",
+ $dt,$dta->{ENSEMBLE}[$e]->{ELAPSED}/60);
+ }
}
$depth += $dta->{ENSEMBLE}[$lastgood]->{REFLR_W} * $dt; # integrate