#!/usr/bin/perl
#======================================================================
# L A D C P I N T S H
# doc: Thu Oct 14 21:22:50 2010
# dlm: Fri Mar 6 15:51:49 2015
# (c) 2010 A.M. Thurnherr & E. Firing
# uE-Info: 8 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
$antsSummary = 'integrate LADCP shear';
# NOTES:
# - the core of this code is a simplified version of avg_sh.m and
# int_sh.m written by Eric Firing
# - comments beginning with ## are taken from Eric's code
# - cubic velocity interpolation across PPI gap from Eric's code has
# not been implemented (yet?)
# - low-pass-filtered shear code has not yet been implemented
# - currently, shear gaps are assumed to have vanishing shear;
# better solutions are possible
# - elapsed time is simply copied from shear elapsed time (i.e.
# it is not interpolated onto the new depth grid)
# WEIRDNESSES IN Eric's CODE:
# - in Eric's [avg_sh.m] the calculation of output shear stddev incorrectly assumes that the 4th column
# in the shear profile is stddev, rather than variance. However, as far as I can tell, this output
# is not used anywhere in Eric's code
# HISTORY:
# Oct 14, 2010: - created
# Oct 20, 2010: - first working version
# Oct 23, 2010: - added support for -b)
# Oct 24, 2010: - fix spuriously small variances that can occur for BT velocities based on very small
# samples (i.e. primarily when chosing a small -r)
# Dec 9, 2010: - allowed for empty BT file
# Feb 16, 2011: - BUG: gaps in shear data were not handled correctly in baroclinic solution
# Jul 7, 2011: - added -m
# Jul 19, 2011: - added shear sigma output for shear inversion
# - BUG: dc/uc v component was not correctly referenced
# Jul 24, 2011: - BUG: BT constraint was erroneously assumed to be available for dc,uc only, if it
# was available for <dc,uc> combo profile
# Jul 25, 2011: - BUG: nan in either of ul/dl u_z.sig caused combined u_z.sig to be nan, too
# Jul 27, 2011: - removed code related to smoothed shear
# - replaced -r by -n
# - removed shear sigma output
# - replaced ndata by nsamp
# - removed -w (Eric's way of dealing with dc/uc temporal variability)
# Feb 19, 2012: - added processing of elapsed time
# - adapted to new shear file layout (nsamp instead of nshear)
# - BUG: uplooker data was not used for downcasts and and only partially for combo data
# May 16, 2012: - adapted to ANTSlib V5.0
# May 25, 2012: - added code to read LDEO_IX bottom-track data
# Jun 14, 2012: - noticed that -b now works also with LDEO SADCP files :-); renamed -b option to -r
# Jun 5, 2013: - BUG: code bombed when either UC or DC was missing
# Jun 28, 2013: - adapated to new :: convention
# - make sure LADCP DUL metadata are dealt with correctly
# Jul 12, 2013: - clarified -u usage with better messages
# Mar 20, 2014: - fiddled while debugging [LADCPproc]
# Jun 7, 2014: - improved error messages
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
require "$ANTS/ants.pl";
require "$ANTS/libstats.pl";
&antsUsage('dm:n:r:s:w:u:',0,
'[-d)ebug]',
'[-r)eference with <BT or SADCP file> [-m)in <samp[10]>]]',
'[secondary -u)plooker <shear file>]',
'[min -n) <shear samp[10]>]',
'[output -s)hear-pro <file>]',
'[LADCP shear file]');
croak("$0: -m meaningless without -r\n")
if defined($opt_m) && !defined($opt_r);
&antsCardOpt(\$opt_n,10); # minimum number of samples for shear
$minBTsamp = &antsCardOpt($opt_m,10); # minimum number of samples for BT data
&antsFileOpt($opt_r); # reference velocity file
&antsFileOpt($opt_u); # UL shear file
if (defined($opt_u)) {
open(ULF,$opt_u) || croak("$opt_u: $!\n");
%UL_P = &antsFileParams(ULF);
}
#======================================================================
# Step 1: Read and Average Shear Data
# - depth bins with less than $opt_n values are blanked out
#======================================================================
sub wavg_sig(@)
{
my($sumSq) = my($n) = 0;
for (my($i)=0; $i<$#_; $i+=2) {
next unless numberp($_[$i+1]);
$sumSq += $_[$i] * $_[$i+1]**2;
$n += $_[$i];
}
return ($n>0) ? sqrt($sumSq/$n) : nan;
}
#--------------------
# Handle Metadata
#--------------------
if (%UL_P) {
croak("$0: inconsistent vertical resolution\n")
unless ($P{LADCPproc::vertical_resolution} == $UL_P{LADCPproc::vertical_resolution});
unless ($P{LADCPproc::bin_length} == $UL_P{LADCPproc::bin_length}) {
&antsInfo("Warnining: different DL/UL bin lengths; derived spectra cannot be corrected");
&antsAddParams('LADCPproc::bin_length','',
'LADCPproc::DL_bin_length',$P{LADCPproc::bin_length},
'LADCPproc::UL_bin_length',$UL_P{LADCPproc::bin_length});
}
}
$depthF = fnrNoErr('depth'); # layout of [LADCPproc] output
unless (defined($depthF)) {
if (defined($opt_u)) {
croak("No 'depth' field in primary shear file (extraneous -u?)\n");
} else {
croak("No 'depth' field in primary shear file\n");
}
}
$dc_nshF = fnrNoErr('dc_nshear');
$dc_nshF = fnr('dc_nsamp') unless defined($dc_nshF);
$dc_uzF = fnr('dc_u_z');
$dc_uzsF = fnrNoErr('dc_u_z.sig');
$dc_uzsF = fnr('dc_u_z_sig') unless defined($dc_uzsF);
$dc_vzF = fnr('dc_v_z');
$dc_vzsF = fnrNoErr('dc_v_z.sig');
$dc_vzsF = fnr('dc_v_z_sig') unless defined($dc_vzsF);
$dc_wzF = fnr('dc_w_z');
$dc_wzsF = fnrNoErr('dc_w_z.sig');
$dc_wzsF = fnr('dc_w_z_sig') unless defined($dc_wzsF);
$dc_elapsedF = fnr('dc_elapsed');
$uc_nshF = fnrNoErr('uc_nshear');
$uc_nshF = fnr('uc_nsamp') unless defined($uc_nshF);
$uc_uzF = fnr('uc_u_z');
$uc_uzsF = fnrNoErr('uc_u_z.sig');
$uc_uzsF = fnr('uc_u_z_sig') unless defined($uc_uzsF);
$uc_vzF = fnr('uc_v_z');
$uc_vzsF = fnrNoErr('uc_v_z.sig');
$uc_vzsF = fnr('uc_v_z_sig') unless defined($uc_vzsF);
$uc_wzF = fnr('uc_w_z');
$uc_wzsF = fnrNoErr('uc_w_z.sig');
$uc_wzsF = fnr('uc_w_z_sig') unless defined($uc_wzsF);
$uc_elapsedF = fnr('uc_elapsed');
my(@gaps); my($curGap) = 0;
for (my($r)=0; &antsIn(); $r++) {
my(@UL_);
if (defined($opt_u)) {
@UL_ = &antsFileIn(ULF); # read UL shear data
undef($opt_u) unless (@UL_); # cheap trick
}
$depth[$r] = $ants_[0][$depthF]; ## depth grid values
croak("$opt_u: inconsistent depth record $r (DL: $depth[$r]; UL: $UL_[$depthF])\n")
if defined($opt_u) && ($UL_[$depthF] != $depth[$r]);
$dc_nsamp = $ants_[0][$dc_nshF]; # number of shear samples
$uc_nsamp = $ants_[0][$uc_nshF];
if (defined($opt_u)) {
$dl_nsamp = $dc_nsamp + $uc_nsamp;
$dc_nsamp += $UL_[$dc_nshF];
$uc_nsamp += $UL_[$uc_nshF];
$ul_nsamp = $dc_nsamp + $uc_nsamp;
}
$dc_nsamp[$r] = $dc_nsamp; # save for each record
$uc_nsamp[$r] = $uc_nsamp;
$nsamp[$r] = $dc_nsamp + $uc_nsamp;
if (defined($opt_u)) {
$ul_nsamp[$r] = $ul_nsamp;
$dl_nsamp[$r] = $dl_nsamp;
}
if (defined($opt_u)) { # dual-head instrument
if ($dc_nsamp > 0) { # downcast shear
my($DLf) = $ants_[0][$dc_nshF] / $dc_nsamp;
my($ULf) = $UL_[$dc_nshF] / $dc_nsamp;
if ($DLf>0 && $ULf>0) {
$dc_uz[$r] = $DLf*$ants_[0][$dc_uzF] + $ULf*$UL_[$dc_uzF];
$dc_vz[$r] = $DLf*$ants_[0][$dc_vzF] + $ULf*$UL_[$dc_vzF];
$dc_wz[$r] = $DLf*$ants_[0][$dc_wzF] + $ULf*$UL_[$dc_wzF];
$dc_elapsed[$r] = $DLf*$ants_[0][$dc_elapsedF] + $ULf*$UL_[$dc_elapsedF];
} elsif ($DLf > 0) {
$dc_uz[$r] = $ants_[0][$dc_uzF];
$dc_vz[$r] = $ants_[0][$dc_vzF];
$dc_wz[$r] = $ants_[0][$dc_wzF];
$dc_elapsed[$r] = $ants_[0][$dc_elapsedF];
} else {
$dc_uz[$r] = $UL_[$dc_uzF];
$dc_vz[$r] = $UL_[$dc_vzF];
$dc_wz[$r] = $UL_[$dc_wzF];
$dc_elapsed[$r] = $UL_[$dc_elapsedF];
}
} else {
$dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = $dc_elapsed[$r] = nan;
}
if ($uc_nsamp > 0) { # upcast shear
my($DLf) = $ants_[0][$uc_nshF] / $uc_nsamp;
my($ULf) = $UL_[$uc_nshF] / $uc_nsamp;
if ($DLf>0 && $Ulf>0) {
$uc_uz[$r] = $DLf*$ants_[0][$uc_uzF] + $ULf*$UL_[$uc_uzF];
$uc_vz[$r] = $DLf*$ants_[0][$uc_vzF] + $ULf*$UL_[$uc_vzF];
$uc_wz[$r] = $DLf*$ants_[0][$uc_wzF] + $ULf*$UL_[$uc_wzF];
$uc_elapsed[$r] = $DLf*$ants_[0][$uc_elapsedF] + $ULf*$UL_[$uc_elapsedF];
} elsif ($DLf > 0) {
$uc_uz[$r] = $ants_[0][$uc_uzF];
$uc_vz[$r] = $ants_[0][$uc_vzF];
$uc_wz[$r] = $ants_[0][$uc_wzF];
$uc_elapsed[$r] = $ants_[0][$uc_elapsedF];
} else {
$uc_uz[$r] = $UL_[$uc_uzF];
$uc_vz[$r] = $UL_[$uc_vzF];
$uc_wz[$r] = $UL_[$uc_wzF];
$uc_elapsed[$r] = $UL_[$uc_elapsedF];
}
} else {
$uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = $uc_elapsed[$r] = nan;
}
} else { # downlooker only
if ($dc_nsamp > 0) { # downcast shear
$dc_uz[$r] = $ants_[0][$dc_uzF];
$dc_vz[$r] = $ants_[0][$dc_vzF];
$dc_wz[$r] = $ants_[0][$dc_wzF];
$dc_elapsed[$r] = $ants_[0][$dc_elapsedF];
} else {
$dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = $dc_elapsed[$r] = nan;
}
if ($uc_nsamp > 0) { # upcast shear
$uc_uz[$r] = $ants_[0][$uc_uzF];
$uc_vz[$r] = $ants_[0][$uc_vzF];
$uc_wz[$r] = $ants_[0][$uc_wzF];
$uc_elapsed[$r] = $ants_[0][$uc_elapsedF];
} else {
$uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = $uc_elapsed[$r] = nan;
}
}
if ($nsamp[$r] > 0) {
my($dcf) = $dc_nsamp / $nsamp[$r];
my($ucf) = $uc_nsamp / $nsamp[$r];
if ($dcf>0 && $ucf>0) {
$uz[$r] = $dcf*$dc_uz[$r] + $ucf*$uc_uz[$r];
$vz[$r] = $dcf*$dc_vz[$r] + $ucf*$uc_vz[$r];
$wz[$r] = $dcf*$dc_wz[$r] + $ucf*$uc_wz[$r];
$elapsed[$r] = $dcf*$dc_elapsed[$r] + $ucf*$uc_elapsed[$r];
} elsif ($dcf > 0) {
$uz[$r] = $dc_uz[$r];
$vz[$r] = $dc_vz[$r];
$wz[$r] = $dc_wz[$r];
$elapsed[$r] = $dc_elapsed[$r];
} else {
$uz[$r] = $uc_uz[$r];
$vz[$r] = $uc_vz[$r];
$wz[$r] = $uc_wz[$r];
$elapsed[$r] = $uc_elapsed[$r];
}
} else {
$uz[$r] = $vz[$r] = $wz[$r] = $elapsed[$r] = nan;
}
if (numberp($uz[$r]) && $curGap>0) { # end of gap
push(@gaps,$curGap) unless ($r == $curGap); # do not report "gap" at beginning of profile
# print(STDERR "$curGap-gap at $depth[$r]m\n");
$curGap = 0;
} elsif (!numberp($uz[$r])) { # currently in gap
# print(STDERR "in gap at $depth[$r]m (nsamp = $nsamp[$r], $dc_nsamp,$uc_nsamp)\n");
$curGap++;
}
if ($nsamp[$r] > 0) {
if (defined($opt_u)) {
$uzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_uzsF],
$UL_[$dc_nshF], $UL_[$dc_uzsF],
$ants_[0][$uc_nshF],$ants_[0][$uc_uzsF],
$UL_[$uc_nshF], $UL_[$uc_uzsF]);
$vzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_vzsF],
$UL_[$dc_nshF], $UL_[$dc_vzsF],
$ants_[0][$uc_nshF],$ants_[0][$uc_vzsF],
$UL_[$uc_nshF], $UL_[$uc_vzsF]);
$wzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_wzsF],
$UL_[$dc_nshF], $UL_[$dc_wzsF],
$ants_[0][$uc_nshF],$ants_[0][$uc_wzsF],
$UL_[$uc_nshF], $UL_[$uc_wzsF]);
} else { # DL only
$uzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_uzsF],
$ants_[0][$uc_nshF],$ants_[0][$uc_uzsF]);
$vzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_vzsF],
$ants_[0][$uc_nshF],$ants_[0][$uc_vzsF]);
$wzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_wzsF],
$ants_[0][$uc_nshF],$ants_[0][$uc_wzsF]);
}
} else {
$uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan;
}
if ($dc_nsamp > 0) { # same calc for downcast only
if (defined($opt_u)) {
$dc_uzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_uzsF],
$UL_[$dc_nshF], $UL_[$dc_uzsF]);
$dc_vzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_vzsF],
$UL_[$dc_nshF], $UL_[$dc_vzsF]);
$dc_wzsig[$r] = wavg_sig($ants_[0][$dc_nshF],$ants_[0][$dc_wzsF],
$UL_[$dc_nshF], $UL_[$dc_wzsF]);
} else {
$dc_uzsig[$r] = $ants_[0][$dc_uzsF];
$dc_vzsig[$r] = $ants_[0][$dc_vzsF];
$dc_wzsig[$r] = $ants_[0][$dc_wzsF];
}
} else {
$uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan;
}
if ($uc_nsamp > 0) { # same calc for upcast only
if (defined($opt_u)) {
$uc_uzsig[$r] = wavg_sig($ants_[0][$uc_nshF],$ants_[0][$uc_uzsF],
$UL_[$uc_nshF], $UL_[$uc_uzsF]);
$uc_vzsig[$r] = wavg_sig($ants_[0][$uc_nshF],$ants_[0][$uc_vzsF],
$UL_[$uc_nshF], $UL_[$uc_vzsF]);
$uc_wzsig[$r] = wavg_sig($ants_[0][$uc_nshF],$ants_[0][$uc_wzsF],
$UL_[$uc_nshF], $UL_[$uc_wzsF]);
} else {
$uc_uzsig[$r] = $ants_[0][$uc_uzsF];
$uc_vzsig[$r] = $ants_[0][$uc_vzsF];
$uc_wzsig[$r] = $ants_[0][$uc_wzsF];
}
} else {
$uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan;
}
}
if (@gaps) {
&antsAddParams('shear_gaps',"@gaps");
print(STDERR "shear gaps: @gaps\n");
} else {
&antsAddParams('shear_gaps',0);
}
#===============================================================================
# Step 2: Low-Pass filter high-quality shear data; not yet implemented
#===============================================================================
#======================================================================
# Step 3: Integrate Shear
# - z(vel) = z(sh) + DZ/2
#======================================================================
my($DZ) = $depth[1] - $depth[0];
for (my($r)=my($u)=my($v)=my($w)=my($dc_u)=my($dc_v)=my($dc_w)=my($uc_u)=my($uc_v)=my($uc_w)=0;
$r<@depth; $r++) {
if ($nsamp[$r] >= $opt_n) {
$u = $u[$r] = $u + $DZ*$uz[$r];
$v = $v[$r] = $v + $DZ*$vz[$r];
$w = $w[$r] = $w + $DZ*$wz[$r];
}
if ($dc_nsamp[$r] >= $opt_n) {
$dc_u = $dc_u[$r] = $dc_u + $DZ*$dc_uz[$r];
$dc_v = $dc_v[$r] = $dc_v + $DZ*$dc_vz[$r];
$dc_w = $dc_w[$r] = $dc_w + $DZ*$dc_wz[$r];
}
if ($uc_nsamp[$r] >= $opt_n) {
$uc_u = $uc_u[$r] = $uc_u + $DZ*$uc_uz[$r];
$uc_v = $uc_v[$r] = $uc_v + $DZ*$uc_vz[$r];
$uc_w = $uc_w[$r] = $uc_w + $DZ*$uc_wz[$r];
}
}
#======================================================================
# Step 4: Reference Velocities
#======================================================================
my($refU,$refV,$refW,$dc_refU,$dc_refV,$dc_refW,$uc_refU,$uc_refV,$uc_refW);
if (defined($opt_r)) { # reference using velocity profile
print(STDERR "Loading reference-velocity data from $opt_r...\n")
if ($opt_d);
open(BTF,$opt_r) || croak("$opt_r: $!\n");
my(@BTL) = &antsFileLayout(BTF);
if (@BTL) { # valid ANTS file
my($BTdF,$BTndF,$BTuF,$BTvF,$BTwF,$BTusF,$BTvsF,$BTwsF);
for (my($f)=0; $f<@BTL; $f++) {
$BTdF = $f if ($BTL[$f] eq 'depth');
$BTuF = $f if ($BTL[$f] eq 'u');
$BTvF = $f if ($BTL[$f] eq 'v');
$BTwF = $f if ($BTL[$f] eq 'w');
$BTusF = $f if ($BTL[$f] eq 'u.sig');
$BTvsF = $f if ($BTL[$f] eq 'v.sig');
$BTwsF = $f if ($BTL[$f] eq 'w.sig');
$BTndF = $f if ($BTL[$f] eq 'nsamp');
$BTerrF= $f if ($BTL[$f] eq 'err');
}
if (defined($BTdF) && defined($BTuF) && # from LADCPproc
defined($BTvF) && defined($BTwF) &&
defined($BTusF) && defined($BTvsF) &&
defined($BTwsF) && defined($BTndF)) {
while (my(@BTr) = &antsFileIn(BTF)) {
my($gi) = int($BTr[$BTdF] / $DZ);
next unless ($BTr[$BTndF] >= $minBTsamp);
$BT_nsamp[$gi] = $BTr[$BTndF];
$BT_u[$gi] = $BTr[$BTuF];
$BT_v[$gi] = $BTr[$BTvF];
$BT_w[$gi] = $BTr[$BTwF];
$BT_u_var[$gi] = $BTr[$BTusF]**2;
$BT_v_var[$gi] = $BTr[$BTvsF]**2;
$BT_w_var[$gi] = $BTr[$BTwsF]**2;
}
&fixLowSampStat(\@BT_u_var,@BT_nsamp); # remove spurious small variances
&fixLowSampStat(\@BT_v_var,@BT_nsamp);
&fixLowSampStat(\@BT_w_var,@BT_nsamp);
} elsif (defined($BTdF) && defined($BTuF) && # LDEO_IX ANTS format
defined($BTvF) && defined($BTerrF)) {
croak("$0: -m not supported for LDEO_IX output\n")
if defined($opt_m);
while (my(@BTr) = &antsFileIn(BTF)) {
my($gi) = int($BTr[$BTdF] / $DZ);
$BT_u[$gi] = $BTr[$BTuF];
$BT_v[$gi] = $BTr[$BTvF];
$BT_u_var[$gi] = $BTr[$BTerrF]**2;
$BT_v_var[$gi] = $BTr[$BTerrF]**2;
# print(STDERR "$gi $BT_u[$gi] $BT_v[$gi] $BT_u_var[$gi] $BT_v_var[$gi]\n");
}
} else {
croak("$opt_r: not a valid reference-velocity file (ANTS format)\n");
}
} else { # non-ANTS file (LDEO_IX assumed)
croak("$0: -m not supported for LDEO_IX output\n")
if defined($opt_m);
while (<BTF>) {
last if /^Columns\s+=\s+z:u:v:err/;
}
croak("$opt_r: not a valid reference-velocity file (non-ANTS format)\n")
unless /^Columns\s+=\s+z:u:v:err/;
while (<BTF>) {
my($depth,$u,$v,$err) = split;
croak("$0: cannot handle non-numeric BT data (implementation restriction)\n")
unless numberp($depth) && numberp($u) && numberp($v) && numberp($err);
my($gi) = int($depth / $DZ);
$BT_u[$gi] = $u;
$BT_v[$gi] = $v;
$BT_u_var[$gi] = $err**2;
$BT_v_var[$gi] = $err**2;
# print(STDERR "$gi $BT_u[$gi] $BT_v[$gi] $BT_u_var[$gi] $BT_v_var[$gi]\n");
}
}
close(BTF);
my($sumU,$sumV,$sumW,$dc_sumU,$dc_sumV,$dc_sumW, # average integrated-shear velocities
$uc_sumU,$uc_sumV,$uc_sumW);
my($nSumVel,$dc_nSumVel,$uc_nSumVel);
my($wSumBTu,$wSumBTv,$wSumBTw); # weighted sums of BT-ref'd velocities
my($dc_wSumBTu,$dc_wSumBTv,$dc_wSumBTw);
my($uc_wSumBTu,$uc_wSumBTv,$uc_wSumBTw);
my($sumVarBTu,$sumVarBTv,$sumVarBTw); # sum of variances of BT-ref'd vels
my($dc_sumVarBTu,$dc_sumVarBTv,$dc_sumVarBTw);
my($uc_sumVarBTu,$uc_sumVarBTv,$uc_sumVarBTw);
for (my($r)=0; $r<@depth; $r++) {
if (numberp($BT_u[$r]) && numberp($u[$r])) {
$nSumVel++;
$sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r];
$wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $sumVarBTu += 1/$BT_u_var[$r];
$wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $sumVarBTv += 1/$BT_v_var[$r];
if (@BT_w) {
$wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $sumVarBTw += 1/$BT_w_var[$r];
}
}
if (numberp($BT_u[$r]) && numberp($dc_u[$r])) {
$dc_nSumVel++;
$dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r];
$dc_wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $dc_sumVarBTu += 1/$BT_u_var[$r];
$dc_wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $dc_sumVarBTv += 1/$BT_v_var[$r];
if (@BT_w) {
$dc_wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $dc_sumVarBTw += 1/$BT_w_var[$r];
}
}
if (numberp($BT_u[$r]) && numberp($uc_u[$r])) {
$uc_nSumVel++;
$uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r];
$uc_wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $uc_sumVarBTu += 1/$BT_u_var[$r];
$uc_wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $uc_sumVarBTv += 1/$BT_v_var[$r];
if (@BT_w) {
$uc_wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $uc_sumVarBTw += 1/$BT_w_var[$r];
}
}
}
if ($nSumVel > 0) {
$refU = $sumU/$nSumVel - $wSumBTu/$sumVarBTu;
$refV = $sumV/$nSumVel - $wSumBTv/$sumVarBTv;
$refW = $sumW/$nSumVel - $wSumBTw/$sumVarBTw if (@BT_w);
if ($dc_nSumVel > 0) {
$dc_refU = $dc_sumU/$dc_nSumVel - $dc_wSumBTu/$dc_sumVarBTu;
$dc_refV = $dc_sumV/$dc_nSumVel - $dc_wSumBTv/$dc_sumVarBTv;
$dc_refW = $dc_sumW/$dc_nSumVel - $dc_wSumBTw/$dc_sumVarBTw if (@BT_w);
} else {
&antsInfo("$opt_r: insufficient reference-velocity data to constrain DC profile --- baroclinic profile only");
}
if ($uc_nSumVel > 0) {
$uc_refU = $uc_sumU/$uc_nSumVel - $uc_wSumBTu/$uc_sumVarBTu;
$uc_refV = $uc_sumV/$uc_nSumVel - $uc_wSumBTv/$uc_sumVarBTv;
$uc_refW = $uc_sumW/$uc_nSumVel - $uc_wSumBTw/$uc_sumVarBTw if (@BT_w);
} else {
&antsInfo("$opt_r: insufficient reference-velocity data to constrain UC profile --- baroclinic profile only");
}
} else {
&antsInfo("$opt_r: no valid reference-velocity data --- baroclinic profiles only");
}
}
unless (defined($refU)) { # no reference velocity => use zero mean
my($sumU,$sumV,$sumW,$nSumVel);
for (my($r)=0; $r<@depth; $r++) {
if (numberp($u[$r])) {
$nSumVel++;
$sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r] if (@BT_w);
}
}
$refU = $sumU / $nSumVel; $refV = $sumV / $nSumVel; $refW = $sumW / $nSumVel if (@BT_w);
}
unless (defined($dc_refU)) {
my($dc_sumU,$dc_sumV,$dc_sumW,$dc_nSumVel);
for (my($r)=0; $r<@depth; $r++) {
if (numberp($dc_u[$r])) {
$dc_nSumVel++;
$dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r] if (@BT_w);
}
}
if ($dc_nSumVel) {
$dc_refU = $dc_sumU / $dc_nSumVel; $dc_refV = $dc_sumV / $dc_nSumVel;
$dc_refW = $dc_sumW / $dc_nSumVel if (@BT_w);
}
}
unless (defined($uc_refU)) {
my($uc_sumU,$uc_sumV,$uc_sumW,$uc_nSumVel);
for (my($r)=0; $r<@depth; $r++) {
if (numberp($uc_u[$r])) {
$uc_nSumVel++;
$uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r] if (@BT_w);
}
}
if ($uc_nSumVel) {
$uc_refU = $uc_sumU / $uc_nSumVel; $uc_refV = $uc_sumV / $uc_nSumVel;
$uc_refW = $uc_sumW / $uc_nSumVel if (@BT_w);
}
}
for (my($r)=0; $r<@depth; $r++) { # reference velocities
$u[$r] -= $refU if defined($u[$r]);
$v[$r] -= $refV if defined($v[$r]);
$w[$r] -= $refW if defined($w[$r]);
$dc_u[$r] -= $dc_refU if defined($dc_u[$r]);
$dc_v[$r] -= $dc_refV if defined($dc_v[$r]);
$dc_w[$r] -= $dc_refW if defined($dc_w[$r]);
$uc_u[$r] -= $uc_refU if defined($uc_u[$r]);
$uc_v[$r] -= $uc_refV if defined($uc_v[$r]);
$uc_w[$r] -= $uc_refW if defined($uc_w[$r]);
}
#======================================================================
# Determine X Factor
#======================================================================
if ($dc_nSumVel && $uc_nSumVel) {
my($first_w,$last_w);
for (my($r)=0; !defined($first_w) || !defined($last_w); $r++) {
$first_w = $dc_w[$r] unless defined($first_w);
$last_w = $uc_w[$r] unless defined($last_w);
}
my($X_Factor) = 100 * abs($last_w-$first_w) / sqrt(@depth / $DZ);
&antsAddParams('X-Factor',$X_Factor);
printf(STDERR "X-Factor = %.1f\n",$X_Factor);
}
#======================================================================
# Output Velocity Profile
#======================================================================
@antsNewLayout = ('depth','elapsed','u','v','w','nsamp',
'dc_elapsed','dc_u','dc_v','dc_w','dc_nsamp',
'uc_elapsed','uc_u','uc_v','uc_w','uc_nsamp');
for (my($r)=0; $r<@depth; $r++) {
&antsOut($depth[$r]+$DZ/2,
$elapsed[$r],$u[$r],$v[$r],$w[$r],$nsamp[$r],
$dc_elapsed[$r],$dc_u[$r],$dc_v[$r],$dc_w[$r],$dc_nsamp[$r],
$uc_elapsed[$r],$uc_u[$r],$uc_v[$r],$uc_w[$r],$uc_nsamp[$r]);
}
#======================================================================
# Output Averaged Shear Profile
#======================================================================
if (defined($opt_s)) {
@antsNewLayout = ('depth','elapsed','u_z','v_z','w_z','u_z.sig','v_z.sig','w_z.sig','nsamp',
'dc_elapsed','dc_u_z','dc_v_z','dc_w_z','dc_u_z.sig','dc_v_z.sig','dc_w_z.sig','dc_nsamp',
'uc_elapsed','uc_u_z','uc_v_z','uc_w_z','uc_u_z.sig','uc_v_z.sig','uc_w_z.sig','uc_nsamp');
&antsOut('EOF');
close(STDOUT);
open(STDOUT,">$opt_s") || croak("$opt_s: $!\n");
for (my($r)=0; $r<@depth; $r++) {
&antsOut($depth[$r],$elapsed[$r],$uz[$r],$vz[$r],$wz[$r],$uzsig[$r],$vzsig[$r],$wzsig[$r],$nsamp[$r],
$dc_elapsed[$r],$dc_uz[$r],$dc_vz[$r],$dc_wz[$r],$dc_uzsig[$r],$dc_vzsig[$r],$dc_wzsig[$r],$dc_nsamp[$r],
$uc_elapsed[$r],$uc_uz[$r],$uc_vz[$r],$uc_wz[$r],$uc_uzsig[$r],$uc_vzsig[$r],$uc_wzsig[$r],$uc_nsamp[$r]);
}
}
&antsExit();