LADCP_VKE
changeset 56 8f120b9f795a
parent 54 828e5466391b
child 57 69e39fcb7f41
equal deleted inserted replaced
55:2d8e1139acd5 56:8f120b9f795a
     1 #!/usr/bin/perl
     1 #!/usr/bin/perl
     2 #======================================================================
     2 #======================================================================
     3 #                    L A D C P _ V K E 
     3 #                    L A D C P _ V K E 
     4 #                    doc: Tue Oct 14 11:05:16 2014 
     4 #                    doc: Tue Oct 14 11:05:16 2014 
     5 #                    dlm: Thu Sep 12 13:49:04 2019
     5 #                    dlm: Sat Jul 24 09:44:26 2021
     6 #                    (c) 2012 A.M. Thurnherr
     6 #                    (c) 2012 A.M. Thurnherr
     7 #                    uE-Info: 111 71 NIL 0 0 72 0 2 4 NIL ofnI
     7 #                    uE-Info: 120 45 NIL 0 0 72 0 2 4 NIL ofnI
     8 #======================================================================
     8 #======================================================================
     9 
     9 
    10 $antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
    10 $antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
    11 
    11 
    12 # NOTES:
    12 # NOTES:
   107 #				  - added eps.ms field to files processed without -k
   107 #				  - added eps.ms field to files processed without -k
   108 #	Dec  9, 2017: - added support for $antsSuppressCommonOptions
   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
   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
   110 #	Apr 25, 2018: - added -y and removed spectral bins from default output
   111 #	Sep 12, 2019: - disabled default -l cut-off (used to be 5e-11 W/kg)
   111 #	Sep 12, 2019: - disabled default -l cut-off (used to be 5e-11 W/kg)
       
   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
       
   121 # HISTORY END
   112 
   122 
   113 ($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
   123 ($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
   114 ($WCALC)   = ($0              =~ m{^(.*)/[^/]*$});
   124 ($WCALC)   = ($0              =~ m{^(.*)/[^/]*$});
   115 $WCALC = '.' if ($WCALC eq '');
   125 $WCALC = '.' if ($WCALC eq '');
   116 $ANTS_TOOLS_AVAILABLE = (`which list` ne '');
   126 $ANTS_TOOLS_AVAILABLE = (`which list` ne '');
   127 #----------------------------------------------------------------------
   137 #----------------------------------------------------------------------
   128 # Empirical constants & defaults
   138 # Empirical constants & defaults
   129 #----------------------------------------------------------------------
   139 #----------------------------------------------------------------------
   130 
   140 
   131 #my($c) = 0.0215;						# Thurnherr et al. (GRL 2015)
   141 #my($c) = 0.0215;						# Thurnherr et al. (GRL 2015)
   132 my($c) = 0.026;							# increased by 20% for V1.2beta7
   142 my($c) = 0.026;							# increased by 21% for V1.2beta7
   133 $opt_q = 3;								# Equatorial band: little more than a guess based on 2015 P16N
   143 $opt_q = 3;								# Equatorial band: little more than a guess based on 2015 P16N
   134 $opt_l = 0;								# [W/kg]; cutoff disabled Sep 12, 2019
   144 $opt_l = 0;								# [W/kg]; cutoff disabled Sep 12, 2019
   135 $opt_a = nan;							# assume background dissipation for samples that pass the tests but have eps below -l
   145 $opt_a = 0;								# assume background dissipation for samples that pass the tests but have eps below -l
   136 $opt_z = 1;								# number of w_ocean samples to require
   146 $opt_z = 50;							# number of w_ocean samples to require (note that the .wprof inputs may have harsher limits)
   137 $opt_o = 0;								# remove mean before calculating spectra
   147 $opt_o = 0;								# remove mean before calculating spectra
   138 $opt_s = 150;							# surface layer to exclude from spectra
   148 $opt_s = 150;							# surface layer to exclude from spectra
   139 $opt_g = 40;							# max gap to interpolate over
   149 $opt_g = 40;							# max gap to interpolate over
   140 $opt_c_default = 100;					# short-wavelength cutoff
   150 $opt_c_default = 100;					# short-wavelength cutoff
   141 
   151 
   194 	for (my($bufi)=0; defined($ARGV[0]); $bufi++) {									# multiple input files: loop until 2nd last
   204 	for (my($bufi)=0; defined($ARGV[0]); $bufi++) {									# multiple input files: loop until 2nd last
   195 		$input_list .= "$P{profile_id}($P{run_label}) ";
   205 		$input_list .= "$P{profile_id}($P{run_label}) ";
   196 		if ($bufi == 0) {															# do once for mulitple files
   206 		if ($bufi == 0) {															# do once for mulitple files
   197 			&antsAddParams('ADCP_bin_length','','ADCP_blanking_distance','',		# delete most %PARAMs, leaving
   207 			&antsAddParams('ADCP_bin_length','','ADCP_blanking_distance','',		# delete most %PARAMs, leaving
   198 						   'ADCP_frequency','','ADCP_orientation','',				# 	only potentially useful ones (from last file):
   208 						   'ADCP_frequency','','ADCP_orientation','',				# 	only potentially useful ones (from last file):
   199 						   'ADCP_pulse_length','','BT_max_bin_range_diff','',		#		%profile_id
   209 						   'ADCP_pulse_length','','BT_bin_range_diff.max','',		#		%profile_id
   200 						   'BT_max_range','','BT_max_w_error','',					#		%water_depth
   210 						   'BT_range.max','','BT_w_error.max','',					#		%water_depth
   201 						   'BT_rms_w_discrepancy','','CTD_time_lags','',			#		%lat
   211 						   'BT_w_discrepancy.rms','','CTD_time_lags','',			#		%lat
   202 						   'LADCP_firstBin','','LADCP_lastBin','',					#		%lon
   212 						   'LADCP_firstBin','','LADCP_lastBin','',					#		%lon
   203 						   'LADCP_w_ocean::version','','SS_min_samp','',			#		%dnXX
   213 						   'LADCP_w_ocean::version','','SS_samp.min','',			#		%dnXX
   204 						   'SS_max_allowed_depth_range','','SS_min_signal','',
   214 						   'SS_allowed_depth_range.max','','SS_signal.min','',
   205 						   'Sv_ref_bin','','TL_max_allowed_three_lag_spread','',
   215 						   'Sv_ref_bin','','TL_allowed_three_lag_spread.max','',
   206 						   'dc_rms_accel_pkg','','dc_rms_tilt','',
   216 						   'dc_pkg_accel.rms','','dc_tilt.rms','',
   207 						   'uc_rms_accel_pkg','','uc_rms_tilt','',
   217 						   'uc_pkg_accel.rms','','uc_tilt.rms','',
   208 						   'max_depth','','max_elapsed','','max_ens','',
   218 						   'depth.max','','elapsed.max','','ens.max','',
   209 						   'min_depth','','min_elapsed','','min_ens','',
   219 						   'depth.min','','elapsed.min','','ens.min','',
   210 						   'out_basename','','outgrid_dz','','run_label','',
   220 						   'out_basename','','outgrid_dz','','run_label','',
   211 						   'outgrid_firstbin','','outgrid_lastbin','',
   221 						   'outgrid_firstbin','','outgrid_lastbin','',
   212 						   'outgrid_minsamp','','per_bin_valid_frac_lim','',
   222 						   'outgrid_minsamp','','per_bin_valid_frac_lim','',
   213 						   'processing_options','','refLr_firstBin','',
   223 						   'processing_options','','refLr_firstBin','',
   214 						   'refLr_lastBin','','rms_w_reflr_err','',
   224 						   'refLr_lastBin','','rms_w_reflr_err','',
   215 						   'rms_w_reflr_err_interior','',
   225 						   'rms_w_reflr_err_interior','',
   216 						   'sidelobe_editing','','surface_layer_depth','',
   226 						   'sidelobe_editing','','surface_layer_depth','',
   217 						   'vessel_draft','','w_max_lim','',
   227 						   'vessel_draft','','w_lim.max','',
   218 						   'water_depth.sig','','water_depth_from','',
   228 						   'water_depth.sig','','water_depth_from','',
   219 			);
   229 			);
   220 		}
   230 		}
   221 		close(TOCLD);																# close LADCP_VKE input
   231 		close(TOCLD);																# close LADCP_VKE input
   222     	my(@specrec);
   232     	my(@specrec);
   265 #	- spectra from previous files are in @specbuf
   275 #	- spectra from previous files are in @specbuf
   266 #----------------------------------------------------------------------
   276 #----------------------------------------------------------------------
   267 
   277 
   268 $n_input_files = 1 + @specbuf;									# number of input files provided
   278 $n_input_files = 1 + @specbuf;									# number of input files provided
   269 
   279 
       
   280 &antsFloatOpt(\$opt_e,$c);										# default parameterization
       
   281 &antsFloatOpt(\$opt_x,1);										# spectral fit stddev scale factor
       
   282 
   270 &antsAddParams('LADCP_VKE::input_files.n',$n_input_files,
   283 &antsAddParams('LADCP_VKE::input_files.n',$n_input_files,
   271 			   'LADCP_VKE::wsamp.min',$opt_z,
   284 			   'LADCP_VKE::wsamp.min',$opt_z,
   272 			   'LADCP_VKE::eps.minlim',$opt_l);
   285 			   'LADCP_VKE::eps.minlim',$opt_l,
   273 
   286 			   'LADCP_VKE::calibration_constant',$opt_e);
   274 &antsFloatOpt(\$opt_e,$c);										# default parameterization
       
   275 &antsFloatOpt(\$opt_x,1);										# spectral fit stddev scale factor
       
   276 
   287 
   277 if (defined($opt_c)) {											# shortwave cutoff supplied
   288 if (defined($opt_c)) {											# shortwave cutoff supplied
   278 	$lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c;
   289 	$lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c;
   279 	&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin);	# ensure eps.VKE is calculated below
   290 	&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin);	# ensure eps.VKE is calculated below
   280 } elsif (defined(antsParam('shortwave_cutoff'))) {				# cutoff already applied
   291 } elsif (defined(antsParam('shortwave_cutoff'))) {				# cutoff already applied
   361 		}
   372 		}
   362 		my($wi) = ($ants_[$r][$widxf]>0) ? $ants_[$r][$widxf] : 0;				# adjust for -1 
   373 		my($wi) = ($ants_[$r][$widxf]>0) ? $ants_[$r][$widxf] : 0;				# adjust for -1 
   363 		for (my($bi)=0; $bi<@specbuf; $bi++) {									# loop over all buffered files
   374 		for (my($bi)=0; $bi<@specbuf; $bi++) {									# loop over all buffered files
   364 			next unless @{$specbuf[$bi][$wi]};									# skip input files w/o valid spectra
   375 			next unless @{$specbuf[$bi][$wi]};									# skip input files w/o valid spectra
   365 			if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) > 0) {			# depth mismatch
   376 			if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) > 0) {			# depth mismatch
   366 				die("assertion failed") unless ($wi == 0);						# only allowed in bottom window
   377 #				die("assertion failed ($specbuf[$bi][$wi][$df] - $ants_[$r][$df] @ wi=$wi)") unless ($wi == 0);						# only allowed in bottom window
   367 				if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) >
   378 				if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) >
   368 					abs($ants_[$r][$maxdf] - $ants_[$r][$mindf])) {
   379 					abs($ants_[$r][$maxdf] - $ants_[$r][$mindf])) {
   369 						printf(STDERR "WARNING: ignoring bottom window from input file #%d because of depth mismatch\n",$bi+1)
   380 						printf(STDERR "WARNING: ignoring window #$wi from input file #%d because of depth mismatch\n",$bi+1)
   370 							if ($f == $fs_fmin);
   381 							if ($f == $fs_fmin);
   371 						next;
   382 						next;
   372 				}
   383 				}
   373 			}
   384 			}
   374 #			print(STDERR "specbuf[$bi][$wi][$f] = $specbuf[$bi][$wi][$f]\n");
       
   375 			if (numberp($specbuf[$bi][$wi][$f])) {
   385 			if (numberp($specbuf[$bi][$wi][$f])) {
   376 				$sum += $specbuf[$bi][$wi][$f]; $ns++;
   386 				$sum += $specbuf[$bi][$wi][$f]; $ns++;
   377             }
   387             }
   378 			$ants_[$r][$doff] += $specbuf[$bi][$wi][$doff]						# update nspec once per input record
   388 			$ants_[$r][$doff] += $specbuf[$bi][$wi][$doff]						# update nspec once per input record
   379 				if ($f == $fs_fmax);											#	... but for all files
   389 				if ($f == $fs_fmax);											#	... but for all files
   560 		$ants_[$r][$fspwrf] = $ants_[$r][$p0f] = $ants_[$r][$rmsf] =
   570 		$ants_[$r][$fspwrf] = $ants_[$r][$p0f] = $ants_[$r][$rmsf] =
   561 			$ants_[$r][$rf] = $ants_[$r][$slpf] = $ants_[$r][$sslpf] = nan;
   571 			$ants_[$r][$rf] = $ants_[$r][$slpf] = $ants_[$r][$sslpf] = nan;
   562 	}
   572 	}
   563 
   573 
   564 	#-----------------------------------------------------------------------------------------------------
   574 	#-----------------------------------------------------------------------------------------------------
   565 	# QC Tests:
   575 	# eps.VKE QC Tests:
   566 	#	- the following limits were independently derived 
   576 	#	- the following limits were independently derived 
   567 	#		p0fit.rms <= 0.4			primary filter used in Thurnherr et al. (GRL 2015)
   577 	#		p0fit.rms <= 0.4			primary filter used in Thurnherr et al. (GRL 2015)
   568 	#		-3 <= p0fit.slope <= -1		based largely on 2016 I08S data with sufficient/insufficient range
   578 	#		-3 <= p0fit.slope <= -1		based largely on 2016 I08S data with sufficient/insufficient range
   569 	#		p0fit.r <= -0.5				based largely on 2016 I08S data with sufficient/insufficient range
   579 	#		p0fit.r <= -0.5				based largely on 2016 I08S data with sufficient/insufficient range
   570 	#		w.nsamp.avg >= 50			based on observations in many data sets with weak backscatter, 
   580 	#		w.nsamp.avg >= 50			based on observations in many data sets with weak backscatter,