diff --git a/LADCPintsh b/LADCPintsh --- a/LADCPintsh +++ b/LADCPintsh @@ -2,9 +2,9 @@ #====================================================================== # L A D C P I N T S H # doc: Thu Oct 14 21:22:50 2010 -# dlm: Sun Oct 24 16:42:59 2010 +# dlm: Fri Dec 10 04:57:50 2010 # (c) 2010 A.M. Thurnherr & E. Firing -# uE-Info: 338 29 NIL 0 0 72 2 2 4 NIL ofnI +# uE-Info: 207 1 NIL 0 0 72 2 2 4 NIL ofnI #====================================================================== $antsSummary = 'integrate LADCP shear'; @@ -27,14 +27,16 @@ # 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 ($ANTS) = (`which list` =~ m{^(.*)/[^/]*$}); require "$ANTS/ants.pl"; require "$ANTS/libstats.pl"; -&antsUsage('b:dr:s:w:',0, +&antsUsage('b:dr:s:w:u:',0, '[-d)ebug]', '[reference with -b)ottom-track ]', + '[-u)plooker ]', '[min -r)aw ] [min -s)moothed ]', '[-w)eighted avg ]', '[LADCP shear file]'); @@ -42,7 +44,12 @@ &antsCardOpt(\$opt_r,10); ## minimum number of samples for shear &antsCardOpt(\$opt_s,20); ## minimum number of samples for smoothed shear &antsCardOpt(\$opt_w,1200); ## probably 1200 for midlatitude-equatorial, 2400 for high lat -&antsFileOpt($opt_b); + +&antsFileOpt($opt_b); # BT file + +&antsFileOpt($opt_u); # UL shear file +open(ULF,$opt_u) || croak("$opt_u: $!\n") + if defined($opt_u); #====================================================================== # Step 1: Read and Average Shear Data @@ -77,61 +84,116 @@ $uc_wzsF = fnr('uc_w_z_sig') unless defined($uc_wzsF); -croak("$0: empty input\n") - unless &antsIn(); - my(@gaps); my($curGap) = 0; -for (my($r)=0; &antsIn(); $r++) { - $depth[$r] = $ants_[0][$depthF]; ## depth grid values - $ndata[$r] = $ants_[0][$dc_nshF] + $ants_[0][$uc_nshF]; ## number of points - if ($ants_[0][$dc_nshF] >= $opt_r) { # downcast shear - $dc_uz[$r] = $ants_[0][$dc_uzF]; - $dc_vz[$r] = $ants_[0][$dc_vzF]; - $dc_wz[$r] = $ants_[0][$dc_wzF]; - } else { - $dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = nan; +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 (DL: $depth[$r]; UL: $UL_[$depthF])\n") + if defined($opt_u) && ($UL_[$depthF] != $depth[$r]); + + $dc_ndata = $ants_[0][$dc_nshF]; # number of shear samples + $uc_ndata = $ants_[0][$uc_nshF]; + if (defined($opt_u)) { + $dc_ndata += $UL_[$dc_nshF]; + $uc_ndata += $UL_[$uc_nshF]; + } + $ndata[$r] = $dc_ndata + $uc_ndata; - if ($ants_[0][$uc_nshF] >= $opt_r) { # upcast shear - $uc_uz[$r] = $ants_[0][$uc_uzF]; - $uc_vz[$r] = $ants_[0][$uc_vzF]; - $uc_wz[$r] = $ants_[0][$uc_wzF]; + if (defined($opt_u)) { + if ($dc_ndata > 0) { # downcast shear + my($DLf) = $ants_[0][$dc_nshF] / $dc_ndata; + my($ULf) = $UL_[$dc_nshF] / $dc_ndata; + 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]; + } elsif ($DLf > 0) { + $dc_uz[$r] = $ants_[0][$dc_uzF]; + $dc_vz[$r] = $ants_[0][$dc_vzF]; + $dc_wz[$r] = $ants_[0][$dc_wzF]; + } else { + $dc_uz[$r] = $UL_[$dc_uzF]; + $dc_vz[$r] = $UL_[$dc_vzF]; + $dc_wz[$r] = $UL_[$dc_wzF]; + } + } else { + $dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = nan; + } + if ($uc_ndata > 0) { # upcast shear + my($DLf) = $ants_[0][$uc_nshF] / $uc_ndata; + my($ULf) = $UL_[$uc_nshF] / $uc_ndata; + 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]; + } elsif ($DLf > 0) { + $uc_uz[$r] = $ants_[0][$uc_uzF]; + $uc_vz[$r] = $ants_[0][$uc_vzF]; + $uc_wz[$r] = $ants_[0][$uc_wzF]; + } else { + $uc_uz[$r] = $UL_[$uc_uzF]; + $uc_vz[$r] = $UL_[$uc_vzF]; + $uc_wz[$r] = $UL_[$uc_wzF]; + } + } else { + $uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = nan; + } } else { - $uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = nan; - } - + if ($dc_ndata >= $opt_r) { # downcast shear + $dc_uz[$r] = $ants_[0][$dc_uzF]; + $dc_vz[$r] = $ants_[0][$dc_vzF]; + $dc_wz[$r] = $ants_[0][$dc_wzF]; + } else { + $dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = nan; + } + if ($uc_ndata >= $opt_r) { # upcast shear + $uc_uz[$r] = $ants_[0][$uc_uzF]; + $uc_vz[$r] = $ants_[0][$uc_vzF]; + $uc_wz[$r] = $ants_[0][$uc_wzF]; + } else { + $uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = nan; + } + } + if ($depth[$r] < $opt_w) { ## use average of up and down casts at top of profile - if ($ants_[0][$dc_nshF]>0 && $ants_[0][$uc_nshF]>0 && $ndata[$r]>=$opt_r) { - $uz[$r] = $ants_[0][$dc_uzF]/2 + $ants_[0][$uc_uzF]/2; - $vz[$r] = $ants_[0][$dc_vzF]/2 + $ants_[0][$uc_vzF]/2; - $wz[$r] = $ants_[0][$dc_wzF]/2 + $ants_[0][$uc_wzF]/2; - } elsif ($ants_[0][$dc_nshF] >= $opt_r) { - $uz[$r] = $ants_[0][$dc_uzF]; - $vz[$r] = $ants_[0][$dc_vzF]; - $wz[$r] = $ants_[0][$dc_wzF]; - } elsif ($ants_[0][$uc_nshF] >= $opt_r) { - $uz[$r] = $ants_[0][$uc_uzF]; - $vz[$r] = $ants_[0][$uc_vzF]; - $wz[$r] = $ants_[0][$uc_wzF]; + if ($dc_ndata>0 && $uc_ndata>0 && $ndata[$r]>=$opt_r) { + $uz[$r] = $dc_uz[$r]/2 + $uc_uz[$r]/2; + $vz[$r] = $dc_vz[$r]/2 + $uc_vz[$r]/2; + $wz[$r] = $dc_wz[$r]/2 + $uc_wz[$r]/2; + } elsif ($dc_ndata >= $opt_r) { + $uz[$r] = $dc_uz[$r]; + $vz[$r] = $dc_vz[$r]; + $wz[$r] = $dc_wz[$r]; + } elsif ($uc_ndata >= $opt_r) { + $uz[$r] = $uc_uz[$r]; + $vz[$r] = $uc_vz[$r]; + $wz[$r] = $uc_wz[$r]; } else { $uz[$r] = $vz[$r] = $wz[$r] = nan; } } else { ## use weighted average of up and down cast data at bottom of profile - if ($ants_[0][$dc_nshF]>0 && $ants_[0][$uc_nshF]>0 && $ndata[$r]>=$opt_r) { - my($dcf) = $ants_[0][$dc_nshF] / $ndata[$r]; - my($ucf) = $ants_[0][$uc_nshF] / $ndata[$r]; - $uz[$r] = $dcf*$ants_[0][$dc_uzF] + $ucf*$ants_[0][$uc_uzF]; - $vz[$r] = $dcf*$ants_[0][$dc_vzF] + $ucf*$ants_[0][$uc_vzF]; - $wz[$r] = $dcf*$ants_[0][$dc_wzF] + $ucf*$ants_[0][$uc_wzF]; - } elsif ($ants_[0][$dc_nshF] >= $opt_r) { - $uz[$r] = $ants_[0][$dc_uzF]; - $vz[$r] = $ants_[0][$dc_vzF]; - $wz[$r] = $ants_[0][$dc_wzF]; - } elsif ($ants_[0][$uc_nshF] >= $opt_r) { - $uz[$r] = $ants_[0][$uc_uzF]; - $vz[$r] = $ants_[0][$uc_vzF]; - $wz[$r] = $ants_[0][$uc_wzF]; + if ($ndata[$r] >= $opt_r) { + my($dcf) = $dc_ndata / $ndata[$r]; + my($ucf) = $uc_ndata / $ndata[$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]; + } elsif ($dcf > 0) { + $uz[$r] = $dc_uz[$r]; + $vz[$r] = $dc_vz[$r]; + $wz[$r] = $dc_wz[$r]; + } else { + $uz[$r] = $uc_uz[$r]; + $vz[$r] = $uc_vz[$r]; + $wz[$r] = $uc_wz[$r]; + } } else { $uz[$r] = $vz[$r] = $wz[$r] = nan; } @@ -142,16 +204,32 @@ # 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 (ndata = $ndata[$r], $dc_ndata,$uc_ndata)\n"); $curGap++; } if ($ndata[$r] >= $opt_r) { ## back to sum of squares - $uzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_uzsF]**2 + - $ants_[0][$uc_nshF] * $ants_[0][$uc_uzsF]**2) / $ndata[$r]); - $vzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_vzsF]**2 + - $ants_[0][$uc_nshF] * $ants_[0][$uc_vzsF]**2) / $ndata[$r]); - $wzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_wzsF]**2 + - $ants_[0][$uc_nshF] * $ants_[0][$uc_wzsF]**2) / $ndata[$r]); + if (defined($opt_u)) { + $uzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_uzsF]**2 + + $UL_[$dc_nshF] * $UL_[$dc_uzsF]**2 + + $ants_[0][$uc_nshF] * $ants_[0][$uc_uzsF]**2 + + $UL_[$uc_nshF] * $UL_[$uc_uzsF]**2) / $ndata[$r]); + $vzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_vzsF]**2 + + $UL_[$dc_nshF] * $UL_[$dc_vzsF]**2 + + $ants_[0][$uc_nshF] * $ants_[0][$uc_vzsF]**2 + + $UL_[$uc_nshF] * $UL_[$uc_vzsF]**2) / $ndata[$r]); + $wzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_wzsF]**2 + + $UL_[$dc_nshF] * $UL_[$dc_wzsF]**2 + + $ants_[0][$uc_nshF] * $ants_[0][$uc_wzsF]**2 + + $UL_[$uc_nshF] * $UL_[$uc_wzsF]**2) / $ndata[$r]); + } else { + $uzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_uzsF]**2 + + $ants_[0][$uc_nshF] * $ants_[0][$uc_uzsF]**2) / $ndata[$r]); + $vzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_vzsF]**2 + + $ants_[0][$uc_nshF] * $ants_[0][$uc_vzsF]**2) / $ndata[$r]); + $wzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_wzsF]**2 + + $ants_[0][$uc_nshF] * $ants_[0][$uc_wzsF]**2) / $ndata[$r]); + } } else { $uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan; } @@ -281,23 +359,24 @@ } } - croak("$opt_b: no valid BT data\n") - unless ($nSumVel > 0); + if ($nSumVel > 0) { + $refU = $sumU/$nSumVel - $wSumBTu/$sumVarBTu; + $refV = $sumV/$nSumVel - $wSumBTv/$sumVarBTv; + $refW = $sumW/$nSumVel - $wSumBTw/$sumVarBTw; + + $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; + + $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; + } else { + &antsInfo("$opt_b: no valid BT data --- calculating baroclinic solution only"); + } +} -# printf(STDERR "v.mu = %f; BT_v.mu = %f(=$wSumBTv/$sumVarBTv); $nSumVel samples\n",$sumV/$nSumVel,$wSumBTv/$sumVarBTv); - $refU = $sumU/$nSumVel - $wSumBTu/$sumVarBTu; - $refV = $sumV/$nSumVel - $wSumBTv/$sumVarBTv; - $refW = $sumW/$nSumVel - $wSumBTw/$sumVarBTw; - - $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; - - $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; - -} else { # baroclinic (zero mean) velocities +unless (defined($refU)) { my($sumU,$sumV,$sumW,$dc_sumU,$dc_sumV,$dc_sumW,$uc_sumU,$uc_sumV,$uc_sumW); for (my($r)=0; $r<@depth; $r++) {