beginning of DIMES US5 cruise
authorA.M. Thurnherr <athurnherr@yahoo.com>
Mon, 11 Nov 2013 11:20:45 -0400
changeset 22 f6635c0384b7
parent 21 ec19ba9622f3
child 23 85c8e2ea2a5b
beginning of DIMES US5 cruise
LADCPproc
LADCPproc.BT
LADCPproc.UHcode
LADCPproc.defaults
LADCPproc.loadCTD
--- a/LADCPproc
+++ b/LADCPproc
@@ -2,9 +2,9 @@
 #======================================================================
 #                    L A D C P P R O C 
 #                    doc: Thu Sep 16 20:36:10 2010
-#                    dlm: Tue Jun 25 14:28:32 2013
+#                    dlm: Wed Sep 25 13:09:21 2013
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 81 54 NIL 0 0 72 10 2 4 NIL ofnI
+#                    uE-Info: 84 61 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # NOTES:
@@ -79,6 +79,9 @@
 #	Jun  5, 2013: - added $bad_beam support
 #	Jun 25, 2013: - added %PARAMS used for spectral correction
 #				  - adapted to new ::-PARAM convention
+#	Sep 25, 2013: - BUG: %PARAM magnetic_declination did not have LADCPproc:: prefix
+#				  - added CTD lat/lon info to most output files (but not BT)
+#				  - BUG: moved %water_depth to common %PARAMs
 
 ($ANTS)    = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
 ($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -138,9 +141,9 @@
 }
 
 if (defined($opt_g)) {
-    ($CTD{lat},$CTD{lon}) = split(',',$opt_g);
+    ($CTD{stn_lat},$CTD{stn_lon}) = split(',',$opt_g);
     croak("$0: cannot decode -g $opt_g\n")
-        unless numberp($CTD{lat}) && numberp($CTD{lon});
+        unless numberp($CTD{stn_lat}) && numberp($CTD{stn_lon});
 }
 
 if (defined($opt_c)) {
@@ -209,12 +212,12 @@
 my($year)  = substr($LADCP{ENSEMBLE}[0]->{DATE},6,4);
 my($month) = substr($LADCP{ENSEMBLE}[0]->{DATE},0,2);
 my($day  ) = substr($LADCP{ENSEMBLE}[0]->{DATE},3,2);
-my($magdec,$maginc,$h_strength,$v_strength) = split('\s+',`magdec $CTD{lon} $CTD{lat} $year $month $day`);
+my($magdec,$maginc,$h_strength,$v_strength) = split('\s+',`magdec $CTD{stn_lon} $CTD{stn_lat} $year $month $day`);
 
 croak("$0: unknown magnetic declination\n")
 	unless defined($magdec);
 
-&antsAddParams('magnetic_declination',$magdec);
+&antsAddParams('LADCPproc::magnetic_declination',$magdec);
 
 #----------------------------------------------------------------------
 # Step 4: Pre-Process CTD & LADCP Data
@@ -424,7 +427,7 @@
 	}
 	my($dr);
 	for ($dr=0; !numberp($CTD{press}[$r+$dr]); $dr--) {}
-	$LADCP{ENSEMBLE}[$ens]->{DEPTH} = depth($CTD{press}[$r+$dr],$CTD{lat});
+	$LADCP{ENSEMBLE}[$ens]->{DEPTH} = depth($CTD{press}[$r+$dr],$CTD{stn_lat});
 	if ($LADCP{ENSEMBLE}[$ens]->{DEPTH} < $min_depth) {
 		$min_depth = $LADCP{ENSEMBLE}[$ens]->{DEPTH};
 		$LADCP_top = $ens;
@@ -444,16 +447,24 @@
 		$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V] *= $sscorr;
 		$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W] *= $sscorr;
     }
+    $LADCP{ENSEMBLE}[$ens]->{CTD_LAT} = $CTD{lat}[$r];
+    $LADCP{ENSEMBLE}[$ens]->{CTD_LON} = $CTD{lon}[$r];
 }
 
 &antsAddParams('LADCPproc::min_depth',round($LADCP{ENSEMBLE}[$LADCP_top]->{DEPTH}),
 			   'LADCPproc::max_depth',round($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH}),
 			   'LADCPproc::start_date',$LADCP{ENSEMBLE}[$LADCP_start]->{DATE},
 			   'LADCPproc::start_time',$LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
+			   'LADCPproc::start_lat',$LADCP{ENSEMBLE}[$LADCP_start]->{CTD_LAT},
+			   'LADCPproc::start_lon',$LADCP{ENSEMBLE}[$LADCP_start]->{CTD_LON},
 			   'LADCPproc::bottom_date',$LADCP{ENSEMBLE}[$LADCP_bottom]->{DATE},
 			   'LADCPproc::bottom_time',$LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
+			   'LADCPproc::bottom_lat',$LADCP{ENSEMBLE}[$LADCP_bottom]->{CTD_LAT},
+			   'LADCPproc::bottom_lon',$LADCP{ENSEMBLE}[$LADCP_bottom]->{CTD_LON},
 			   'LADCPproc::end_date',$LADCP{ENSEMBLE}[$LADCP_end]->{DATE},
-			   'LADCPproc::end_time',$LADCP{ENSEMBLE}[$LADCP_end]->{TIME});
+			   'LADCPproc::end_time',$LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
+			   'LADCPproc::end_lat',$LADCP{ENSEMBLE}[$LADCP_end]->{CTD_LAT},
+			   'LADCPproc::end_lon',$LADCP{ENSEMBLE}[$LADCP_end]->{CTD_LON});
 
 print(STDERR "\n");
 
@@ -513,13 +524,16 @@
 		if ($opt_d);
 }
 
+&antsAddParams('LADCPproc::water_depth',round($water_depth),
+			   'LADCPproc::water_depth.sig',round($sig_water_depth));
+
 #-----------------------------------------------------------------------
 # Step 8: Edit Data & also produce depth-time-series output via callback
 #-----------------------------------------------------------------------
 
 print(STDERR "Calculating shear profiles & producing time-depth-series (.tds) output...");
 
-@antsNewLayout = ('ens','elapsed','CTD_depth','downcast','depth','u_z','v_z','w_z');
+@antsNewLayout = ('ens','elapsed','CTD_depth','downcast','CTD_lat','CTD_lon','depth','u_z','v_z','w_z');
 
 	#--------------------------------------------------------------------------------
 	# callback routine to output .tds data, called once each for down-/upcasts after
@@ -544,6 +558,8 @@
 					&antsOut($ens_vals[$gi][$i],
 							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
 							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
+							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
+							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
 							 depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
 				}
 	        }
@@ -553,6 +569,8 @@
 					&antsOut($ens_vals[$gi][$i],
 							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
 							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{DEPTH},$downcast,
+							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LAT},
+							 $LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{CTD_LON},
 							 depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
 				}
 	        }
@@ -574,14 +592,13 @@
 			   'LADCPproc::min_wake_w',$min_wake_w,'LADCPproc::n_wake_bins',$n_wake_bins,
 			   'LADCPproc::e_max',$e_max,'LADCPproc::min_cor',$min_cor,
 			   'LADCPproc::max_shdev',$max_shdev,'LADCPproc::max_shdev_sum',$max_shdev_sum,
-			   'LADCPproc::water_depth',round($water_depth),'LADCPproc::water_depth.sig',round($sig_water_depth),
 			   'LADCPproc::min_hab',round($min_hab),'LADCPproc::PPI_editing_enabled',$PPI_editing_enabled,
 			   'LADCPproc::clip_margin',$clip_margin,'LADCPproc::first_clip_bin',$first_clip_bin,
 			   'LADCPproc::Svbin_start',$Svbin_start,'LADCPproc::Svbin_end',$Svbin_end,
 			   'LADCPproc::BT_bin_start',$BT_bin_start,'LADCPproc::BT_bin_search_above',$BT_bin_search_above,
 			   'LADCPproc::BT_max_bin_spread',$BT_max_bin_spread,'LADCPproc::BT_max_w_difference',$BT_max_w_difference,
 );
-if (defined($BT_min_depth)) {
+if (defined($BT_min_depth)) {										# BT-related params
 	&antsAddParams('LADCPproc::BT_min_depth',$BT_min_depth,'LADCPproc::BT_max_depth',$BT_max_depth);
 } else {
 	&antsAddParams('LADCPproc::BT_max_depth_error',$BT_max_depth_error);
@@ -600,6 +617,7 @@
 @dc_vsh_mu = @vsh_mu; @dc_vsh_sig = @vsh_sig;
 @dc_wsh_mu = @wsh_mu; @dc_wsh_sig = @wsh_sig;
 @dc_esh_mu = @esh_mu;
+@dc_lash_mu = @lash_mu; @dc_losh_mu = @losh_mu;
 
 print(STDERR "\n\tupcast...") if ($opt_d);
 edit_velocity($LADCP_end,$LADCP_bottom);							# upcast
@@ -611,6 +629,7 @@
 @uc_vsh_mu = @vsh_mu; @uc_vsh_sig = @vsh_sig;
 @uc_wsh_mu = @wsh_mu; @uc_wsh_sig = @wsh_sig;
 @uc_esh_mu = @esh_mu;
+@uc_lash_mu = @lash_mu; @uc_losh_mu = @losh_mu;
 
 print(STDERR "\n\tcombined...") if ($opt_d);
 my($nsh) = (@dc_ush_mu > @uc_ush_mu) ? scalar(@dc_ush_mu) : scalar(@uc_ush_mu);
@@ -659,8 +678,8 @@
 
 	print(STDERR "Writing shear profiles...");
 	
-	@antsNewLayout = ('depth','dc_elapsed','dc_nsamp','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
-							  'uc_elapsed','uc_nsamp','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
+	@antsNewLayout = ('depth','dc_elapsed','dc_lat','dc_lon','dc_nsamp','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
+							  'uc_elapsed','uc_lat','uc_lon','uc_nsamp','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
 							  'elapsed','nsamp','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.nsamp');
 							  
 	&antsOut('EOF');
@@ -671,12 +690,12 @@
 
 	for (my($gi)=0; $gi<$nsh; $gi++) {
 		&antsOut(depthOfGI($gi),										# depth in center of bin
-				 $dc_esh_mu[$gi],										# downcast
+				 $dc_esh_mu[$gi],$dc_lash_mu[$gi],$dc_losh_mu[$gi],		# downcast
 				 numberp($dc_sh_n[$gi])?$dc_sh_n[$gi]:0,
 				 $dc_ush_mu[$gi],$dc_ush_sig[$gi],
 				 $dc_vsh_mu[$gi],$dc_vsh_sig[$gi],
 				 $dc_wsh_mu[$gi],$dc_wsh_sig[$gi],
-				 $uc_esh_mu[$gi],										# upcast
+				 $uc_esh_mu[$gi],$uc_lash_mu[$gi],$uc_losh_mu[$gi],		# upcast
 				 numberp($uc_sh_n[$gi])?$uc_sh_n[$gi]:0,
 				 $uc_ush_mu[$gi],$uc_ush_sig[$gi],
 				 $uc_vsh_mu[$gi],$uc_vsh_sig[$gi],
@@ -701,7 +720,7 @@
 	print(STDERR "Writing acoustic backscatter depth-time-series to <$opt_a>...");
 
 	
-	@antsNewLayout = ('ens','elapsed','CTD_depth','depth','bin','downcast',
+	@antsNewLayout = ('ens','elapsed','CTD_depth','CTD_lat','CTD_lon','depth','bin','downcast',
 					  'Sv_beam1','Sv_beam2','Sv_beam3','Sv_beam4','Sv.median');
 	&antsOut('EOF');
     close(STDOUT);
@@ -721,6 +740,7 @@
 			&antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
 					 $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
 					 $LADCP{ENSEMBLE}[$ens]->{DEPTH},
+					 $LADCP{ENSEMBLE}[$ens]->{CTD_LAT},$LADCP{ENSEMBLE}[$ens]->{CTD_LON},
 					 &depthOfBin($ens,$bin),$bin+1,
 					 ($ens <= $LADCP_bottom) ? 1 : 0,
 					 @{$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]},
@@ -737,7 +757,7 @@
 if (defined($opt_t)) {
 	print(STDERR "Writing time series to <$opt_t>...");
 	
-	@antsNewLayout = ('ens','elapsed','depth','CTD_w','LADCP_w');
+	@antsNewLayout = ('ens','elapsed','depth','CTD_lat','CTD_lon','CTD_w','LADCP_w');
 	&antsOut('EOF');
 	$antsCurParams = $commonParams;
 	close(STDOUT);
@@ -747,6 +767,8 @@
 		&antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
 				 $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l,
 				 $LADCP{ENSEMBLE}[$ens]->{DEPTH},
+				 $LADCP{ENSEMBLE}[$ens]->{CTD_LAT},
+				 $LADCP{ENSEMBLE}[$ens]->{CTD_LON},
 				 $LADCP{ENSEMBLE}[$ens]->{CTD_W},
 				 $LADCP{ENSEMBLE}[$ens]->{W});
 	}
--- a/LADCPproc.BT
+++ b/LADCPproc.BT
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L A D C P P R O C . B T 
 #                    doc: Wed Oct 20 21:05:37 2010
-#                    dlm: Wed May 16 18:45:15 2012
+#                    dlm: Thu Sep 19 13:17:41 2013
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 16 51 NIL 0 0 72 10 2 4 NIL ofnI
+#                    uE-Info: 119 96 NIL 0 0 72 10 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -14,6 +14,7 @@
 #				  - added $BT processing parameters
 #				  - changed from echo amplitude to Sv
 #	May 16, 2012: - added support for -r)DI BT data
+#	Sep 19, 2013: - added support for $BT_range_method
 
 my($BEAM1) = 0;
 my($BEAM2) = 1;
@@ -62,36 +63,66 @@
 		unless defined($BT_min_depth) ||
 			   (abs($water_depth-depthOfBin($ens,$range_bin)) < $sig_water_depth + $BT_max_depth_error);
 
-	# try bin of max plus one above and below
-	# this does not really work because, often, only one of the bins has valid velocities
-	my($w1) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];
-	my($w2) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
-	my($w3) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
-
-	printf(STDERR "w123 = %.1f,%.1f,%.1f\n",$w1,$w2,$w3)
-		if ($DEBUG);
-
-	$w1 = 9e99 unless numberp($w1);									# no valid velocities
-	$w2 = 9e99 unless numberp($w1);
-	$w3 = 9e99 unless numberp($w1);
-
 	my($CTD_u,$CTD_v,$CTD_w);
+	
+	if ($BT_range_method == 0) {										# take BT vel from bin with Sv max
+		
+		$nBTvalidVelFlag++,return unless numberp($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W]);
+		$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$U];
+		$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$V];
+		$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
+		
+	} elsif ($BT_range_method == 1) {									# take "best-fit" with w_reflr
 
-	if (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2) &&
-		abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w3)) {
-			$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$U];
-			$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$V];
-			$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];
-	} elsif (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2)) {
-			$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$U];
-			$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$V];
-			$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
+		# try bin of max plus one above and below
+		# this does not really work because, often, only one of the bins has valid velocities
+		my($w1) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];
+		my($w2) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
+		my($w3) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
+	
+		printf(STDERR "w123 = %.1f,%.1f,%.1f\n",$w1,$w2,$w3)
+			if ($DEBUG);
+	
+		$w1 = 9e99 unless numberp($w1); 								# no valid velocities
+		$w2 = 9e99 unless numberp($w1);
+		$w3 = 9e99 unless numberp($w1);
+	
+		if (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2) &&
+			abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w3)) {
+				$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$U];
+				$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$V];
+				$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];
+		} elsif (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2)) {
+				$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$U];
+				$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$V];
+				$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
+		} else {
+				$nBTvalidVelFlag++,return if ($w2 == 9e99); 			# none of 3 bins has valid velocity
+				$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$U];
+				$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$V];
+				$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
+	    }
+
+	} elsif ($BT_range_method == 2) {									# Visbeck method (median from 3 bins)
+		croak("$0: need \$BT_range_Visbeck_center\n")
+			unless defined($BT_range_Visbeck_center);
+
+		$CTD_u = median($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center+1][$U],
+						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center][$U],
+						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center-1][$U]);
+		$nBTvalidVelFlag++,return unless numberp($CTD_u);
+		$CTD_v = median($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center+1][$V],
+						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center][$V],
+						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center-1][$V]);
+		$CTD_w = median($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center+1][$W],
+						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center][$W],
+						$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+$BT_range_Visbeck_center-1][$W]);
+		
 	} else {
-			$nBTvalidVelFlag++,return if ($w2 == 9e99);				# none of 3 bins has valid velocity
-			$CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$U];
-			$CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$V];
-			$CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
+		
+		croak("$0: unknown \$BT_range_method == $BT_range_method\n");
 	}
+			
 
 	return ($CTD_u,$CTD_v,$CTD_w);
 }
--- a/LADCPproc.UHcode
+++ b/LADCPproc.UHcode
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L A D C P P R O C . U H C O D E 
 #                    doc: Fri Sep 17 20:27:53 2010
-#                    dlm: Wed Apr 11 19:22:44 2012
+#                    dlm: Wed Sep 25 12:29:46 2013
 #                    (c) 2010 A.M. Thurnherr & E. Firing
-#                    uE-Info: 38 46 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 39 47 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # PERLified functions from Eric's [merge.c]; with mods
@@ -36,6 +36,7 @@
 #	Jul 12, 2011: - replaced -p by $PPI_editing_enabled flag
 #	Feb 19, 2012: - added elapsed time to binned shear output
 #	Apr 11, 2012: - added MISSING_CTD_DATA_BIT
+#	Sep 25, 2013: - added code to calc gridded lat/lon info
 
 #======================================================================
 # VELOCITY EDITING
@@ -286,6 +287,7 @@
 # CALCULATE VELOCITY SHEAR
 #	- final output in @ush_mu,@vsh_mu,@wsh_mu,@ush_sig,@vsh_sig,@wsh_sig
 #		NEW (ant): elapsed time output in @esh_mu
+#				   lat/lon output in @lash_mu, @losh_mu
 # 	- @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh are defined "local" in calc_shear
 #======================================================================
 
@@ -408,11 +410,11 @@
 
 	outTDseries($De==1) if ($edit_shear);								# output depth-time time series
 
-	@ush_mu  = @vsh_mu  = @wsh_mu  = @esh_mu  = ();
+	@ush_mu  = @vsh_mu  = @wsh_mu  = @esh_mu = @lash_mu = @losh_mu = ();
 	@ush_sig = @vsh_sig = @wsh_sig = ();
 
 	for (my($gi)=0; $gi<@ush_vals; $gi++) {								# calc grid means & stddev
-		my($sum_ush,$sum_vsh,$sum_wsh,$sum_esh);
+		my($sum_ush,$sum_vsh,$sum_wsh,$sum_esh,$sum_lash,$sum_losh);
 
 		$sh_n[$gi] = @{$ush_vals[$gi]};
 		
@@ -421,11 +423,15 @@
 			$sum_vsh += $vsh_vals[$gi][$vi];
 			$sum_wsh += $wsh_vals[$gi][$vi];
 			$sum_esh += $LADCP{ENSEMBLE}[$ens_vals[$gi][$vi]]->{ELAPSED_TIME}+$CTD{first_elapsed}-$opt_l;
+			$sum_lash += $LADCP{ENSEMBLE}[$ens_vals[$gi][$vi]]->{CTD_LAT};
+			$sum_losh += $LADCP{ENSEMBLE}[$ens_vals[$gi][$vi]]->{CTD_LON};
 		}
 		$ush_mu[$gi] = $sh_n[$gi] ? $sum_ush/$sh_n[$gi] : nan;
 		$vsh_mu[$gi] = $sh_n[$gi] ? $sum_vsh/$sh_n[$gi] : nan;
 		$wsh_mu[$gi] = $sh_n[$gi] ? $sum_wsh/$sh_n[$gi] : nan;
 		$esh_mu[$gi] = $sh_n[$gi] ? $sum_esh/$sh_n[$gi] : nan;
+		$lash_mu[$gi] = $sh_n[$gi] ? $sum_lash/$sh_n[$gi] : nan;
+		$losh_mu[$gi] = $sh_n[$gi] ? $sum_losh/$sh_n[$gi] : nan;
 	}
 
 	for (my($gi)=0; $gi<@ush_vals; $gi++) {								# calc & grid stddevs
--- a/LADCPproc.defaults
+++ b/LADCPproc.defaults
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L A D C P P R O C . D E F A U L T S 
 #                    doc: Fri Sep 17 09:44:21 2010
-#                    dlm: Tue Aug  6 20:57:30 2013
+#                    dlm: Thu Sep 19 13:34:34 2013
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 33 66 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 142 75 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # default parameters for [LADCPproc]
@@ -31,6 +31,7 @@
 #	Jun  5, 2013: - added $bad_beam
 #	Sep  6, 2013: - BUG: BT_begin_search_above value of 300m was correct; 
 #						 the original bug was in the documentation
+#	Sep 19, 2013: - added support for $BT_range_method
 
 #----------------------------------------------------------------------
 # Data editing
@@ -117,7 +118,7 @@
 
 # Maximum difference between water depth and average distance of echo max.,
 # if $BT_min_depth and $BT_max_depth are not set. The stddev of the detected
-# water depth is added To this number.
+# water depth is added to this number.
 
 $BT_max_depth_error = 20;
 
@@ -127,6 +128,21 @@
 
 $BT_max_w_difference = 0.03;
 
+# When the BT velocity is taken from the wrong bin, the measurement is not
+# taken from the main acoustic beam but, rather, from a side lobe. As a
+# result, the BT measurement is biased, causing biased near-bottom LADCP
+# velocities when the instrument is moving horizontally.
+# Methods:
+#	0: Take BT velocity from bin with max(Sv). Not recommended.
+#	1: Chose either bin at max(Sv) or one of its neighbors, 
+#	   depending which shows the smallest discrepancy between w_BT and w_reflr. 
+#	   Causes underestimation of instrument speed in P403 tow-yos.
+#	2: Visbeck (2002): use median from 3 bins. Requires $BT_range_Visbeck_center
+#	   to define where the 3 bins should be centered wrt. the max(Sv) bin.
+#	   Visbeck suggests +1. Based on P403 tow-yo 030 it looks more like -2.
+
+$BT_range_method = 1;		
+
 #----------------------------------------------------------------------
 # Shear Processing Parameters
 #----------------------------------------------------------------------
--- a/LADCPproc.loadCTD
+++ b/LADCPproc.loadCTD
@@ -1,9 +1,9 @@
 #======================================================================
 #                    L A D C P P R O C . L O A D C T D 
 #                    doc: Thu Dec  9 18:39:01 2010
-#                    dlm: Tue Jun 25 14:42:50 2013
+#                    dlm: Wed Sep 25 12:59:57 2013
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 25 0 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 210 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -22,6 +22,8 @@
 #				  - BUG: CNV format error was not detected correctly any more
 #   Jan  8, 2013: - added CTD_ASCII_header_lines
 #   Jun 25, 2013: - adapted to :: %PARAM convention
+#	Sep 25, 2013: - renamed "std" %PARAMs %lat, %lon, %ITS to conform to new convention
+#				  - added support for carry-through of lat/lon info
 
 sub readCTD_ASCII($$)
 {
@@ -30,7 +32,7 @@
 	croak("$fn: unknown pressure field\n")    unless defined($CTD_ASCII_press_field);
 	croak("$fn: unknown temperature field\n") unless defined($CTD_ASCII_temp_field);
 	croak("$fn: unknown salinity field\n")    unless defined($CTD_ASCII_salin_field);
-	unless (numberp($dtaR->{lat})) {
+	unless (numberp($dtaR->{stn_lat})) {
 		croak("$fn: unknown latitude field\n")    unless defined($CTD_ASCII_lat_field);
 		croak("$fn: unknown longitude field\n")   unless defined($CTD_ASCII_lon_field);
 	}
@@ -55,16 +57,21 @@
 		if (defined($CTD_ASCII_lat_field) &&
 		    numberp($rec[$CTD_ASCII_lat_field-1]) &&
 			$rec[$CTD_ASCII_lat_field-1] != $CTD_ASCII_badval) {
+				push(@{$dtaR->{lat}},$rec[$CTD_ASCII_lat_field-1]);
+				push(@{$dtaR->{lon}},$rec[$CTD_ASCII_lon_field-1]);
 				$nPos++;
 				$sumLat += $rec[$CTD_ASCII_lat_field-1];
 				$sumLon += $rec[$CTD_ASCII_lon_field-1];
+		} else {
+			push(@{$dtaR->{lat}},nan);
+			push(@{$dtaR->{lon}},nan);
 		}
 	}
 	close(F);
 	
 	if ($nPos > 0) {
-		$dtaR->{lat} = $sumLat / $nPos;
-		$dtaR->{lon} = $sumLon / $nPos;
+		$dtaR->{stn_lat} = $sumLat / $nPos;
+		$dtaR->{stn_lon} = $sumLon / $nPos;
 	}
 
 	$dtaR->{sampint} = 1 / $CTD_ASCII_sampfreq;
@@ -95,6 +102,8 @@
 			$tempF	= $1,next if ($hdr =~ /name (\d+) = t090C:/);
 			$salinF = $1,next if ($hdr =~ /name (\d+) = sal00:/);
 		}
+		$latF = $1,next if ($hdr =~ /name (\d+) = latitude:/);
+		$lonF = $1,next if ($hdr =~ /name (\d+) = longitude:/);
 	
 		&antsAddParams('LADCPproc::CTD_start_time',$1),next				# selected metadata
 			if ($hdr =~ /start_time = (.*)/);
@@ -112,14 +121,14 @@
 	
 		if ($hdr =~ /Latitude\s*[=:]\s*/) {
 			($deg,$min,$NS) = split(/ /,$');
-			$dtaR->{lat} = $deg + $min/60;
-			$dtaR->{lat} *= -1 if ($NS eq 'S');
+			$dtaR->{stn_lat} = $deg + $min/60;
+			$dtaR->{stn_lat} *= -1 if ($NS eq 'S');
 			next;
 		}
 		if ($hdr =~ /Longitude\s*[=:]\s*/) {
 			($deg,$min,$EW) = split(/ /,$');
-			$dtaR->{lon} = $deg + $min/60;
-			$dtaR->{lon} *= -1 if ($EW eq 'W');
+			$dtaR->{stn_lon} = $deg + $min/60;
+			$dtaR->{stn_lon} *= -1 if ($EW eq 'W');
 			next;
 		}
 	    
@@ -155,6 +164,8 @@
 			push(@{$dtaR->{press}},($rec[$pressF] == $CTD_badval) ? nan : $rec[$pressF]);
 			push(@{$dtaR->{temp}}, ($rec[$tempF]  == $CTD_badval) ? nan : $rec[$tempF]);
 			push(@{$dtaR->{salin}},($rec[$salinF] == $CTD_badval) ? nan : $rec[$salinF]);
+			push(@{$dtaR->{press}},(!defined($latF) || ($rec[$latF] == $CTD_badval)) ? nan : $rec[$latF]);
+			push(@{$dtaR->{press}},(!defined($lonF) || ($rec[$lonF] == $CTD_badval)) ? nan : $rec[$lonF]);
 		}
 	} elsif ($CTD_file_type eq 'binary') {
 	
@@ -177,6 +188,8 @@
 			push(@{$dtaR->{press}},($dta[$r*$CTD_nfields+$pressF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$pressF]);
 			push(@{$dtaR->{temp}}, ($dta[$r*$CTD_nfields+$tempF]  == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$tempF]);
 			push(@{$dtaR->{salin}},($dta[$r*$CTD_nfields+$salinF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$salinF]);
+			push(@{$dtaR->{lat}},(!defined($latF) || ($dta[$r*$CTD_nfields+$latF] == $CTD_badval)) ? nan : $dta[$r*$CTD_nfields+$latF]);
+			push(@{$dtaR->{lon}},(!defined($lonF) || ($dta[$r*$CTD_nfields+$lonF] == $CTD_badval)) ? nan : $dta[$r*$CTD_nfields+$lonF]);
 		}
 	} else {
 		croak("$fn: unknown CTD file type $CTD_file_type\n");
@@ -194,13 +207,13 @@
 		readCTD_CNV($fn,$dtaR);
 	}
 
-	croak("$0: unknown latitude\n") unless defined($dtaR->{lat});
-	&antsAddParams('lat',$dtaR->{lat});
-	croak("$0: unknown longitude\n") unless defined($dtaR->{lon});
-	&antsAddParams('lon',$dtaR->{lon});
+	croak("$0: unknown latitude\n") unless defined($dtaR->{stn_lat});
+	&antsAddParams('LADCPproc::CTD_lat',$dtaR->{stn_lat});
+	croak("$0: unknown longitude\n") unless defined($dtaR->{stn_lon});
+	&antsAddParams('LADCPproc::CTD_lon',$dtaR->{stn_lon});
 
 	&antsAddParams('LADCPproc::CTD_sampfreq',1/$dtaR->{sampint});
-	&antsAddParams('ITS',$P{ITS} = 90);
+	&antsAddParams('LADCPproc::CTD_ITS',$P{ITS} = 90);
 }
 
 1;