author | Andreas Thurnherr <ant@ldeo.columbia.edu> |
Mon, 27 Jun 2022 20:07:11 -1000 | |
changeset 59 | 4118a8e880de |
parent 57 | 69e39fcb7f41 |
permissions | -rwxr-xr-x |
18 | 1 |
#!/usr/bin/perl |
2 |
#====================================================================== |
|
3 |
# L A D C P _ V K E |
|
4 |
# doc: Tue Oct 14 11:05:16 2014 |
|
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
5 |
# dlm: Tue May 17 16:46:00 2022 |
18 | 6 |
# (c) 2012 A.M. Thurnherr |
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
7 |
# uE-Info: 126 87 NIL 0 0 72 2 2 4 NIL ofnI |
18 | 8 |
#====================================================================== |
9 |
||
10 |
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles'; |
|
11 |
||
43 | 12 |
# NOTES: |
13 |
# - LADCP_VKE takes either one .wspec file (produced by [LADCP_wspec] |
|
14 |
# and, possibly, post-processed) or any number of .wprof files (the |
|
15 |
# normal case) as input |
|
16 |
# - when .wprof files are supplied, [LADCP_wspec] is used to calculate |
|
17 |
# the corresponding VKE spectra |
|
18 |
# - when multiple .wprof files are supplied, the resulting spectra |
|
19 |
# are averaged |
|
20 |
# - the averaged spectra will only cover the windows that are present in |
|
21 |
# the final(!) .wprof input files |
|
22 |
||
23 |
# TODO: |
|
42 | 24 |
# ! verify that p0fit.slope.sig is correct (-x scale factor) |
25 |
||
18 | 26 |
# HISTORY: |
27 |
# Oct 14, 2014: - created from [LADCPfs] |
|
28 |
# Oct 15, 2014: - added parameterization output |
|
29 |
# Oct 17, 2014: - changed parameterization constant $c to 0.021 |
|
30 |
# Nov 6, 2014: - restored from backup and adapted to ANTS V6.0 |
|
31 |
# Nov 7, 2014: - changed parameterization constant $c to 0.0215 |
|
26
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
32 |
# Apr 16, 2015: - disabled output activation unless ANTS tools are available |
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
33 |
# - removed superfluous $ANTSBIN definition |
29 | 34 |
# May 18, 2015: - added -p)ulse <length> |
35 |
# Jun 11, 2015: - removed w_z code (yfname param requirement) |
|
36 |
# Jun 12, 2015: - adapted to &antsParam() |
|
37 |
# - BUG: %k.0 and %lambda.0 had been required erroneously (.1 are first ones to be used) |
|
38 |
# - made finescale limits optional |
|
39 |
# - renamed -b=>-l, -c=>-o, -s=>-i |
|
40 |
# Jun 14, 2015: - renamed -p=>-a |
|
41 |
# - added -p |
|
42 |
# - removed weird evals |
|
43 |
# Jun 15, 2015: - added plot & other mods |
|
44 |
# Jun 16, 2015: - define default outputs when STDOUT is a tty |
|
45 |
# - adapted to re-added pwrdens.0 in LADCP_wspec output |
|
46 |
# - modified plot label |
|
47 |
# Jun 17, 2015: - added eps.w to plot |
|
31 | 48 |
# Jul 30, 2015: - moved main label and made it blue |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
49 |
# Sep 28, 2015: - added a p0fit r test to ensure that correlation has the correct sign |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
50 |
# Oct 12, 2015: - BUG: plot was in landscape mode |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
51 |
# - added input run label(s) to plot |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
52 |
# - require ANTSlibs V6.2 for release |
33 | 53 |
# Oct 13, 2015: - adapted to [version.pl] |
54 |
# Nov 30, 2015: - BUG: -a was not allowed in usage |
|
55 |
# - added -f |
|
56 |
# Dec 27, 2015: - reduced minlim for eps to 1e-13 W/kg |
|
34 | 57 |
# Dec 29, 2015: - added 3rd consistency check (p0 limit) |
58 |
# Jan 25, 2016: - added software version %PARAM |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
59 |
# Mar 7, 2016: - removed unused spectral normalization code |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
60 |
# - added latitude constraint for calculation of eps.w |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
61 |
# Mar 8, 2016: - renamed fs_pwr to pwr.fs for consistency |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
62 |
# - renamed eps.w to eps.VKE for consistency |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
63 |
# - changed default output filenames for -d and -u |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
64 |
# - removed ./ from figure label |
36 | 65 |
# - update ANTSlib to 6.3 |
41 | 66 |
# Mar 16, 2016: - adapted to gmt5 |
42 | 67 |
# Mar 17, 2016: - added version to plot |
68 |
# - added -k)e dissipation |
|
69 |
# - reduced eps scale |
|
70 |
# Mar 26, 2016: - added slope |
|
71 |
# Mar 27, 2016: - added slope stddev |
|
72 |
# - added opt_x |
|
73 |
# - re-designed QC checks |
|
74 |
# - added -z)ap |
|
75 |
# Mar 28, 2016: - made QC tests consistent with I08S observations based on p0fit.rms == 0.4 |
|
76 |
# - added support for w.nsamp.avg |
|
77 |
# - removed -z from LDCP_wspec |
|
78 |
# - reduced -z default to 1, because 2 noise should affect only shotest scales |
|
79 |
# Mar 29, 2016: - changed default of -x to 1 & removed from usage |
|
80 |
# - BUG: -l/-a did not work |
|
81 |
# - removed p0fit.nsamp from output (is constant) |
|
43 | 82 |
# Mar 30, 2016: - added support for multiple wprof input files |
83 |
# Mar 31, 2016: - standardized %PARAMs (got rid of old binpgrams:: params) |
|
84 |
# - added spectral averaging |
|
85 |
# - disabled -l/-a code (override bin size & pulse length) |
|
86 |
# - added removal of params on multiple files |
|
87 |
# - disabled code to create different plot/output names on -d/-u |
|
88 |
# - added -q)uatorial cutoff lat |
|
89 |
# Apr 1, 2016: - cosmetics |
|
90 |
# Apr 2, 2016: - added low-p0 power-law fit as alternative |
|
91 |
# Apr 3, 2016: - removed low-p0 power-law fit as it was not well thought out |
|
92 |
# - changed -l default from 1e-7 to 5e-8, which agrees better with GHP |
|
93 |
# solution without latiudinal dependency in case of 2016 I08S |
|
94 |
# - define own LADCP_wspec defaults here |
|
95 |
# - changed default of -q from 5 to 3 |
|
96 |
# Apr 6, 2016: - cosmetics |
|
44 | 97 |
# Apr 16, 2016: - assertion failed on I98S#153 => modify code to allow nans in LADCP_wspec |
98 |
# output |
|
99 |
# Apr 27, 2016: - cosmetics (usage message) |
|
48 | 100 |
# Aug 16, 2016: - increased -l default to 1.2e-7 based on UK2.5 SR1 repeat section |
101 |
# Sep 1, 2016: - changed -l to mean epsilon, and increased value to 1e-10 |
|
102 |
# - added %eps.minlim |
|
103 |
# Mar 13, 2017: - added -a)mbient <eps> |
|
104 |
# Mar 14, 2017: - disabled -a) by default, because -a 0 is clearly bad and |
|
105 |
# I have no evidence yet that -a something is better than -l 0 |
|
49 | 106 |
# Oct 17, 2017: - added default 'eps' field on -k |
107 |
# - added eps.ms field to files processed without -k |
|
108 |
# Dec 9, 2017: - added support for $antsSuppressCommonOptions |
|
109 |
# Apr 24, 2018: - BUG: output was one field too wide (filled with nans) because antsBufNFields was not reset |
|
110 |
# Apr 25, 2018: - added -y and removed spectral bins from default output |
|
54 | 111 |
# Sep 12, 2019: - disabled default -l cut-off (used to be 5e-11 W/kg) |
56 | 112 |
# Jul 1, 2021: - made %PARAMs more standard |
113 |
# Jul 23, 2021: - returned to published empirical calibration constant (20% difference is not significant) |
|
114 |
# - added calibration constant to output metadata |
|
115 |
# - changed opt_a default from nan to 0 (ambient mixing) |
|
116 |
# - BUG: -z default was only 1, which means that -k20 from LADCP_w_ocean applies; in |
|
117 |
# practice I found with A20 that there are the effective value is >60; which means |
|
118 |
# that this bug did not matter much in practice |
|
119 |
# - changed calibration constant back, because it makes sense, and because it is now |
|
120 |
# recorded in the meta-data |
|
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
121 |
# Oct 19, 2021: - added -h) (min samples) |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
122 |
# - disabled spectral tests based on data editing implemented for GO-SHIP A20 |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
123 |
# - this allowed return to original $c value for VKE parameterization |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
124 |
# - doubled default surface layer from 150 to 300m |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
125 |
# May 17. 2022: - BUG: -h was wrong (used # of files instead of # of spectra) |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
126 |
# - changed semantics to make output file name based on input file name |
56 | 127 |
# HISTORY END |
18 | 128 |
|
129 |
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$}); |
|
33 | 130 |
($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
131 |
$WCALC = '.' if ($WCALC eq ''); |
|
26
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
132 |
$ANTS_TOOLS_AVAILABLE = (`which list` ne ''); |
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
133 |
|
33 | 134 |
require "$WCALC/version.pl"; |
18 | 135 |
require "$ANTSLIB/ants.pl"; |
136 |
require "$ANTSLIB/libLADCP.pl"; |
|
29 | 137 |
require "$ANTSLIB/libGMT.pl"; |
43 | 138 |
&antsAddParams('LADCP_VKE::version',$VERSION); |
29 | 139 |
|
140 |
use FileHandle; |
|
141 |
use IPC::Open2; |
|
18 | 142 |
|
143 |
#---------------------------------------------------------------------- |
|
43 | 144 |
# Empirical constants & defaults |
42 | 145 |
#---------------------------------------------------------------------- |
146 |
||
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
147 |
my($c) = 0.0215; # Thurnherr et al. (GRL 2015) |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
148 |
#my($c) = 0.026; # increased by 21% for V1.2beta7 when spectral filters were introduced; disabled Oct 2021 when filters were disabled vor V2.1 |
43 | 149 |
$opt_q = 3; # Equatorial band: little more than a guess based on 2015 P16N |
54 | 150 |
$opt_l = 0; # [W/kg]; cutoff disabled Sep 12, 2019 |
56 | 151 |
$opt_a = 0; # assume background dissipation for samples that pass the tests but have eps below -l |
152 |
$opt_z = 50; # number of w_ocean samples to require (note that the .wprof inputs may have harsher limits) |
|
43 | 153 |
$opt_o = 0; # remove mean before calculating spectra |
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
154 |
$opt_s = 300; # surface layer to exclude from spectra |
43 | 155 |
$opt_g = 40; # max gap to interpolate over |
156 |
$opt_c_default = 100; # short-wavelength cutoff |
|
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
157 |
$opt_h = 2; # require 2 spectra for estimate |
42 | 158 |
|
159 |
#---------------------------------------------------------------------- |
|
18 | 160 |
# Usage |
161 |
#---------------------------------------------------------------------- |
|
162 |
||
49 | 163 |
$antsSuppressCommonOptions = 1; |
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
164 |
&antsUsage('a:bc:de:f:g:h:i:k:l:mno:p:q:r:S:s:tuw:x:yz:',0, |
43 | 165 |
"[poly-o)rder <n[$opt_o]> to de-mean data; -1 to disable>] [apply cosine-t)aper]", |
29 | 166 |
'[-d)own/-u)pcast-only] [exclude -b)ottom window]', # LADCP_wspec options |
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
167 |
"[require spectral -s)amples <min[${opt_h}]>", |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
168 |
"[-S)urface <layer depth to exclude[${opt_s}m]>", |
43 | 169 |
"[-g)ap <max depth layer to fill with interpolation[${opt_g}m]>]", |
29 | 170 |
'[-w)indow <power-of-two input-records[32]>]', |
44 | 171 |
"[shortwave -c)utoff <kz or lambda[${opt_c_default}m]>]", # LADCP_VKE options |
48 | 172 |
"[e-q)uatorial cutoff <latitude[${opt_q}deg]>]", |
173 |
"[-l)ow-eps <cutoff[${opt_l} W/kg]>] [-a)mbient <eps[${opt_a} W/kg]>]", |
|
43 | 174 |
"[-z) ignore velocities derived from fewer than <N[$opt_z]> samples]", |
29 | 175 |
'[o-m)it spectral correction] [spectral-tilt-correction -r)ange <max[0m]>]', |
43 | 176 |
"[-e)ps-parameterization <constant[${c}s^-0.5]>", |
48 | 177 |
'[include microstructure -k)e dissipation <file:field> in _VKE plot]', |
49 | 178 |
'[-y) record spectra in output file]', |
48 | 179 |
'[write output -f)iles to <directory>]', |
180 |
'[write output filed with -i)ndividual spectra <basename>]', |
|
181 |
'[output -p)lot <ps-file[#_VKE.ps]>]', |
|
43 | 182 |
'[file...]'); |
29 | 183 |
|
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
184 |
#---------------------------------------------------------------------- |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
185 |
# Determine Output File Name From Input File Names |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
186 |
#---------------------------------------------------------------------- |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
187 |
|
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
188 |
unless ($ARGV eq '-') { |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
189 |
($basename) = ($ARGV =~ m{([^/]*)\.[^\.]*$}); |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
190 |
for (my($i)=0; $i<@ARGV; $i++) { |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
191 |
my($basename2) = ($ARGV[$i] =~ m{([^/]*)\.[^\.]*$}); |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
192 |
next if ($basename2 eq $basename); |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
193 |
undef($basename); |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
194 |
last; |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
195 |
} |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
196 |
} |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
197 |
|
43 | 198 |
#---------------------------------------------------------------------------- |
199 |
# Calculate VKE spectra with [LADCP_wspec] if input is a set of w_ocean files |
|
200 |
#---------------------------------------------------------------------------- |
|
42 | 201 |
|
43 | 202 |
my($widx_min) = 99; # sentinels |
203 |
my($widx_max) = -99; |
|
18 | 204 |
|
42 | 205 |
if (defined(fnrNoErr('dc_w'))) { # pre-process with LADCP_wspec when handed vertical-velocity input |
43 | 206 |
&antsInstallBufFull('eof'); # read first file |
207 |
&antsIn(); |
|
208 |
||
209 |
my($opts); # set up options to pass |
|
210 |
$opts .= ' -d' if ($opt_d); |
|
29 | 211 |
$opts .= ' -u' if ($opt_u); |
212 |
$opts .= ' -b' if ($opt_b); |
|
213 |
$opts .= ' -t' if ($opt_t); |
|
214 |
$opts .= " -s $opt_s" if defined($opt_s); |
|
215 |
$opts .= " -g $opt_g" if defined($opt_g); |
|
216 |
$opts .= " -w $opt_w" if defined($opt_w); |
|
217 |
$opts .= " -o $opt_o" if defined($opt_o); |
|
43 | 218 |
|
42 | 219 |
open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # spawn sub-process |
29 | 220 |
croak("LADCP_wspec $opts: $!\n"); |
43 | 221 |
print(TOCLD $antsOldHeaders); # feed already gobbled header |
222 |
for (my($r)=0; $r<@ants_; $r++) { # feed all .wprof records to LADCP_wspec |
|
223 |
print(TOCLD "@{$ants_[$r]}\n"); |
|
224 |
} |
|
225 |
||
226 |
for (my($bufi)=0; defined($ARGV[0]); $bufi++) { # multiple input files: loop until 2nd last |
|
227 |
$input_list .= "$P{profile_id}($P{run_label}) "; |
|
228 |
if ($bufi == 0) { # do once for mulitple files |
|
229 |
&antsAddParams('ADCP_bin_length','','ADCP_blanking_distance','', # delete most %PARAMs, leaving |
|
230 |
'ADCP_frequency','','ADCP_orientation','', # only potentially useful ones (from last file): |
|
56 | 231 |
'ADCP_pulse_length','','BT_bin_range_diff.max','', # %profile_id |
232 |
'BT_range.max','','BT_w_error.max','', # %water_depth |
|
233 |
'BT_w_discrepancy.rms','','CTD_time_lags','', # %lat |
|
43 | 234 |
'LADCP_firstBin','','LADCP_lastBin','', # %lon |
56 | 235 |
'LADCP_w_ocean::version','','SS_samp.min','', # %dnXX |
236 |
'SS_allowed_depth_range.max','','SS_signal.min','', |
|
237 |
'Sv_ref_bin','','TL_allowed_three_lag_spread.max','', |
|
238 |
'dc_pkg_accel.rms','','dc_tilt.rms','', |
|
239 |
'uc_pkg_accel.rms','','uc_tilt.rms','', |
|
240 |
'depth.max','','elapsed.max','','ens.max','', |
|
241 |
'depth.min','','elapsed.min','','ens.min','', |
|
43 | 242 |
'out_basename','','outgrid_dz','','run_label','', |
243 |
'outgrid_firstbin','','outgrid_lastbin','', |
|
244 |
'outgrid_minsamp','','per_bin_valid_frac_lim','', |
|
245 |
'processing_options','','refLr_firstBin','', |
|
246 |
'refLr_lastBin','','rms_w_reflr_err','', |
|
247 |
'rms_w_reflr_err_interior','', |
|
248 |
'sidelobe_editing','','surface_layer_depth','', |
|
56 | 249 |
'vessel_draft','','w_lim.max','', |
43 | 250 |
'water_depth.sig','','water_depth_from','', |
251 |
); |
|
252 |
} |
|
253 |
close(TOCLD); # close LADCP_VKE input |
|
254 |
my(@specrec); |
|
255 |
while (@specrec = &antsFileIn(FROMCLD)) { |
|
256 |
my($i) = ($specrec[$widxf]>0) ? $specrec[$widxf] : 0; |
|
257 |
@{$specbuf[$bufi][$i]} = @specrec; |
|
44 | 258 |
# print(STDERR "specrec[$i] = @specrec\n"); |
43 | 259 |
} |
260 |
close(FROMCLD); |
|
44 | 261 |
# printf(STDERR "windows up to %d buffered\n",$#{$specbuf[$bufi]}); |
43 | 262 |
|
263 |
$antsBufSkip = @ants_; |
|
264 |
&antsIn(); # read next .wprof file |
|
265 |
open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # process it |
|
266 |
croak("LADCP_wspec $opts: $!\n"); |
|
267 |
print(TOCLD $antsOldHeaders); |
|
268 |
for (my($r)=0; $r<@ants_; $r++) { |
|
269 |
print(TOCLD "@{$ants_[$r]}\n"); |
|
270 |
} |
|
271 |
} |
|
272 |
||
273 |
close(TOCLD); # connect stdout from LADCP_VKE to stdin |
|
29 | 274 |
open(STDIN,"<&FROMCLD") || croak("dup(FROMCLD): $!\n"); |
275 |
close(FROMCLD); |
|
43 | 276 |
<>; # clear EOF condition |
277 |
||
42 | 278 |
undef(%P); # shouldn't matter, because we'll get the same %PARAMs back |
279 |
undef(@antsLayout); # shouldn't matter, because it will get overwritten |
|
43 | 280 |
undef($antsOldHeaders); # forget those |
281 |
undef(@ants_); |
|
49 | 282 |
$antsBufNFields = 0; |
43 | 283 |
|
284 |
} elsif (defined(fnrNoErr('pwrdens.0'))) { |
|
285 |
croak("$0: -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n") |
|
286 |
if ($opt_d || $opt_u || $opt_b || defined($opt_w) || defined($opt_s) || defined($opt_g)); |
|
29 | 287 |
} else { |
288 |
if ($ARGV[0]) { |
|
289 |
croak("$ARGV[0]: no such file or directory\n"); |
|
290 |
} else { |
|
291 |
croak("$0: empyt input\n"); |
|
292 |
} |
|
293 |
} |
|
18 | 294 |
|
29 | 295 |
#---------------------------------------------------------------------- |
43 | 296 |
# Handle LADCP_VKE usage & read spectra from final file |
297 |
# - spectra from previous files are in @specbuf |
|
29 | 298 |
#---------------------------------------------------------------------- |
299 |
||
43 | 300 |
$n_input_files = 1 + @specbuf; # number of input files provided |
301 |
||
56 | 302 |
&antsFloatOpt(\$opt_e,$c); # default parameterization |
303 |
&antsFloatOpt(\$opt_x,1); # spectral fit stddev scale factor |
|
304 |
||
43 | 305 |
&antsAddParams('LADCP_VKE::input_files.n',$n_input_files, |
48 | 306 |
'LADCP_VKE::wsamp.min',$opt_z, |
56 | 307 |
'LADCP_VKE::eps.minlim',$opt_l, |
308 |
'LADCP_VKE::calibration_constant',$opt_e); |
|
18 | 309 |
|
42 | 310 |
if (defined($opt_c)) { # shortwave cutoff supplied |
29 | 311 |
$lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c; |
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
312 |
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below |
42 | 313 |
} elsif (defined(antsParam('shortwave_cutoff'))) { # cutoff already applied |
29 | 314 |
$lmin = 2*$PI/antsParam('shortwave_cutoff'); |
42 | 315 |
} else { # use 100m default cutoff |
43 | 316 |
$lmin = $opt_c_default; |
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
317 |
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below |
29 | 318 |
} |
42 | 319 |
$lmax = 9e99; # no longwave cutoff implemented yet |
29 | 320 |
|
42 | 321 |
&antsInstallBufFull(0); # load entire file |
29 | 322 |
&antsIn(); |
43 | 323 |
$P{run_label} = "$input_list$P{profile_id}($P{run_label})" |
324 |
if ($n_input_files > 1); |
|
325 |
||
326 |
my($Hbuf) = $antsOldHeaders; # save for later (used on -i) |
|
29 | 327 |
|
328 |
&antsRequireParam('profile_id'); |
|
329 |
&antsRequireParam('lambda.1'); |
|
330 |
&antsRequireParam('k.1'); |
|
331 |
&antsRequireParam('resolution_bandwidth'); |
|
332 |
&antsRequireParam('input_depth_resolution'); |
|
333 |
&antsRequireParam('output_depth_resolution'); |
|
334 |
&antsRequireParam('ADCP_bin_length'); |
|
335 |
&antsAddParams('ADCP_pulse_length',antsParam('ADCP_bin_length')) |
|
336 |
unless defined(antsParam('ADCP_pulse_length')); |
|
337 |
||
42 | 338 |
$imin = 0; # find frequency bin limits |
29 | 339 |
for ($nfreq=1; defined(antsParam("lambda.$nfreq")); $nfreq++) { |
340 |
$imin = $nfreq if ($imin==0 && antsParam("lambda.$nfreq")<=$lmax); |
|
341 |
$imax = $nfreq if (antsParam("lambda.$nfreq") >= $lmin); |
|
18 | 342 |
} |
343 |
croak("$0: <lambda.min=$lmin> < min(lambda)") |
|
344 |
unless defined($imax); |
|
345 |
||
346 |
$pg_fmin = fnr('pwrdens.0'); # first power field in spectra |
|
347 |
$fs_fmin = $pg_fmin + $imin; # first power field in finescale range |
|
348 |
$fs_fmax = $pg_fmin + $imax; # last power field in finescale range |
|
349 |
||
29 | 350 |
$widxf = fnr('widx'); |
43 | 351 |
$df = fnr('depth'); |
18 | 352 |
$mindf = fnr('depth.min'); |
353 |
$maxdf = fnr('depth.max'); |
|
42 | 354 |
$wsf = fnr('w.nsamp.avg'); |
43 | 355 |
$doff = fnr('nspec'); |
18 | 356 |
|
357 |
#---------------------------------------------------------------------- |
|
29 | 358 |
# Redirect STDOUT & create plot if STDOUT is a tty |
18 | 359 |
#---------------------------------------------------------------------- |
360 |
||
29 | 361 |
if (-t STDOUT) { |
33 | 362 |
$opt_f = '.' unless defined($opt_f); # output directory |
363 |
croak("$opt_f: not a directory\n") unless (-d $opt_f); |
|
364 |
||
29 | 365 |
my($id) = &antsRequireParam('profile_id'); |
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
366 |
|
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
367 |
if (defined($basename)) { |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
368 |
$opt_p = sprintf('%s/%s_VKE.ps',$opt_f,$basename) unless defined($opt_p); |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
369 |
$outfile = sprintf('%s/%s.VKE',$opt_f,$basename); |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
370 |
} else { |
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
371 |
$opt_p = sprintf('%s/%03d_VKE.ps',$opt_f,$id) unless defined($opt_p); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
372 |
$outfile = sprintf('%s/%03d.VKE',$opt_f,$id); |
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
373 |
} |
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
374 |
$outfile =~ s@^./@@; |
29 | 375 |
open(STDOUT,">$outfile") || die("$outfile: $!\n"); |
33 | 376 |
} elsif (defined($opt_f)) { |
377 |
croak("-f can only be used without STDOUT redirection\n"); |
|
29 | 378 |
} |
379 |
||
380 |
#---------------------------------------------------------------------- |
|
381 |
# Library |
|
382 |
#---------------------------------------------------------------------- |
|
383 |
||
43 | 384 |
sub average_spectra($) |
385 |
{ |
|
386 |
my($r) = @_; |
|
387 |
for (my($f)=$fs_fmin-1; $f<=$fs_fmax; $f++) { # average w.nsamp.avg & spectral densities |
|
388 |
my($sum) = my($ns) = 0; |
|
389 |
if (numberp($ants_[$r][$f])) { # final file has valid spectra |
|
390 |
$sum = $ants_[$r][$f]; $ns = 1; # NB: nspec is correct even if it doesn't |
|
391 |
} |
|
392 |
my($wi) = ($ants_[$r][$widxf]>0) ? $ants_[$r][$widxf] : 0; # adjust for -1 |
|
393 |
for (my($bi)=0; $bi<@specbuf; $bi++) { # loop over all buffered files |
|
394 |
next unless @{$specbuf[$bi][$wi]}; # skip input files w/o valid spectra |
|
395 |
if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) > 0) { # depth mismatch |
|
56 | 396 |
# die("assertion failed ($specbuf[$bi][$wi][$df] - $ants_[$r][$df] @ wi=$wi)") unless ($wi == 0); # only allowed in bottom window |
43 | 397 |
if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) > |
398 |
abs($ants_[$r][$maxdf] - $ants_[$r][$mindf])) { |
|
56 | 399 |
printf(STDERR "WARNING: ignoring window #$wi from input file #%d because of depth mismatch\n",$bi+1) |
43 | 400 |
if ($f == $fs_fmin); |
401 |
next; |
|
402 |
} |
|
403 |
} |
|
44 | 404 |
if (numberp($specbuf[$bi][$wi][$f])) { |
405 |
$sum += $specbuf[$bi][$wi][$f]; $ns++; |
|
406 |
} |
|
43 | 407 |
$ants_[$r][$doff] += $specbuf[$bi][$wi][$doff] # update nspec once per input record |
408 |
if ($f == $fs_fmax); # ... but for all files |
|
409 |
} |
|
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
410 |
$ants_[$r][$f] = ($ants_[$r][$doff] >= $opt_h) ? $sum/$ns : nan; # update averaged spectral density |
43 | 411 |
} |
412 |
} |
|
413 |
||
414 |
sub integrate_fs_power($) # integrate fs spectrum |
|
18 | 415 |
{ |
416 |
my($r) = @_; |
|
417 |
||
418 |
$ants_[$r][$fspwrf] = 0; |
|
419 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
|
29 | 420 |
$ants_[$r][$fspwrf] += $ants_[$r][$f]; |
18 | 421 |
} |
29 | 422 |
$ants_[$r][$fspwrf] *= antsParam('resolution_bandwidth'); |
18 | 423 |
} |
424 |
||
425 |
||
42 | 426 |
sub fit_universal_w_spec($) # vertical velocity => p0 |
18 | 427 |
{ |
428 |
my($r) = @_; |
|
429 |
my($nsamp) = $fs_fmax - $fs_fmin + 1; |
|
430 |
||
431 |
#--------------------------------------------------- |
|
432 |
# fit slope-2 line in log-log space (main estimator) |
|
433 |
#--------------------------------------------------- |
|
434 |
||
42 | 435 |
if ($nsamp >= 2) { # require min 2 wavenumber samples |
18 | 436 |
|
43 | 437 |
my($DOF) = 0; |
438 |
||
42 | 439 |
my($sumd,$sumx,$sumy) = (0,0,0); # fit kz^-2 power law |
440 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
|
18 | 441 |
my($i) = $f - $pg_fmin; |
29 | 442 |
$sumx += log10(antsParam("k.$i")); |
18 | 443 |
$sumy += log10($ants_[$r][$f]); |
43 | 444 |
$sumd += log10($ants_[$r][$f]) + 2*log10(antsRequireParam("k.$i")); |
18 | 445 |
} |
446 |
my($p0) = $sumd/$nsamp; |
|
447 |
$ants_[$r][$p0f] = 10**$p0; |
|
448 |
||
42 | 449 |
my($avgx) = $sumx/$nsamp; # avg for r calc |
18 | 450 |
my($avgy) = $sumy/$nsamp; |
42 | 451 |
my($sumsqerr,$sxx,$syy,$sxy,$sumsqxt,$sumx,$sumy,$sumwt) = # r, rms error, pwrlaw slope |
452 |
(0,0,0,0,0,0,0,0); |
|
18 | 453 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
42 | 454 |
my($i) = $f - $pg_fmin; |
455 |
my($x) = log10(&antsParam("k.$i")); |
|
456 |
my($y) = log10($ants_[$r][$f]); |
|
43 | 457 |
my($ysig) = $opt_x * $y / sqrt($ants_[$r][$doff]); |
42 | 458 |
my($xt) = $x - $avgx; $sxx += &SQR($xt); # correlation coeff (r) |
459 |
my($yt) = $y - $avgy; $syy += &SQR($yt); $sxy += $xt * $yt; |
|
460 |
my($wt) = 1 / &SQR($ysig); $sumwt += $wt; # slope (linear fit in log-log space) |
|
461 |
$sumx += $x * $wt; |
|
462 |
$sumy += $y * $wt; |
|
463 |
$sumsqerr += ($p0 - 2*$x - $y)**2; # rms error |
|
18 | 464 |
} |
42 | 465 |
my($midx) = $sumx / $sumwt; |
466 |
# print(STDERR "$sumx:$sumy:$sumwt\n"); |
|
467 |
# print(STDERR "$midx\n"); |
|
468 |
my($sumsqdx,$sumslp) = (0,0); |
|
469 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
|
470 |
my($i) = $f - $pg_fmin; |
|
471 |
my($x) = log10(&antsParam("k.$i")); |
|
472 |
my($y) = log10($ants_[$r][$f]); |
|
43 | 473 |
my($ysig) = $opt_x * $y / sqrt($ants_[$r][$doff]); |
42 | 474 |
my($dx) = ($x - $midx) / $ysig; |
475 |
$sumsqdx += &SQR($dx); |
|
476 |
$sumslp += $dx * $y / $ysig; |
|
477 |
} |
|
478 |
$ants_[$r][$rmsf] = sqrt($sumsqerr/$nsamp); |
|
479 |
$ants_[$r][$rf] = $sxy/(sqrt($sxx * $syy) + $SMALL_AMOUNT); |
|
480 |
$ants_[$r][$slpf] = $sumslp / $sumsqdx; |
|
481 |
$ants_[$r][$sslpf] = sqrt(1 / $sumsqdx); |
|
18 | 482 |
|
483 |
} else { |
|
484 |
&antsInfo("WARNING: no fit --- need min 2 samples"); |
|
485 |
} |
|
42 | 486 |
} |
487 |
||
488 |
#---------------------------------------------------------------------- |
|
489 |
# Load Dissipation Data |
|
490 |
#---------------------------------------------------------------------- |
|
491 |
||
492 |
my(@eps_ms,@depth_ms); # output variables |
|
493 |
if (defined($opt_k)) { |
|
494 |
my($file,$field) = split(':',$opt_k); |
|
49 | 495 |
$field = 'eps' unless defined($field); |
42 | 496 |
open(ADDF,"$file") || croak("$file: $!\n"); # open file |
497 |
my(@afl) = &antsFileLayout(ADDF); # read layout |
|
498 |
my($akf,$aef); |
|
499 |
for (my($f)=0; $f<=$#afl; $f++) { # find depth & eps fields |
|
500 |
$akf = $f if ($afl[$f] eq 'depth'); |
|
501 |
$aef = $f if ($afl[$f] eq $field); |
|
502 |
} |
|
503 |
croak("$file: fields 'depth' and '$field' required\n") |
|
504 |
unless defined($akf) && defined($aef); |
|
505 |
while (1) { # load entire profile |
|
506 |
my(@ar); |
|
507 |
last unless @ar = &antsFileIn(ADDF); |
|
508 |
next unless numberp($ar[$akf]) && numberp($ar[$aef]); |
|
509 |
push(@depth_ms,$ar[$akf]); push(@eps_ms,$ar[$aef]); |
|
510 |
} |
|
511 |
close(ADDF); |
|
18 | 512 |
} |
513 |
||
514 |
#---------------------------------------------------------------------- |
|
515 |
# Process File |
|
516 |
#---------------------------------------------------------------------- |
|
517 |
||
43 | 518 |
$fspwrf = &antsNewField('pwr.fs'); # derived fields |
18 | 519 |
|
42 | 520 |
$p0f = &antsNewField('p0'); # VKE density |
521 |
$rmsf = &antsNewField('p0fit.rms'); # rms misfit |
|
522 |
$rf = &antsNewField('p0fit.r'); # correlation coefficient |
|
523 |
$slpf = &antsNewField('p0fit.slope'); # power-law slope |
|
524 |
$sslpf = &antsNewField('p0fit.slope.sig'); # power-law slope stddev |
|
525 |
$wepsf = &antsNewField('eps.VKE'); # epsilon from VKE |
|
49 | 526 |
$msepsf = &antsNewField('eps.ms'); # externally supplied microstructure eps if available |
18 | 527 |
|
528 |
my(@outLayout) = @antsNewLayout; # save for later |
|
529 |
for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input |
|
530 |
$totf = $f if ($outLayout[$f] eq 'pwr.tot'); |
|
531 |
$tnsf = $f if ($outLayout[$f] eq 'pwr.tot.nsamp'); |
|
532 |
} |
|
533 |
croak("$0: cannot find fields 'pwr.tot' or 'pwr.tot.nsamp' in input\n") |
|
534 |
unless defined($totf) || defined($tnsf); |
|
29 | 535 |
$lsf = defined($tnsf) ? $tnsf : $totf; |
18 | 536 |
|
29 | 537 |
if ($opt_p) { # begin plot |
538 |
$xmin = (&antsParam('output_depth_resolution')>=350) ? 0.012 : 0.018; |
|
539 |
$xmax = 0.2; $ymin = 1; $ymax = 1e4; |
|
540 |
$plotsize = 13; |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
541 |
GMT_begin($opt_p,"-JX${plotsize}l","-R$xmin/$xmax/$ymin/$ymax",'-P -X6 -Y4'); # init plot |
29 | 542 |
GMT_psxy('-G25/255/25 -L'); |
543 |
printf(GMT "%f %f\n",$xmin,$xmin**(-2)*sqrt(2)); |
|
544 |
printf(GMT "%f %f\n",$xmax,$xmax**(-2)*sqrt(2)); |
|
545 |
printf(GMT "%f %f\n",$xmax,$xmax**(-2)/sqrt(2)); |
|
546 |
printf(GMT "%f %f\n",$xmin,$xmin**(-2)/sqrt(2)); |
|
547 |
} |
|
548 |
||
549 |
my(@sumPwr,@nPwr); |
|
550 |
my(@sumGoodPwr,@nGoodPwr); |
|
551 |
my($min_depth) = 9e99; |
|
552 |
my($max_depth) = -9e99; |
|
18 | 553 |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
554 |
my($latM) = abs(&antsRequireParam('lat')); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
555 |
&antsInfo("WARNING: low latitude-profile no epsilon estimated") |
43 | 556 |
unless ($latM > $opt_q); |
18 | 557 |
|
43 | 558 |
for (my($r)=0; $r<@ants_; $r++) { # loop over all windows |
559 |
||
560 |
if (numberp($ants_[$r][$pg_fmin])) { # there is a spectrum |
|
18 | 561 |
|
43 | 562 |
#-------------------------- |
563 |
# apply spectral correction |
|
564 |
#-------------------------- |
|
565 |
||
566 |
unless ($opt_m) { |
|
567 |
for (my($i)=0; $i<$nfreq; $i++) { # loop over wavenumbers |
|
568 |
$ants_[$r][$i+$pg_fmin] *= |
|
569 |
T_w(antsParam("k.$i"),antsParam('ADCP_bin_length'), |
|
570 |
antsParam('ADCP_pulse_length'),antsParam('input_depth_resolution'), |
|
571 |
$opt_r); |
|
572 |
} |
|
573 |
} |
|
574 |
||
575 |
#------------------------ |
|
576 |
# calculate fs quantities |
|
577 |
#------------------------ |
|
578 |
||
579 |
average_spectra($r); # average all avaiable spectra |
|
580 |
integrate_fs_power($r); # calculate total finescale power |
|
581 |
fit_universal_w_spec($r); # fit kz^-2 spectrum & calc stats |
|
582 |
||
583 |
if (numberp($ants_[$r][$p0f])) { # update min/max depth |
|
584 |
$min_depth = $ants_[$r][$mindf] if ($ants_[$r][$mindf] < $min_depth); |
|
585 |
$max_depth = $ants_[$r][$maxdf] if ($ants_[$r][$maxdf] > $max_depth); |
|
18 | 586 |
} |
587 |
||
43 | 588 |
} else { # no spectrum |
589 |
$ants_[$r][$fspwrf] = $ants_[$r][$p0f] = $ants_[$r][$rmsf] = |
|
590 |
$ants_[$r][$rf] = $ants_[$r][$slpf] = $ants_[$r][$sslpf] = nan; |
|
42 | 591 |
} |
592 |
||
593 |
#----------------------------------------------------------------------------------------------------- |
|
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
594 |
# Except for -z and the low-latitude test, all others were disabled in October 2021 for version 2.1. |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
595 |
# This decision was based on two main reasons: 1) The checks were not used for the estimates |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
596 |
# in the Lele paper, which is the first major publiction involving the parameterization. |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
597 |
# 2) While working on the 2021 QC of the GO-SHIP A20 section, which passes through a region |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
598 |
# of extremlely weakbackscatter, I found that regional mean profiles of the filtered |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
599 |
# eps.VKE estimates agree very well with (p0/0.022)**2, but the latter are less noisy |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
600 |
# because there are signficantly more smaples. Importantly, the filtered eps.VKE were |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
601 |
# derived with c=0.026, so as an imporant side effect of disabling the filters I returned |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
602 |
# to the c value from the original publication. |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
603 |
# |
56 | 604 |
# eps.VKE QC Tests: |
42 | 605 |
# - the following limits were independently derived |
606 |
# p0fit.rms <= 0.4 primary filter used in Thurnherr et al. (GRL 2015) |
|
607 |
# -3 <= p0fit.slope <= -1 based largely on 2016 I08S data with sufficient/insufficient range |
|
608 |
# p0fit.r <= -0.5 based largely on 2016 I08S data with sufficient/insufficient range |
|
609 |
# w.nsamp.avg >= 50 based on observations in many data sets with weak backscatter, |
|
610 |
# including DoMORE, GO-SHIP P16S, GO-SHIP I08S |
|
611 |
# - then, I plotted slope & r vs. rms and found that rms = 0.4 corresponds to, on average |
|
612 |
# -3 <= p0fit.slope <= -1 |
|
613 |
# p0fit.r <= -0.4 |
|
614 |
# - in a plot of rms vs nsamp the limiting value of 0.4 is hit at 50 samples |
|
615 |
# |
|
54 | 616 |
# => SET OF MUTUALLY CONSISTENT CRITERIA |
42 | 617 |
# |
618 |
# Additional Empirical Filters: |
|
43 | 619 |
# - latitude > 3deg guess based on Thurnherr et al., 2015, Gregg et al., 2003, 2015 GO-SHIP P16N |
54 | 620 |
# - eps >= 5e-11 W/kg based on DIMES data, cutting where errors become > factor 2; DISABLED SEPTEMBER 2019 |
42 | 621 |
#----------------------------------------------------------------------------------------------------- |
622 |
||
48 | 623 |
if ($latM > $opt_q && # 1) not (too) equatorial |
57
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
624 |
# $ants_[$r][$rmsf] <= 0.4 && # 2) rms spectra misfit <= 0.4 (as in Thurnherr et al., GRL 2015) |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
625 |
# $ants_[$r][$slpf]>=-3 && $ants_[$r][$slpf]<=-1 && # 3) slope consistent with -2 power law |
69e39fcb7f41
whoosher version pre P2
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
56
diff
changeset
|
626 |
# $ants_[$r][$rf] <= -0.4 && # 4) p and k_z are well correlated |
48 | 627 |
$ants_[$r][$wsf] >= $opt_z) { # 5) minimum # of samples |
628 |
if (($ants_[$r][$p0f]/$opt_e)**2 >= $opt_l) { # level is above eps.minlim (-l) |
|
629 |
$ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2; # => Thurnherr et al. (GRL 2015) |
|
630 |
} else { # level is below eps.minlim |
|
631 |
$ants_[$r][$wepsf] = $opt_a; # => set to arbitrary background value |
|
632 |
} |
|
633 |
} else { # failed any of checks 1-5 => nan |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
634 |
$ants_[$r][$wepsf] = nan; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
635 |
} |
29 | 636 |
|
42 | 637 |
#------------------------------------------------- |
638 |
# average external microstructure eps, if supplied |
|
639 |
#------------------------------------------------- |
|
640 |
||
641 |
if (defined($opt_k)) { |
|
642 |
my($sum,$n) = (0,0); |
|
643 |
for (my($i)=0; $i<@eps_ms; $i++) { # linearly search all eps records |
|
644 |
next unless $depth_ms[$i] >= $ants_[$r][$mindf] && |
|
645 |
$depth_ms[$i] <= $ants_[$r][$maxdf]; |
|
646 |
$sum += $eps_ms[$i]; $n++; |
|
647 |
} |
|
648 |
$ants_[$r][$msepsf] = $n ? $sum / $n : nan; |
|
49 | 649 |
} else { |
650 |
$ants_[$r][$msepsf] = nan; |
|
651 |
} |
|
18 | 652 |
|
653 |
#--------------- |
|
654 |
# produce output |
|
655 |
#--------------- |
|
656 |
||
29 | 657 |
if (defined($opt_p)) { # plot current spectrum on -p & calc mean |
658 |
my($R) = 0; # RGB map |
|
659 |
my($G) = int(200*(1-$r/@ants_)); |
|
660 |
my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_)); |
|
41 | 661 |
GMT_psxy("-W1,$R/$G/$B"); |
29 | 662 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { # avg & plot high-quality spectra only |
663 |
next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]); # - omit zero wavenumber ($pg_fmin) |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
664 |
next unless numberp($ants_[$r][$wepsf]); |
29 | 665 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
666 |
printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]); |
|
667 |
$nGoodPwr[$f]++; |
|
668 |
$sumGoodPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f]; |
|
669 |
} |
|
41 | 670 |
GMT_psxy("-W1,$R/$G/$B,."); # avg & plot all spectra with dots |
29 | 671 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { |
672 |
next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]); |
|
673 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
674 |
printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]); |
|
675 |
$nPwr[$f]++; |
|
676 |
$sumPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f]; |
|
677 |
} |
|
678 |
} |
|
679 |
||
680 |
||
681 |
if (defined($opt_i)) { # output current spectrum on -i |
|
18 | 682 |
open(STDOUT_DUP,">&",STDOUT) || croak("$0: cannot dup STDOUT\n"); |
683 |
@antsNewLayout = ('k','lambda','pwrdens','finescale','pwrdens_fit'); |
|
29 | 684 |
if ($ants_[$r][$widxf] > 0) { |
685 |
$ofname = sprintf("${opt_i}_%02d.wspec",$r+1); |
|
686 |
} else { |
|
687 |
$ofname = sprintf("${opt_i}_btm.wspec"); |
|
688 |
} |
|
18 | 689 |
open(STDOUT,">$ofname") || croak("$ofname: $!\n"); |
690 |
undef($antsOldHeaders); |
|
26
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
691 |
&antsActivateOut() if ($ANTS_TOOLS_AVAILABLE); |
18 | 692 |
|
693 |
my($saveParams) = $antsCurParams; # add all extra input fields as %PARAMs |
|
43 | 694 |
&antsAddParams('LADCP_VKE::depth.min',$ants_[$r][$mindf], |
695 |
'LADCP_VKE::depth.max',$ants_[$r][$maxdf]); |
|
18 | 696 |
for (my($f)=$lsf+1; $f<@outLayout; $f++) { |
43 | 697 |
&antsAddParams("LADCP_VKE::$outLayout[$f]",$ants_[$r][$f]); |
18 | 698 |
} |
699 |
||
29 | 700 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { |
701 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
702 |
my($l) = antsParam(sprintf("lambda.%d",$f-$pg_fmin)); |
|
18 | 703 |
&antsOut($k,$l,$ants_[$r][$f], |
704 |
($f>=$fs_fmin && $f<=$fs_fmax), |
|
705 |
numberp($ants_[$r][$p0f]) ? $ants_[$r][$p0f] * $k**(-2) : nan); |
|
706 |
} |
|
707 |
||
708 |
&antsOut('EOF'); |
|
709 |
open(STDOUT,">&",STDOUT_DUP) || croak("$0: cannot restore STDOUT\n"); |
|
710 |
close(STDOUT_DUP); |
|
711 |
$antsCurParams = $saveParams; |
|
712 |
} |
|
713 |
} |
|
714 |
||
29 | 715 |
if (defined($opt_p)) { # complete plot |
716 |
if (@nPwr) { |
|
41 | 717 |
GMT_psxy('-W2,255/100/100,.'); # plot mean spectrum; dotted; entire range |
29 | 718 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
719 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
720 |
printf(GMT "$k %g\n",$sumPwr[$f]/$nPwr[$f]); |
|
721 |
} |
|
722 |
} |
|
723 |
||
724 |
if (@nGoodPwr) { |
|
41 | 725 |
GMT_psxy('-W1,255/100/100'); # plot mean spectrum; entire range |
29 | 726 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { |
727 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
728 |
printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]); |
|
729 |
} |
|
41 | 730 |
GMT_psxy('-W4,255/100/100'); # plot mean fs spectrum (heavy) |
29 | 731 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
732 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
733 |
printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]); |
|
734 |
} |
|
735 |
} |
|
736 |
||
737 |
GMT_psbasemap('-Bf3a2:"Vertical Wavenumber [rad/m]":/' . # annotate axes |
|
738 |
'f3a1:"VKE/p@-0@- [m@+3@+s@+-1/2@+/(rad m@+-1@+)]":WeS'); |
|
739 |
GMT_setJ("-JX-${plotsize}l"); |
|
740 |
GMT_setR(sprintf('-R%f/%f/%f/%f',2*$PI/$xmax,2*$PI/$xmin,$ymin,$ymax)); |
|
741 |
GMT_psbasemap('-B3:"Vertical Wavelength [m]:N"'); |
|
742 |
||
30 | 743 |
GMT_unitcoords_logscale(); # print profile number |
43 | 744 |
GMT_pstext('-F+f14,Helvetica,blue+jTL -Gwhite'); |
745 |
if (defined($outfile)) { printf(GMT "0.02 0.98 $outfile [$P{run_label}] %s\n",antsParam('input_data')); } |
|
746 |
else { printf(GMT "0.02 0.98 %03d [$P{run_label}] %s\n",antsParam('profile_id'),antsParam('input_data')); } |
|
41 | 747 |
GMT_pstext('-F+f9,Helvetica,orange+jTR -N -Gwhite'); |
748 |
print(GMT "0.99 0.99 V$VERSION\n"); |
|
29 | 749 |
|
41 | 750 |
GMT_set('FONT_ANNOT_PRIMARY 10','MAP_ANNOT_OFFSET_PRIMARY 0.01c', # eps profile inset |
751 |
'FONT_LABEL 10','MAP_LABEL_OFFSET 0.2c'); |
|
29 | 752 |
$min_depth = 400 if ($min_depth > 1e4); |
753 |
$max_depth = 1500 if ($max_depth < 0); |
|
42 | 754 |
GMT_setR(sprintf("-R1e-11/1e-7/%d/%d",round($min_depth-250,500),round($max_depth+250,500))); |
29 | 755 |
GMT_setJ(sprintf('-JX%fl/-%f',$plotsize/2.2,$plotsize/2.2)); |
756 |
GMT_psbasemap('-X0.3 -Y0.3 -Bg1a1f3p:"@~e@~@-VKE@- [W kg@+-1@+]":/g1000a500f100:"Depth [m]":wEsN'); |
|
41 | 757 |
GMT_psxy(); |
29 | 758 |
for (my($r)=0; $r<@ants_; $r++) { |
42 | 759 |
my($R) = 0; my($G) = int(200*(1-$r/@ants_)); # calculate color ramp |
29 | 760 |
my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_)); |
42 | 761 |
if (numberp($ants_[$r][$msepsf])) { # microstructure eps |
762 |
print(GMT "> -W1.5,DarkOrange2\n"); |
|
763 |
print(GMT "$ants_[$r][$msepsf] $ants_[$r][$mindf]\n"); |
|
764 |
print(GMT "$ants_[$r][$msepsf] $ants_[$r][$maxdf]\n"); |
|
765 |
} |
|
766 |
if (numberp($ants_[$r][$wepsf])) { |
|
767 |
print(GMT "> -W2,$R/$G/$B\n"); # plot eps.w in blue |
|
768 |
print(GMT "$ants_[$r][$wepsf] $ants_[$r][$mindf]\n"); |
|
769 |
print(GMT "$ants_[$r][$wepsf] $ants_[$r][$maxdf]\n"); |
|
770 |
} |
|
29 | 771 |
} |
772 |
GMT_end(); # finish plot |
|
773 |
} |
|
774 |
||
42 | 775 |
@antsNewLayout = @outLayout; # restore layout |
18 | 776 |
$antsOldHeaders = $Hbuf; |
777 |
$antsHeadersPrinted = 0; |
|
778 |
||
49 | 779 |
unless (defined($opt_y)) { # remove spectral bins from output |
780 |
splice(@antsNewLayout,$pg_fmin,$nfreq); |
|
781 |
for (my($r)=0; $r<@ants_; $r++) { |
|
782 |
splice(@{$ants_[$r]},$pg_fmin,$nfreq); |
|
783 |
} |
|
784 |
} |
|
785 |
||
786 |
&antsFlush(); # output results |
|
18 | 787 |
&antsExit(); |