LADCP_w_ocean
changeset 41 6bddb82924e3
parent 40 408fc95bcff8
child 42 f7690c7b92e0
equal deleted inserted replaced
40:408fc95bcff8 41:6bddb82924e3
     1 #!/usr/bin/perl
     1 #!/usr/bin/perl
     2 #======================================================================
     2 #======================================================================
     3 #                    L A D C P _ W _ O C E A N 
     3 #                    L A D C P _ W _ O C E A N 
     4 #                    doc: Fri Dec 17 18:11:13 2010
     4 #                    doc: Fri Dec 17 18:11:13 2010
     5 #                    dlm: Sun Mar 13 15:08:59 2016
     5 #                    dlm: Thu Mar 17 05:55:25 2016
     6 #                    (c) 2010 A.M. Thurnherr
     6 #                    (c) 2010 A.M. Thurnherr
     7 #                    uE-Info: 473 122 NIL 0 3 72 0 2 4 NIL ofnI
     7 #                    uE-Info: 229 64 NIL 0 0 72 2 2 4 NIL ofnI
     8 #======================================================================
     8 #======================================================================
     9 
     9 
    10 # TODO:
    10 # TODO:
    11 #	! plots:
    11 #	! plots:
    12 #		- make blue labels consistent
       
    13 #		- avoid over-plotting axis labels
    12 #		- avoid over-plotting axis labels
    14 #		- allow for different nsamp magnitudes
    13 #		- allow for different nsamp magnitudes
    15 #		- add "dc" "uc" labels
    14 #		- add "dc" "uc" labels
    16 #		- add software version
    15 #		- add software version
    17 #		- eps.VKE limits
    16 #		- eps.VKE limits
   223 #	Mar  8, 2016: - removed L0 water-depth-difference warning
   222 #	Mar  8, 2016: - removed L0 water-depth-difference warning
   224 #				  - added test for 1500m/s sound speed
   223 #				  - added test for 1500m/s sound speed
   225 #	Mar  9, 2016: - added hab field to .wprof output
   224 #	Mar  9, 2016: - added hab field to .wprof output
   226 #	Mar 13, 2016: - cleaned up warnings created before LADCP_file is defined
   225 #	Mar 13, 2016: - cleaned up warnings created before LADCP_file is defined
   227 #				  - added sanity checks and warnings
   226 #				  - added sanity checks and warnings
       
   227 #	Mar 17, 2016: - added {dc,uc}_rms_{tilt,w_pkg}
       
   228 #				  - replaced a couple of **2 by &SQR()
       
   229 #				  - replaced %ADCP_orientation values by DL & UL
   228 # HISTORY END
   230 # HISTORY END
   229 
   231 
   230 # CTD REQUIREMENTS
   232 # CTD REQUIREMENTS
   231 #	- elapsed		elapsed seconds; see note below
   233 #	- elapsed		elapsed seconds; see note below
   232 #	- depth
   234 #	- depth
   763 	next unless ($per_bin_nsamp[$bin]/($lastGoodEns-$firstGoodEns) < $per_bin_valid_frac_lim);
   765 	next unless ($per_bin_nsamp[$bin]/($lastGoodEns-$firstGoodEns) < $per_bin_valid_frac_lim);
   764 	$first_bad_bin = $bin;
   766 	$first_bad_bin = $bin;
   765 	last;
   767 	last;
   766 }
   768 }
   767 
   769 
       
   770 my($dc_sumsq,$dc_n,$uc_sumsq,$uc_n) = (0,0,0,0);
   768 $fprm = $pte = 0;
   771 $fprm = $pte = 0;
   769 for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
   772 for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
   770 	$pte += editTilt($ens,$opt_t);
   773 	my($ne) = editTilt($ens,$opt_t);
       
   774 	if ($ne == 0) {
       
   775 		if ($ens <= $LADCP_atbottom) {
       
   776 			$dc_sumsq += &SQR($LADCP{ENSEMBLE}[$ens]->{TILT});
       
   777 			$dc_n++;
       
   778         } else {
       
   779 			$uc_sumsq += &SQR($LADCP{ENSEMBLE}[$ens]->{TILT});
       
   780 			$uc_n++;
       
   781         }
       
   782     } else {
       
   783 		$pte += $ne;
       
   784     }
   771 	$fprm += editFarBins($ens,$first_bad_bin) if defined($first_bad_bin);
   785 	$fprm += editFarBins($ens,$first_bad_bin) if defined($first_bad_bin);
   772 }
   786 }
       
   787 my($dc_rms_tilt) = ($dc_n > 0) ? sqrt($dc_sumsq/$dc_n) : nan;
       
   788 my($uc_rms_tilt) = ($uc_n > 0) ? sqrt($uc_sumsq/$uc_n) : nan;
       
   789 &antsAddParams('dc_rms_tilt',$dc_rms_tilt,'uc_rms_tilt',$uc_rms_tilt);
       
   790 
   773 progress("\tattitude threshold (max_tilt = %d deg): %d velocites removed (%d%% of total)\n",
   791 progress("\tattitude threshold (max_tilt = %d deg): %d velocites removed (%d%% of total)\n",
   774 	$opt_t,$pte,round(100*$pte/$nvv));
   792 	$opt_t,$pte,round(100*$pte/$nvv));
   775 progress("\tvelocities beyond bin $first_bad_bin (<%d%% valid values): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
   793 progress("\tvelocities beyond bin $first_bad_bin (<%d%% valid values): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",
   776 	round(100*$per_bin_valid_frac_lim),$fprm,round(100*$fprm/$nvw));
   794 	round(100*$per_bin_valid_frac_lim),$fprm,round(100*$fprm/$nvw));
   777 
   795 
   812 
   830 
   813 $CTD{W_t}[0]  = $CTD{W_t} [@{$CTD{ELAPSED}}] = nan;							# calculate w-derivatives
   831 $CTD{W_t}[0]  = $CTD{W_t} [@{$CTD{ELAPSED}}] = nan;							# calculate w-derivatives
   814 $CTD{W_tt}[0] = $CTD{W_tt}[@{$CTD{ELAPSED}}] = nan;
   832 $CTD{W_tt}[0] = $CTD{W_tt}[@{$CTD{ELAPSED}}] = nan;
   815 for (my($s)=1; $s<@{$CTD{ELAPSED}}-1; $s++) {								# centered differences for both
   833 for (my($s)=1; $s<@{$CTD{ELAPSED}}-1; $s++) {								# centered differences for both
   816 	$CTD{W_t} [$s] = ($CTD{W}[$s+1] - $CTD{W}[$s-1]) / (2*$CTD{DT});
   834 	$CTD{W_t} [$s] = ($CTD{W}[$s+1] - $CTD{W}[$s-1]) / (2*$CTD{DT});
   817 	$CTD{W_tt}[$s] = ($CTD{W}[$s+1] + $CTD{W}[$s-1] - 2*$CTD{W}[$s]) / $CTD{DT}**2;
   835 	$CTD{W_tt}[$s] = ($CTD{W}[$s+1] + $CTD{W}[$s-1] - 2*$CTD{W}[$s]) / &SQR($CTD{DT});
   818 }
   836 }
   819 
   837 
   820 error("$0: CTD start depth must be numeric\n")
   838 error("$0: CTD start depth must be numeric\n")
   821 	unless numberp($CTD{DEPTH}[0]);
   839 	unless numberp($CTD{DEPTH}[0]);
   822 
   840 
   956 &antsAddParams('w_max_lim',$w_max_lim);
   974 &antsAddParams('w_max_lim',$w_max_lim);
   957 
   975 
   958 my($cli) = 0;																	# current-lag index
   976 my($cli) = 0;																	# current-lag index
   959 my($lag) = $CTD_time_lag[$cli];													# current lag
   977 my($lag) = $CTD_time_lag[$cli];													# current lag
   960 
   978 
       
   979 my($dc_sumsq,$dc_n,$uc_sumsq,$uc_n) = (0,0,0,0);								# for w_pkg calcs
   961 for (my($skipped)=0,my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
   980 for (my($skipped)=0,my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) {
   962 
   981 
   963 	if ($ens > $CTD_tl_toEns[$cli]) {											# use correct lag piece
   982 	if ($ens > $CTD_tl_toEns[$cli]) {											# use correct lag piece
   964 		$cli++;
   983 		$cli++;
   965 		die("assertion failed!\n\ttest: \$cli(=$cli) <= \$#CTD_time_lag\n")
   984 		die("assertion failed!\n\ttest: \$cli(=$cli) <= \$#CTD_time_lag\n")
   997 		error(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs\n",
  1016 		error(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs\n",
   998 			$CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan]))
  1017 			$CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan]))
   999 				unless ($CTD{DEPTH}[$scan] >= 0);
  1018 				unless ($CTD{DEPTH}[$scan] >= 0);
  1000 		$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan];
  1019 		$LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan];
  1001 		$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan;
  1020 		$LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan;
       
  1021 		if ($ens <= $LADCP_atbottom) {
       
  1022 			$dc_sumsq += &SQR($CTD{W}[$scan]); $dc_n++;
       
  1023         } else {
       
  1024 			$uc_sumsq += &SQR($CTD{W}[$scan]); $uc_n++;
       
  1025 		}
  1002 		my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan];
  1026 		my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan];
  1003 		if (abs($reflr_ocean_w) <= $w_max_lim) {
  1027 		if (abs($reflr_ocean_w) <= $w_max_lim) {
  1004 			$sumWsq += &SQR($reflr_ocean_w);
  1028 			$sumWsq += &SQR($reflr_ocean_w);
  1005 			$nWsq++;
  1029 			$nWsq++;
  1006 			if ($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} > 100 &&
  1030 			if ($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} > 100 &&
  1014 	} else{
  1038 	} else{
  1015 	    undef($LADCP{ENSEMBLE}[$ens]->{REFLR_W});						# don't output in time-series file
  1039 	    undef($LADCP{ENSEMBLE}[$ens]->{REFLR_W});						# don't output in time-series file
  1016 	    undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});						# old DEPTH from calcLADCPts() must be removed
  1040 	    undef($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});						# old DEPTH from calcLADCPts() must be removed
  1017 	}
  1041 	}
  1018 }
  1042 }
       
  1043 
       
  1044 my($dc_rms_wpkg) = ($dc_n > 0) ? sqrt($dc_sumsq / $dc_n) : nan;			# rms package velocity
       
  1045 my($uc_rms_wpkg) = ($uc_n > 0) ? sqrt($uc_sumsq / $uc_n) : nan;
       
  1046 &antsAddParams('dc_rms_w_pkg',$dc_rms_wpkg,'uc_rms_w_pkg',$uc_rms_wpkg);
  1019 	
  1047 	
  1020 if ($nWsq > 0 && $nWsqI > 0) {
  1048 if ($nWsq > 0 && $nWsqI > 0) {
  1021 	&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI));
  1049 	&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI));
  1022 	progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n",
  1050 	progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n",
  1023 						100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI));
  1051 						100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI));
  1109 	warning(2,"unknown water depth --- cannot edit sidelobes or PPI near the seabed\n");
  1137 	warning(2,"unknown water depth --- cannot edit sidelobes or PPI near the seabed\n");
  1110 	&antsAddParams('water_depth','unknown','water_depth.sig','nan');
  1138 	&antsAddParams('water_depth','unknown','water_depth.sig','nan');
  1111 }
  1139 }
  1112 	
  1140 	
  1113 if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {				# DOWNLOOKER
  1141 if ($LADCP{ENSEMBLE}[$LADCP_atbottom]->{XDUCER_FACING_DOWN}) {				# DOWNLOOKER
  1114 	&antsAddParams('ADCP_orientation','downlooker');
  1142 	&antsAddParams('ADCP_orientation','DL');
  1115 
  1143 
  1116 	if (defined($water_depth)) {
  1144 	if (defined($water_depth)) {
  1117 		progress("Editing data to remove sidelobe interference from seabed...\n");
  1145 		progress("Editing data to remove sidelobe interference from seabed...\n");
  1118 		($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,$water_depth);
  1146 		($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,$water_depth);
  1119 		progress("\t$nvrm velocities from $nerm ensembles removed\n");
  1147 		progress("\t$nvrm velocities from $nerm ensembles removed\n");
  1147 	        progress("\t$nvrm velocities from $nerm ensembles removed\n");
  1175 	        progress("\t$nvrm velocities from $nerm ensembles removed\n");
  1148 	    }
  1176 	    }
  1149 	}
  1177 	}
  1150 	
  1178 	
  1151 } else {																	# UPLOOKER
  1179 } else {																	# UPLOOKER
  1152 	&antsAddParams('ADCP_orientation','uplooker');
  1180 	&antsAddParams('ADCP_orientation','UL');
  1153 
  1181 
  1154 	progress("Editing data to remove sidelobe interference from sea surface...\n");
  1182 	progress("Editing data to remove sidelobe interference from sea surface...\n");
  1155 	($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,undef);
  1183 	($nvrm,$nerm) = editSideLobes($firstGoodEns,$lastGoodEns,undef);
  1156 	progress("\t$nvrm velocities from $nerm ensembles removed\n");
  1184 	progress("\t$nvrm velocities from $nerm ensembles removed\n");
  1157 
  1185