V2.0
- lots of bug fixes
- major new features:
- dropped CTD scans handled correctly (no more apparent clock drifts)
- support for data files collected with Nortek Signature instruments
- much improved data editing
- significant changes:
- no minimum limit for eps_VKE
- updated for GMT6
- much better data quality information in summary plots
#!/usr/bin/perl
#======================================================================
# L A D C P _ V K E
# doc: Tue Oct 14 11:05:16 2014
# dlm: Sat Jul 24 09:44:26 2021
# (c) 2012 A.M. Thurnherr
# uE-Info: 120 45 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
# NOTES:
# - LADCP_VKE takes either one .wspec file (produced by [LADCP_wspec]
# and, possibly, post-processed) or any number of .wprof files (the
# normal case) as input
# - when .wprof files are supplied, [LADCP_wspec] is used to calculate
# the corresponding VKE spectra
# - when multiple .wprof files are supplied, the resulting spectra
# are averaged
# - the averaged spectra will only cover the windows that are present in
# the final(!) .wprof input files
# TODO:
# ! verify that p0fit.slope.sig is correct (-x scale factor)
# HISTORY:
# Oct 14, 2014: - created from [LADCPfs]
# Oct 15, 2014: - added parameterization output
# Oct 17, 2014: - changed parameterization constant $c to 0.021
# Nov 6, 2014: - restored from backup and adapted to ANTS V6.0
# Nov 7, 2014: - changed parameterization constant $c to 0.0215
# Apr 16, 2015: - disabled output activation unless ANTS tools are available
# - removed superfluous $ANTSBIN definition
# May 18, 2015: - added -p)ulse <length>
# Jun 11, 2015: - removed w_z code (yfname param requirement)
# Jun 12, 2015: - adapted to &antsParam()
# - BUG: %k.0 and %lambda.0 had been required erroneously (.1 are first ones to be used)
# - made finescale limits optional
# - renamed -b=>-l, -c=>-o, -s=>-i
# Jun 14, 2015: - renamed -p=>-a
# - added -p
# - removed weird evals
# Jun 15, 2015: - added plot & other mods
# Jun 16, 2015: - define default outputs when STDOUT is a tty
# - adapted to re-added pwrdens.0 in LADCP_wspec output
# - modified plot label
# Jun 17, 2015: - added eps.w to plot
# Jul 30, 2015: - moved main label and made it blue
# Sep 28, 2015: - added a p0fit r test to ensure that correlation has the correct sign
# Oct 12, 2015: - BUG: plot was in landscape mode
# - added input run label(s) to plot
# - require ANTSlibs V6.2 for release
# Oct 13, 2015: - adapted to [version.pl]
# Nov 30, 2015: - BUG: -a was not allowed in usage
# - added -f
# Dec 27, 2015: - reduced minlim for eps to 1e-13 W/kg
# Dec 29, 2015: - added 3rd consistency check (p0 limit)
# Jan 25, 2016: - added software version %PARAM
# Mar 7, 2016: - removed unused spectral normalization code
# - added latitude constraint for calculation of eps.w
# Mar 8, 2016: - renamed fs_pwr to pwr.fs for consistency
# - renamed eps.w to eps.VKE for consistency
# - changed default output filenames for -d and -u
# - removed ./ from figure label
# - update ANTSlib to 6.3
# Mar 16, 2016: - adapted to gmt5
# Mar 17, 2016: - added version to plot
# - added -k)e dissipation
# - reduced eps scale
# Mar 26, 2016: - added slope
# Mar 27, 2016: - added slope stddev
# - added opt_x
# - re-designed QC checks
# - added -z)ap
# Mar 28, 2016: - made QC tests consistent with I08S observations based on p0fit.rms == 0.4
# - added support for w.nsamp.avg
# - removed -z from LDCP_wspec
# - reduced -z default to 1, because 2 noise should affect only shotest scales
# Mar 29, 2016: - changed default of -x to 1 & removed from usage
# - BUG: -l/-a did not work
# - removed p0fit.nsamp from output (is constant)
# Mar 30, 2016: - added support for multiple wprof input files
# Mar 31, 2016: - standardized %PARAMs (got rid of old binpgrams:: params)
# - added spectral averaging
# - disabled -l/-a code (override bin size & pulse length)
# - added removal of params on multiple files
# - disabled code to create different plot/output names on -d/-u
# - added -q)uatorial cutoff lat
# Apr 1, 2016: - cosmetics
# Apr 2, 2016: - added low-p0 power-law fit as alternative
# Apr 3, 2016: - removed low-p0 power-law fit as it was not well thought out
# - changed -l default from 1e-7 to 5e-8, which agrees better with GHP
# solution without latiudinal dependency in case of 2016 I08S
# - define own LADCP_wspec defaults here
# - changed default of -q from 5 to 3
# Apr 6, 2016: - cosmetics
# Apr 16, 2016: - assertion failed on I98S#153 => modify code to allow nans in LADCP_wspec
# output
# Apr 27, 2016: - cosmetics (usage message)
# Aug 16, 2016: - increased -l default to 1.2e-7 based on UK2.5 SR1 repeat section
# Sep 1, 2016: - changed -l to mean epsilon, and increased value to 1e-10
# - added %eps.minlim
# Mar 13, 2017: - added -a)mbient <eps>
# Mar 14, 2017: - disabled -a) by default, because -a 0 is clearly bad and
# I have no evidence yet that -a something is better than -l 0
# Oct 17, 2017: - added default 'eps' field on -k
# - added eps.ms field to files processed without -k
# Dec 9, 2017: - added support for $antsSuppressCommonOptions
# Apr 24, 2018: - BUG: output was one field too wide (filled with nans) because antsBufNFields was not reset
# Apr 25, 2018: - added -y and removed spectral bins from default output
# Sep 12, 2019: - disabled default -l cut-off (used to be 5e-11 W/kg)
# Jul 1, 2021: - made %PARAMs more standard
# Jul 23, 2021: - returned to published empirical calibration constant (20% difference is not significant)
# - added calibration constant to output metadata
# - changed opt_a default from nan to 0 (ambient mixing)
# - BUG: -z default was only 1, which means that -k20 from LADCP_w_ocean applies; in
# practice I found with A20 that there are the effective value is >60; which means
# that this bug did not matter much in practice
# - changed calibration constant back, because it makes sense, and because it is now
# recorded in the meta-data
# HISTORY END
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
$WCALC = '.' if ($WCALC eq '');
$ANTS_TOOLS_AVAILABLE = (`which list` ne '');
require "$WCALC/version.pl";
require "$ANTSLIB/ants.pl";
require "$ANTSLIB/libLADCP.pl";
require "$ANTSLIB/libGMT.pl";
&antsAddParams('LADCP_VKE::version',$VERSION);
use FileHandle;
use IPC::Open2;
#----------------------------------------------------------------------
# Empirical constants & defaults
#----------------------------------------------------------------------
#my($c) = 0.0215; # Thurnherr et al. (GRL 2015)
my($c) = 0.026; # increased by 21% for V1.2beta7
$opt_q = 3; # Equatorial band: little more than a guess based on 2015 P16N
$opt_l = 0; # [W/kg]; cutoff disabled Sep 12, 2019
$opt_a = 0; # assume background dissipation for samples that pass the tests but have eps below -l
$opt_z = 50; # number of w_ocean samples to require (note that the .wprof inputs may have harsher limits)
$opt_o = 0; # remove mean before calculating spectra
$opt_s = 150; # surface layer to exclude from spectra
$opt_g = 40; # max gap to interpolate over
$opt_c_default = 100; # short-wavelength cutoff
#----------------------------------------------------------------------
# Usage
#----------------------------------------------------------------------
$antsSuppressCommonOptions = 1;
&antsUsage('a:bc:de:f:g:i:k:l:mno:p:q:r:s:tuw:x:yz:',0,
"[poly-o)rder <n[$opt_o]> to de-mean data; -1 to disable>] [apply cosine-t)aper]",
'[-d)own/-u)pcast-only] [exclude -b)ottom window]', # LADCP_wspec options
"[-s)urface <layer depth to exclude[${opt_s}m]>",
"[-g)ap <max depth layer to fill with interpolation[${opt_g}m]>]",
'[-w)indow <power-of-two input-records[32]>]',
"[shortwave -c)utoff <kz or lambda[${opt_c_default}m]>]", # LADCP_VKE options
"[e-q)uatorial cutoff <latitude[${opt_q}deg]>]",
"[-l)ow-eps <cutoff[${opt_l} W/kg]>] [-a)mbient <eps[${opt_a} W/kg]>]",
"[-z) ignore velocities derived from fewer than <N[$opt_z]> samples]",
'[o-m)it spectral correction] [spectral-tilt-correction -r)ange <max[0m]>]',
"[-e)ps-parameterization <constant[${c}s^-0.5]>",
'[include microstructure -k)e dissipation <file:field> in _VKE plot]',
'[-y) record spectra in output file]',
'[write output -f)iles to <directory>]',
'[write output filed with -i)ndividual spectra <basename>]',
'[output -p)lot <ps-file[#_VKE.ps]>]',
'[file...]');
#----------------------------------------------------------------------------
# Calculate VKE spectra with [LADCP_wspec] if input is a set of w_ocean files
#----------------------------------------------------------------------------
my($widx_min) = 99; # sentinels
my($widx_max) = -99;
if (defined(fnrNoErr('dc_w'))) { # pre-process with LADCP_wspec when handed vertical-velocity input
&antsInstallBufFull('eof'); # read first file
&antsIn();
my($opts); # set up options to pass
$opts .= ' -d' if ($opt_d);
$opts .= ' -u' if ($opt_u);
$opts .= ' -b' if ($opt_b);
$opts .= ' -t' if ($opt_t);
$opts .= " -s $opt_s" if defined($opt_s);
$opts .= " -g $opt_g" if defined($opt_g);
$opts .= " -w $opt_w" if defined($opt_w);
$opts .= " -o $opt_o" if defined($opt_o);
open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # spawn sub-process
croak("LADCP_wspec $opts: $!\n");
print(TOCLD $antsOldHeaders); # feed already gobbled header
for (my($r)=0; $r<@ants_; $r++) { # feed all .wprof records to LADCP_wspec
print(TOCLD "@{$ants_[$r]}\n");
}
for (my($bufi)=0; defined($ARGV[0]); $bufi++) { # multiple input files: loop until 2nd last
$input_list .= "$P{profile_id}($P{run_label}) ";
if ($bufi == 0) { # do once for mulitple files
&antsAddParams('ADCP_bin_length','','ADCP_blanking_distance','', # delete most %PARAMs, leaving
'ADCP_frequency','','ADCP_orientation','', # only potentially useful ones (from last file):
'ADCP_pulse_length','','BT_bin_range_diff.max','', # %profile_id
'BT_range.max','','BT_w_error.max','', # %water_depth
'BT_w_discrepancy.rms','','CTD_time_lags','', # %lat
'LADCP_firstBin','','LADCP_lastBin','', # %lon
'LADCP_w_ocean::version','','SS_samp.min','', # %dnXX
'SS_allowed_depth_range.max','','SS_signal.min','',
'Sv_ref_bin','','TL_allowed_three_lag_spread.max','',
'dc_pkg_accel.rms','','dc_tilt.rms','',
'uc_pkg_accel.rms','','uc_tilt.rms','',
'depth.max','','elapsed.max','','ens.max','',
'depth.min','','elapsed.min','','ens.min','',
'out_basename','','outgrid_dz','','run_label','',
'outgrid_firstbin','','outgrid_lastbin','',
'outgrid_minsamp','','per_bin_valid_frac_lim','',
'processing_options','','refLr_firstBin','',
'refLr_lastBin','','rms_w_reflr_err','',
'rms_w_reflr_err_interior','',
'sidelobe_editing','','surface_layer_depth','',
'vessel_draft','','w_lim.max','',
'water_depth.sig','','water_depth_from','',
);
}
close(TOCLD); # close LADCP_VKE input
my(@specrec);
while (@specrec = &antsFileIn(FROMCLD)) {
my($i) = ($specrec[$widxf]>0) ? $specrec[$widxf] : 0;
@{$specbuf[$bufi][$i]} = @specrec;
# print(STDERR "specrec[$i] = @specrec\n");
}
close(FROMCLD);
# printf(STDERR "windows up to %d buffered\n",$#{$specbuf[$bufi]});
$antsBufSkip = @ants_;
&antsIn(); # read next .wprof file
open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # process it
croak("LADCP_wspec $opts: $!\n");
print(TOCLD $antsOldHeaders);
for (my($r)=0; $r<@ants_; $r++) {
print(TOCLD "@{$ants_[$r]}\n");
}
}
close(TOCLD); # connect stdout from LADCP_VKE to stdin
open(STDIN,"<&FROMCLD") || croak("dup(FROMCLD): $!\n");
close(FROMCLD);
<>; # clear EOF condition
undef(%P); # shouldn't matter, because we'll get the same %PARAMs back
undef(@antsLayout); # shouldn't matter, because it will get overwritten
undef($antsOldHeaders); # forget those
undef(@ants_);
$antsBufNFields = 0;
} elsif (defined(fnrNoErr('pwrdens.0'))) {
croak("$0: -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n")
if ($opt_d || $opt_u || $opt_b || defined($opt_w) || defined($opt_s) || defined($opt_g));
} else {
if ($ARGV[0]) {
croak("$ARGV[0]: no such file or directory\n");
} else {
croak("$0: empyt input\n");
}
}
#----------------------------------------------------------------------
# Handle LADCP_VKE usage & read spectra from final file
# - spectra from previous files are in @specbuf
#----------------------------------------------------------------------
$n_input_files = 1 + @specbuf; # number of input files provided
&antsFloatOpt(\$opt_e,$c); # default parameterization
&antsFloatOpt(\$opt_x,1); # spectral fit stddev scale factor
&antsAddParams('LADCP_VKE::input_files.n',$n_input_files,
'LADCP_VKE::wsamp.min',$opt_z,
'LADCP_VKE::eps.minlim',$opt_l,
'LADCP_VKE::calibration_constant',$opt_e);
if (defined($opt_c)) { # shortwave cutoff supplied
$lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c;
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below
} elsif (defined(antsParam('shortwave_cutoff'))) { # cutoff already applied
$lmin = 2*$PI/antsParam('shortwave_cutoff');
} else { # use 100m default cutoff
$lmin = $opt_c_default;
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below
}
$lmax = 9e99; # no longwave cutoff implemented yet
&antsInstallBufFull(0); # load entire file
&antsIn();
$P{run_label} = "$input_list$P{profile_id}($P{run_label})"
if ($n_input_files > 1);
my($Hbuf) = $antsOldHeaders; # save for later (used on -i)
&antsRequireParam('profile_id');
&antsRequireParam('lambda.1');
&antsRequireParam('k.1');
&antsRequireParam('resolution_bandwidth');
&antsRequireParam('input_depth_resolution');
&antsRequireParam('output_depth_resolution');
&antsRequireParam('ADCP_bin_length');
&antsAddParams('ADCP_pulse_length',antsParam('ADCP_bin_length'))
unless defined(antsParam('ADCP_pulse_length'));
$imin = 0; # find frequency bin limits
for ($nfreq=1; defined(antsParam("lambda.$nfreq")); $nfreq++) {
$imin = $nfreq if ($imin==0 && antsParam("lambda.$nfreq")<=$lmax);
$imax = $nfreq if (antsParam("lambda.$nfreq") >= $lmin);
}
croak("$0: <lambda.min=$lmin> < min(lambda)")
unless defined($imax);
$pg_fmin = fnr('pwrdens.0'); # first power field in spectra
$fs_fmin = $pg_fmin + $imin; # first power field in finescale range
$fs_fmax = $pg_fmin + $imax; # last power field in finescale range
$widxf = fnr('widx');
$df = fnr('depth');
$mindf = fnr('depth.min');
$maxdf = fnr('depth.max');
$wsf = fnr('w.nsamp.avg');
$doff = fnr('nspec');
#----------------------------------------------------------------------
# Redirect STDOUT & create plot if STDOUT is a tty
#----------------------------------------------------------------------
if (-t STDOUT) {
$opt_f = '.' unless defined($opt_f); # output directory
croak("$opt_f: not a directory\n") unless (-d $opt_f);
my($id) = &antsRequireParam('profile_id');
# if ($opt_d) {
# $opt_p = sprintf('%s/%03d_dc_VKE.ps',$opt_f,$id) unless defined($opt_p);
# $outfile = sprintf('%s/%03d_dc.VKE',$opt_f,$id);
# } elsif ($opt_u) {
# $opt_p = sprintf('%s/%03d_uc_VKE.ps',$opt_f,$id) unless defined($opt_p);
# $outfile = sprintf('%s/%03d_uc.VKE',$opt_f,$id);
# } else {
$opt_p = sprintf('%s/%03d_VKE.ps',$opt_f,$id) unless defined($opt_p);
$outfile = sprintf('%s/%03d.VKE',$opt_f,$id);
# }
$outfile =~ s@^./@@;
open(STDOUT,">$outfile") || die("$outfile: $!\n");
} elsif (defined($opt_f)) {
croak("-f can only be used without STDOUT redirection\n");
}
#----------------------------------------------------------------------
# Library
#----------------------------------------------------------------------
sub average_spectra($)
{
my($r) = @_;
for (my($f)=$fs_fmin-1; $f<=$fs_fmax; $f++) { # average w.nsamp.avg & spectral densities
my($sum) = my($ns) = 0;
if (numberp($ants_[$r][$f])) { # final file has valid spectra
$sum = $ants_[$r][$f]; $ns = 1; # NB: nspec is correct even if it doesn't
}
my($wi) = ($ants_[$r][$widxf]>0) ? $ants_[$r][$widxf] : 0; # adjust for -1
for (my($bi)=0; $bi<@specbuf; $bi++) { # loop over all buffered files
next unless @{$specbuf[$bi][$wi]}; # skip input files w/o valid spectra
if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) > 0) { # depth mismatch
# die("assertion failed ($specbuf[$bi][$wi][$df] - $ants_[$r][$df] @ wi=$wi)") unless ($wi == 0); # only allowed in bottom window
if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) >
abs($ants_[$r][$maxdf] - $ants_[$r][$mindf])) {
printf(STDERR "WARNING: ignoring window #$wi from input file #%d because of depth mismatch\n",$bi+1)
if ($f == $fs_fmin);
next;
}
}
if (numberp($specbuf[$bi][$wi][$f])) {
$sum += $specbuf[$bi][$wi][$f]; $ns++;
}
$ants_[$r][$doff] += $specbuf[$bi][$wi][$doff] # update nspec once per input record
if ($f == $fs_fmax); # ... but for all files
}
$ants_[$r][$f] = ($ns > 0) ? $sum/$ns : nan; # update averaged spectral density
}
}
sub integrate_fs_power($) # integrate fs spectrum
{
my($r) = @_;
$ants_[$r][$fspwrf] = 0;
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
$ants_[$r][$fspwrf] += $ants_[$r][$f];
}
$ants_[$r][$fspwrf] *= antsParam('resolution_bandwidth');
}
sub fit_universal_w_spec($) # vertical velocity => p0
{
my($r) = @_;
my($nsamp) = $fs_fmax - $fs_fmin + 1;
#---------------------------------------------------
# fit slope-2 line in log-log space (main estimator)
#---------------------------------------------------
if ($nsamp >= 2) { # require min 2 wavenumber samples
my($DOF) = 0;
my($sumd,$sumx,$sumy) = (0,0,0); # fit kz^-2 power law
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($i) = $f - $pg_fmin;
$sumx += log10(antsParam("k.$i"));
$sumy += log10($ants_[$r][$f]);
$sumd += log10($ants_[$r][$f]) + 2*log10(antsRequireParam("k.$i"));
}
my($p0) = $sumd/$nsamp;
$ants_[$r][$p0f] = 10**$p0;
my($avgx) = $sumx/$nsamp; # avg for r calc
my($avgy) = $sumy/$nsamp;
my($sumsqerr,$sxx,$syy,$sxy,$sumsqxt,$sumx,$sumy,$sumwt) = # r, rms error, pwrlaw slope
(0,0,0,0,0,0,0,0);
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($i) = $f - $pg_fmin;
my($x) = log10(&antsParam("k.$i"));
my($y) = log10($ants_[$r][$f]);
my($ysig) = $opt_x * $y / sqrt($ants_[$r][$doff]);
my($xt) = $x - $avgx; $sxx += &SQR($xt); # correlation coeff (r)
my($yt) = $y - $avgy; $syy += &SQR($yt); $sxy += $xt * $yt;
my($wt) = 1 / &SQR($ysig); $sumwt += $wt; # slope (linear fit in log-log space)
$sumx += $x * $wt;
$sumy += $y * $wt;
$sumsqerr += ($p0 - 2*$x - $y)**2; # rms error
}
my($midx) = $sumx / $sumwt;
# print(STDERR "$sumx:$sumy:$sumwt\n");
# print(STDERR "$midx\n");
my($sumsqdx,$sumslp) = (0,0);
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($i) = $f - $pg_fmin;
my($x) = log10(&antsParam("k.$i"));
my($y) = log10($ants_[$r][$f]);
my($ysig) = $opt_x * $y / sqrt($ants_[$r][$doff]);
my($dx) = ($x - $midx) / $ysig;
$sumsqdx += &SQR($dx);
$sumslp += $dx * $y / $ysig;
}
$ants_[$r][$rmsf] = sqrt($sumsqerr/$nsamp);
$ants_[$r][$rf] = $sxy/(sqrt($sxx * $syy) + $SMALL_AMOUNT);
$ants_[$r][$slpf] = $sumslp / $sumsqdx;
$ants_[$r][$sslpf] = sqrt(1 / $sumsqdx);
} else {
&antsInfo("WARNING: no fit --- need min 2 samples");
}
}
#----------------------------------------------------------------------
# Load Dissipation Data
#----------------------------------------------------------------------
my(@eps_ms,@depth_ms); # output variables
if (defined($opt_k)) {
my($file,$field) = split(':',$opt_k);
$field = 'eps' unless defined($field);
open(ADDF,"$file") || croak("$file: $!\n"); # open file
my(@afl) = &antsFileLayout(ADDF); # read layout
my($akf,$aef);
for (my($f)=0; $f<=$#afl; $f++) { # find depth & eps fields
$akf = $f if ($afl[$f] eq 'depth');
$aef = $f if ($afl[$f] eq $field);
}
croak("$file: fields 'depth' and '$field' required\n")
unless defined($akf) && defined($aef);
while (1) { # load entire profile
my(@ar);
last unless @ar = &antsFileIn(ADDF);
next unless numberp($ar[$akf]) && numberp($ar[$aef]);
push(@depth_ms,$ar[$akf]); push(@eps_ms,$ar[$aef]);
}
close(ADDF);
}
#----------------------------------------------------------------------
# Process File
#----------------------------------------------------------------------
$fspwrf = &antsNewField('pwr.fs'); # derived fields
$p0f = &antsNewField('p0'); # VKE density
$rmsf = &antsNewField('p0fit.rms'); # rms misfit
$rf = &antsNewField('p0fit.r'); # correlation coefficient
$slpf = &antsNewField('p0fit.slope'); # power-law slope
$sslpf = &antsNewField('p0fit.slope.sig'); # power-law slope stddev
$wepsf = &antsNewField('eps.VKE'); # epsilon from VKE
$msepsf = &antsNewField('eps.ms'); # externally supplied microstructure eps if available
my(@outLayout) = @antsNewLayout; # save for later
for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input
$totf = $f if ($outLayout[$f] eq 'pwr.tot');
$tnsf = $f if ($outLayout[$f] eq 'pwr.tot.nsamp');
}
croak("$0: cannot find fields 'pwr.tot' or 'pwr.tot.nsamp' in input\n")
unless defined($totf) || defined($tnsf);
$lsf = defined($tnsf) ? $tnsf : $totf;
if ($opt_p) { # begin plot
$xmin = (&antsParam('output_depth_resolution')>=350) ? 0.012 : 0.018;
$xmax = 0.2; $ymin = 1; $ymax = 1e4;
$plotsize = 13;
GMT_begin($opt_p,"-JX${plotsize}l","-R$xmin/$xmax/$ymin/$ymax",'-P -X6 -Y4'); # init plot
GMT_psxy('-G25/255/25 -L');
printf(GMT "%f %f\n",$xmin,$xmin**(-2)*sqrt(2));
printf(GMT "%f %f\n",$xmax,$xmax**(-2)*sqrt(2));
printf(GMT "%f %f\n",$xmax,$xmax**(-2)/sqrt(2));
printf(GMT "%f %f\n",$xmin,$xmin**(-2)/sqrt(2));
}
my(@sumPwr,@nPwr);
my(@sumGoodPwr,@nGoodPwr);
my($min_depth) = 9e99;
my($max_depth) = -9e99;
my($latM) = abs(&antsRequireParam('lat'));
&antsInfo("WARNING: low latitude-profile no epsilon estimated")
unless ($latM > $opt_q);
for (my($r)=0; $r<@ants_; $r++) { # loop over all windows
if (numberp($ants_[$r][$pg_fmin])) { # there is a spectrum
#--------------------------
# apply spectral correction
#--------------------------
unless ($opt_m) {
for (my($i)=0; $i<$nfreq; $i++) { # loop over wavenumbers
$ants_[$r][$i+$pg_fmin] *=
T_w(antsParam("k.$i"),antsParam('ADCP_bin_length'),
antsParam('ADCP_pulse_length'),antsParam('input_depth_resolution'),
$opt_r);
}
}
#------------------------
# calculate fs quantities
#------------------------
average_spectra($r); # average all avaiable spectra
integrate_fs_power($r); # calculate total finescale power
fit_universal_w_spec($r); # fit kz^-2 spectrum & calc stats
if (numberp($ants_[$r][$p0f])) { # update min/max depth
$min_depth = $ants_[$r][$mindf] if ($ants_[$r][$mindf] < $min_depth);
$max_depth = $ants_[$r][$maxdf] if ($ants_[$r][$maxdf] > $max_depth);
}
} else { # no spectrum
$ants_[$r][$fspwrf] = $ants_[$r][$p0f] = $ants_[$r][$rmsf] =
$ants_[$r][$rf] = $ants_[$r][$slpf] = $ants_[$r][$sslpf] = nan;
}
#-----------------------------------------------------------------------------------------------------
# eps.VKE QC Tests:
# - the following limits were independently derived
# p0fit.rms <= 0.4 primary filter used in Thurnherr et al. (GRL 2015)
# -3 <= p0fit.slope <= -1 based largely on 2016 I08S data with sufficient/insufficient range
# p0fit.r <= -0.5 based largely on 2016 I08S data with sufficient/insufficient range
# w.nsamp.avg >= 50 based on observations in many data sets with weak backscatter,
# including DoMORE, GO-SHIP P16S, GO-SHIP I08S
# - then, I plotted slope & r vs. rms and found that rms = 0.4 corresponds to, on average
# -3 <= p0fit.slope <= -1
# p0fit.r <= -0.4
# - in a plot of rms vs nsamp the limiting value of 0.4 is hit at 50 samples
#
# => SET OF MUTUALLY CONSISTENT CRITERIA
#
# Additional Empirical Filters:
# - latitude > 3deg guess based on Thurnherr et al., 2015, Gregg et al., 2003, 2015 GO-SHIP P16N
# - eps >= 5e-11 W/kg based on DIMES data, cutting where errors become > factor 2; DISABLED SEPTEMBER 2019
#-----------------------------------------------------------------------------------------------------
if ($latM > $opt_q && # 1) not (too) equatorial
$ants_[$r][$rmsf] <= 0.4 && # 2) rms spectra misfit <= 0.4 (as in Thurnherr et al., GRL 2015)
$ants_[$r][$slpf]>=-3 && $ants_[$r][$slpf]<=-1 && # 3) slope consistent with -2 power law
$ants_[$r][$rf] <= -0.4 && # 4) p and k_z are well correlated
$ants_[$r][$wsf] >= $opt_z) { # 5) minimum # of samples
if (($ants_[$r][$p0f]/$opt_e)**2 >= $opt_l) { # level is above eps.minlim (-l)
$ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2; # => Thurnherr et al. (GRL 2015)
} else { # level is below eps.minlim
$ants_[$r][$wepsf] = $opt_a; # => set to arbitrary background value
}
} else { # failed any of checks 1-5 => nan
$ants_[$r][$wepsf] = nan;
}
#-------------------------------------------------
# average external microstructure eps, if supplied
#-------------------------------------------------
if (defined($opt_k)) {
my($sum,$n) = (0,0);
for (my($i)=0; $i<@eps_ms; $i++) { # linearly search all eps records
next unless $depth_ms[$i] >= $ants_[$r][$mindf] &&
$depth_ms[$i] <= $ants_[$r][$maxdf];
$sum += $eps_ms[$i]; $n++;
}
$ants_[$r][$msepsf] = $n ? $sum / $n : nan;
} else {
$ants_[$r][$msepsf] = nan;
}
#---------------
# produce output
#---------------
if (defined($opt_p)) { # plot current spectrum on -p & calc mean
my($R) = 0; # RGB map
my($G) = int(200*(1-$r/@ants_));
my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_));
GMT_psxy("-W1,$R/$G/$B");
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { # avg & plot high-quality spectra only
next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]); # - omit zero wavenumber ($pg_fmin)
next unless numberp($ants_[$r][$wepsf]);
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]);
$nGoodPwr[$f]++;
$sumGoodPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f];
}
GMT_psxy("-W1,$R/$G/$B,."); # avg & plot all spectra with dots
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) {
next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]);
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]);
$nPwr[$f]++;
$sumPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f];
}
}
if (defined($opt_i)) { # output current spectrum on -i
open(STDOUT_DUP,">&",STDOUT) || croak("$0: cannot dup STDOUT\n");
@antsNewLayout = ('k','lambda','pwrdens','finescale','pwrdens_fit');
if ($ants_[$r][$widxf] > 0) {
$ofname = sprintf("${opt_i}_%02d.wspec",$r+1);
} else {
$ofname = sprintf("${opt_i}_btm.wspec");
}
open(STDOUT,">$ofname") || croak("$ofname: $!\n");
undef($antsOldHeaders);
&antsActivateOut() if ($ANTS_TOOLS_AVAILABLE);
my($saveParams) = $antsCurParams; # add all extra input fields as %PARAMs
&antsAddParams('LADCP_VKE::depth.min',$ants_[$r][$mindf],
'LADCP_VKE::depth.max',$ants_[$r][$maxdf]);
for (my($f)=$lsf+1; $f<@outLayout; $f++) {
&antsAddParams("LADCP_VKE::$outLayout[$f]",$ants_[$r][$f]);
}
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) {
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
my($l) = antsParam(sprintf("lambda.%d",$f-$pg_fmin));
&antsOut($k,$l,$ants_[$r][$f],
($f>=$fs_fmin && $f<=$fs_fmax),
numberp($ants_[$r][$p0f]) ? $ants_[$r][$p0f] * $k**(-2) : nan);
}
&antsOut('EOF');
open(STDOUT,">&",STDOUT_DUP) || croak("$0: cannot restore STDOUT\n");
close(STDOUT_DUP);
$antsCurParams = $saveParams;
}
}
if (defined($opt_p)) { # complete plot
if (@nPwr) {
GMT_psxy('-W2,255/100/100,.'); # plot mean spectrum; dotted; entire range
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
printf(GMT "$k %g\n",$sumPwr[$f]/$nPwr[$f]);
}
}
if (@nGoodPwr) {
GMT_psxy('-W1,255/100/100'); # plot mean spectrum; entire range
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) {
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]);
}
GMT_psxy('-W4,255/100/100'); # plot mean fs spectrum (heavy)
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]);
}
}
GMT_psbasemap('-Bf3a2:"Vertical Wavenumber [rad/m]":/' . # annotate axes
'f3a1:"VKE/p@-0@- [m@+3@+s@+-1/2@+/(rad m@+-1@+)]":WeS');
GMT_setJ("-JX-${plotsize}l");
GMT_setR(sprintf('-R%f/%f/%f/%f',2*$PI/$xmax,2*$PI/$xmin,$ymin,$ymax));
GMT_psbasemap('-B3:"Vertical Wavelength [m]:N"');
GMT_unitcoords_logscale(); # print profile number
GMT_pstext('-F+f14,Helvetica,blue+jTL -Gwhite');
if (defined($outfile)) { printf(GMT "0.02 0.98 $outfile [$P{run_label}] %s\n",antsParam('input_data')); }
else { printf(GMT "0.02 0.98 %03d [$P{run_label}] %s\n",antsParam('profile_id'),antsParam('input_data')); }
GMT_pstext('-F+f9,Helvetica,orange+jTR -N -Gwhite');
print(GMT "0.99 0.99 V$VERSION\n");
GMT_set('FONT_ANNOT_PRIMARY 10','MAP_ANNOT_OFFSET_PRIMARY 0.01c', # eps profile inset
'FONT_LABEL 10','MAP_LABEL_OFFSET 0.2c');
$min_depth = 400 if ($min_depth > 1e4);
$max_depth = 1500 if ($max_depth < 0);
GMT_setR(sprintf("-R1e-11/1e-7/%d/%d",round($min_depth-250,500),round($max_depth+250,500)));
GMT_setJ(sprintf('-JX%fl/-%f',$plotsize/2.2,$plotsize/2.2));
GMT_psbasemap('-X0.3 -Y0.3 -Bg1a1f3p:"@~e@~@-VKE@- [W kg@+-1@+]":/g1000a500f100:"Depth [m]":wEsN');
GMT_psxy();
for (my($r)=0; $r<@ants_; $r++) {
my($R) = 0; my($G) = int(200*(1-$r/@ants_)); # calculate color ramp
my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_));
if (numberp($ants_[$r][$msepsf])) { # microstructure eps
print(GMT "> -W1.5,DarkOrange2\n");
print(GMT "$ants_[$r][$msepsf] $ants_[$r][$mindf]\n");
print(GMT "$ants_[$r][$msepsf] $ants_[$r][$maxdf]\n");
}
if (numberp($ants_[$r][$wepsf])) {
print(GMT "> -W2,$R/$G/$B\n"); # plot eps.w in blue
print(GMT "$ants_[$r][$wepsf] $ants_[$r][$mindf]\n");
print(GMT "$ants_[$r][$wepsf] $ants_[$r][$maxdf]\n");
}
}
GMT_end(); # finish plot
}
@antsNewLayout = @outLayout; # restore layout
$antsOldHeaders = $Hbuf;
$antsHeadersPrinted = 0;
unless (defined($opt_y)) { # remove spectral bins from output
splice(@antsNewLayout,$pg_fmin,$nfreq);
for (my($r)=0; $r<@ants_; $r++) {
splice(@{$ants_[$r]},$pg_fmin,$nfreq);
}
}
&antsFlush(); # output results
&antsExit();