#!/usr/bin/perl
#======================================================================
# L A D C P _ W
# doc: Fri Dec 17 18:11:13 2010
# dlm: Sat Dec 25 03:03:35 2010
# (c) 2010 A.M. Thurnherr
# uE-Info: 11 91 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
$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
# CTD REQUIREMENTS
# - elapsed elapsed seconds; see note below
# - depth
# - ss sound speed
# - w ddepth/dt
# - winch_w low-frequency w
# - wave_w high-frequency w
# - w_t dw/dt --- OPTIONAL and not used during processing
# 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_LADCP_data.pl";
require "$WCALC/calc_LADCP_tis.pl";
require "$WCALC/calc_timelags.pl";
require "$WCALC/find_seabed.pl";
require "$PERL_TOOLS/RDI_BB_Read.pl";
require "$PERL_TOOLS/RDI_Coords.pl";
$antsParseHeader = 0;
&antsUsage('3:4a:b:c:d:e:f:g:h:i:m:n:o:p:r:s:t:v:w:x:',1,
'[-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>]',
'[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_correction_params/uc_correction_param>]',
'[use LADCP -b)ins for w profiles <bin[2],bin[*]>',
'[-o)utput bin <resolution[10m]>] [require -s)amples <min[20]>]',
'[-f) write time-series <file>] [-d)ump depth-bins to <basename>] [-p)ackage effect <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,5); # 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
$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
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_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();
#----------------------------------------------------------------------
# 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);
}
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));
} 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);
}
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));
}
#-------------------------------------------------------------------
# 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
#-------------------------------------------------------------------
if ($LADCP{BEAM_COORDINATES}) {
progress("Calculating earth-coordinate velocities...\n");
$nvw = 0;
for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
($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});
}
}
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)=0; $bin<$LADCP{N_BINS}; $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]);
}
}
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());
($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w,$CTD_winch_w,$CTD_wave_w) =
&fnr('elapsed','depth','ss','w','winch_w','wave_w');
$CTD_w_t = &fnrNoErr('w_t');
$CTD_maxdepth = -1;
do {
croak("$0: cannot deal with non-numeric CTD fields\n")
unless &antsNumbers($CTD_elapsed,$CTD_depth,$CTD_w);
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{WINCH_W}},$ants_[0][$CTD_winch_w]);
push(@{$CTD{WAVE_W}}, $ants_[0][$CTD_wave_w]);
push(@{$CTD{W_T}}, defined($CTD_W_T) ? $ants_[0][$CTD_W_T] : nan);
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]);
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);
$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});
}
$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($warned),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
info("CTD data begin with LADCP in water\n") unless ($warned++);
$firstGoodEns = $ens+1;
next;
}
unless ($scan <= $#{$CTD{ELAPSED}}) {
info("CTD data end with LADCP in water\n");
$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
$LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $CTD{SVEL}[$scan]/1500; # correct for sound-speed variations
$LADCP{ENSEMBLE}[$ens]->{ELAPSED_MISMATCH} =
$CTD{ELAPSED}[$scan] - ($LADCP{ENSEMBLE}[$ens]->{ELAPSED} + $CTD{TIME_LAG});
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]->{DEPTH} = $CTD{DEPTH}[$scan];
$LADCP{ENSEMBLE}[$ens]->{SOUND_SPEED} = $CTD{SVEL}[$scan];
$LADCP{ENSEMBLE}[$ens]->{CTD_W} = $CTD{W}[$scan];
$LADCP{ENSEMBLE}[$ens]->{CTD_W_T} = $CTD{W_T}[$scan];
$LADCP{ENSEMBLE}[$ens]->{WINCH_W} = $CTD{WINCH_W}[$scan];
$LADCP{ENSEMBLE}[$ens]->{WAVE_W} = $CTD{WAVE_W}[$scan];
my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $LADCP{ENSEMBLE}[$ens]->{CTD_W};
if (abs($reflr_ocean_w) <= $opt_m) {
$sumWsq += &SQR($reflr_ocean_w);
$nWsq++;
if ($LADCP{ENSEMBLE}[$ens]->{DEPTH} > 100 &&
$LADCP{ENSEMBLE}[$ens]->{DEPTH} < $LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH}-100) {
$sumWsqI += &SQR($reflr_ocean_w);
$nWsqI++;
}
} else {
undef($LADCP{ENSEMBLE}[$ens]->{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]->{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");
}
#-------------------
# Find or set seabed
#--------------------
if (numberp($opt_h)) {
progress("Setting water depth (-h)\n");
$water_depth = $opt_h;
$sig_water_depth = 0;
} else {
progress("Finding seabed...\n");
($water_depth,$sig_water_depth) =
find_seabed(\%LADCP,$LADCP_atbottom,$LADCP{BEAM_COORDINATES});
}
&antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth);
progress("\t%d(%d) m water depth\n",$water_depth,$sig_water_depth)
if defined($water_depth);
#----------------------------------------------------------------------------
# Remove data contaminated by sidelobe reflection from seabed and sea surface
#----------------------------------------------------------------------------
if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {
&antsAddParams('ADCP_orientation','downlooker');
if (defined($water_depth)) {
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("could not find seabed --- 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,$water_depth);
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
#---------------------------------------------------------------------------
my($tanSqBeamAngle) = tan(rad($LADCP{BEAM_ANGLE}))**2;
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 []
$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 sscorr_bindz($) # see RDI Coord Trans manual sec. 4.2
{
my($ens) = @_;
my(@bindz);
my($curdz) = 0; # calc avg sndspeed btw transducer & 1st bin
$curdz-- until numberp($sVelProf[int($LADCP{ENSEMBLE}[$ens]->{DEPTH}+$curdz)]);
my($avgss) = my($ADCPss) = $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{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]->{DEPTH}+$curdz)])) {
$sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{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]->{DEPTH}+$curdz)])) {
$sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{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]->{DEPTH}+$curdz)])) {
$sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{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]->{DEPTH}+$curdz)])) {
$sumss += $sVelProf[int($LADCP{ENSEMBLE}[$ens]->{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}));
}
return @bindz;
}
progress("Binning velocities...\n");
progress("\tdowncast...\n");
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
next unless numberp($LADCP{ENSEMBLE}[$ens]->{DEPTH});
my(@dz) = sscorr_bindz($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($depth) = $LADCP{ENSEMBLE}[$ens]->{DEPTH} + $dz[$bin]; # depth-bin velocities
$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{CTD_W},
$LADCP{ENSEMBLE}[$ens]->{DEPTH},
$depth);
my($bi) = $depth/$opt_o;
push(@{$DNCAST{ENS}[$bi]},$ens);
push(@{$DNCAST{ELAPSED}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED});
push(@{$DNCAST{CTD_W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CTD_W});
push(@{$DNCAST{WINCH_W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{WINCH_W});
push(@{$DNCAST{WAVE_W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{WAVE_W});
push(@{$DNCAST{BIN}[$bi]},$bin);
push(@{$DNCAST{DEPTH}[$bi]},$depth);
push(@{$DNCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin]);
push(@{$DNCAST{PKGCORR_W}[10*round($LADCP{ENSEMBLE}[$ens]->{WAVE_W},0.1)+50]},
$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin])
if defined($opt_p) && abs($LADCP{ENSEMBLE}[$ens]->{WAVE_W}) < 5;
}
}
if (defined($opt_x)) {
progress("\t\tapplying wave_w correction...\n");
for (my($bi)=0; $bi<=$#{$DNCAST{ENS}}; $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]*$DNCAST{WAVE_W}[$bi][$i]**$e;
}
}
}
}
for (my($bi)=0; $bi<=$#{$DNCAST{ENS}}; $bi++) {
$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]->{DEPTH});
my(@dz) = sscorr_bindz($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
my($depth) = $LADCP{ENSEMBLE}[$ens]->{DEPTH} + $dz[$bin];
$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin] =
sscorr_w($LADCP{ENSEMBLE}[$ens]->{W}[$bin],
$LADCP{ENSEMBLE}[$ens]->{CTD_W},
$LADCP{ENSEMBLE}[$ens]->{DEPTH},
$depth);
my($bi) = $depth/$opt_o;
push(@{$UPCAST{ENS}[$bi]},$ens);
push(@{$UPCAST{ELAPSED}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED});
push(@{$UPCAST{CTD_W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CTD_W});
push(@{$UPCAST{WINCH_W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{WINCH_W});
push(@{$UPCAST{WAVE_W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{WAVE_W});
push(@{$UPCAST{BIN}[$bi]},$bin);
push(@{$UPCAST{DEPTH}[$bi]},$depth);
push(@{$UPCAST{W}[$bi]},$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin]);
push(@{$UPCAST{PKGCORR_W}[10*round($LADCP{ENSEMBLE}[$ens]->{WAVE_W},0.1)+50]},
$LADCP{ENSEMBLE}[$ens]->{CORRECTED_OCEAN_W}[$bin])
if defined($opt_p) && abs($LADCP{ENSEMBLE}[$ens]->{WAVE_W}) < 5;
}
}
if (defined($opt_x)) {
progress("\t\tapplying wave_w correction...\n");
for (my($bi)=0; $bi<=$#{$UPCAST{ENS}}; $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]*$UPCAST{WAVE_W}[$bi][$i]**$e;
}
}
}
}
for (my($bi)=0; $bi<=$#{$UPCAST{ENS}}; $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]};
}
#---------------
# 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');
for (my($bi)=0; $bi<=max($#{$DNCAST{ENS}},$#{$UPCAST{ENS}}); $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]);
}
#-------------------------
# surface-wave effect file
#-------------------------
if (defined($opt_p)) {
progress("Writing surface-wave-correction data to <$opt_p>...\n");
@antsNewLayout = defined($opt_x) ?
('wave_w','dc_w','dc_w.mad','dc_w.N','dc_w_corr','uc_w','uc_w.mad','uc_w.N','uc_w_corr') :
('wave_w','dc_w','dc_w.mad','dc_w.N','uc_w','uc_w.mad','uc_w.N');
&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($wave_w) = ($bi-50) / 10;
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]*$wave_w**$e;
}
for (my($e)=0; $e<@uc_corr_poly; $e++) {
$uc_corr += $uc_corr_poly[$e]*$wave_w**$e;
}
&antsOut($wave_w,$dc_w,mad2($dc_w,@{$DNCAST{PKGCORR_W}[$bi]}),$dc_N,$dc_corr,
$uc_w,mad2($uc_w,@{$UPCAST{PKGCORR_W}[$bi]}),$uc_N,$uc_corr);
} else {
&antsOut($wave_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','elapsed_mismatch',
'depth','sound_speed','pitch','gimbal_pitch','roll','tilt','heading',
'CTD_w','CTD_w_t','CTD_winch_w','CTD_wave_w','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} - $LADCP{ENSEMBLE}[$ens]->{CTD_W}
: undef;
&antsOut($ens,
$LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED},
$LADCP{ENSEMBLE}[$ens]->{ELAPSED_MISMATCH},
$LADCP{ENSEMBLE}[$ens]->{DEPTH},
$LADCP{ENSEMBLE}[$ens]->{SOUND_SPEED},
$LADCP{ENSEMBLE}[$ens]->{PITCH},
$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH},
$LADCP{ENSEMBLE}[$ens]->{ROLL},
$LADCP{ENSEMBLE}[$ens]->{TILT},
$LADCP{ENSEMBLE}[$ens]->{HEADING},
$LADCP{ENSEMBLE}[$ens]->{CTD_W}, # NB: sound-speed corrected above
$LADCP{ENSEMBLE}[$ens]->{CTD_W_T},
$LADCP{ENSEMBLE}[$ens]->{WINCH_W},
$LADCP{ENSEMBLE}[$ens]->{WAVE_W},
$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) = $LADCP{ENSEMBLE}[$ens]->{SOUND_SPEED}/1500;
&antsOut($ens,
$bin,
$LADCP{ENSEMBLE}[$ens]->{CTD_ELAPSED},
$LADCP{ENSEMBLE}[$ens]->{ELAPSED_MISMATCH},
$depth,
$LADCP{ENSEMBLE}[$ens]->{SOUND_SPEED},
$LADCP{ENSEMBLE}[$ens]->{PITCH},
$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH},
$LADCP{ENSEMBLE}[$ens]->{ROLL},
$LADCP{ENSEMBLE}[$ens]->{TILT},
$LADCP{ENSEMBLE}[$ens]->{HEADING},
$LADCP{ENSEMBLE}[$ens]->{CTD_W},
$LADCP{ENSEMBLE}[$ens]->{WINCH_W},
$LADCP{ENSEMBLE}[$ens]->{WAVE_W},
$LADCP{ENSEMBLE}[$ens]->{CTD_W_T},
$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','elapsed_mismatch','depth','sound_speed','pitch','gimbal_pitch',
'roll','tilt','heading','CTD_w','CTD_winch_w','CTD_wave_w','CTD_w_t','LADCP_w','LADCP_errvel',
'LADCP_reflr_w','LADCP_reflr_w_err','w');
for (my($bi)=0; $bi<=$#{$DNCAST{ENS}}; $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{ENS}[$bi]}; $eii++) {
&outProfBinRec($DNCAST{ENS}[$bi][$eii],$DNCAST{BIN}[$bi][$eii],$DNCAST{DEPTH}[$bi][$eii]);
}
}
progress("Writing profile-bin data of upcast...\n");
for (my($bi)=0; $bi<=$#{$UPCAST{ENS}}; $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{ENS}[$bi]}; $eii++) {
&outProfBinRec($UPCAST{ENS}[$bi][$eii],$UPCAST{BIN}[$bi][$eii],$UPCAST{DEPTH}[$bi][$eii]);
}
close(STDOUT);
}
}
&antsExit();