function [ds,p,dr] = calc_shear3(d,p,ps,dr)
% function [ds,p,dr] = calc_shear3(d,p,ps,dr)
%
% - compute shear profiles
% - use only central difference
% - use 2*std editing
%
% version 8 last change 06.08.2019
% Martin Visbeck, LDEO, 3/7/97
% some modifications and code cleanup GK, 16.05.2015 2-->3
% reinstated the use of d.weight for the identification of shear pairs
% GK, 03.12.2018 3-->4
% reorganization of the dz handling GK, 12.12.2018 4-->5
% added error output, catch different shear methods GK, 22.07.2019 5-->6
% more text output
% added new parameter ps.shear_throw_out_percent GK, 05.08.2019 6-->7
% fixed bug in shear calculation GK, 06.08.2019 7-->8
%======================================================================
% C A L C _ S H E A R 3 . M
% doc: Wed Sep 4 17:55:22 2019
% dlm: Wed Sep 4 18:55:36 2019
% (c) 2019 G. Krahmann, M. Visbeck
% uE-Info: 400 0 NIL 0 0 72 0 2 4 NIL ofnI
%======================================================================
% CHANGES BY ANT:
% Sep 4, 2019: - replaced messages by LDEO warn mechanism
%
% general function info
%
disp(' ')
disp(['CALC_SHEAR3: calculate a baroclinic velocity profile based on shears only'])
%
% decide whether central differences or simple differences are to be used in
% the shear calculation
%
% 1: simple differences
% 2: central differences
%
diff_type = 2;
%
% resolution of final shear profile in meter
%
dz = ps.dz;
shear_dz = dz * diff_type;
disp([' Averaging shear profile over ',num2str(shear_dz),' m intervals'])
%
% Discard a certain amount of data as suspected outliers that
% for relatively small numbers of values will skew the mean.
% In the old calc_shear2.m the allowed range was the std of the shears
% times the set factor stdf. Default was stdf=2 .
% In calc_shear3.m the outlier determination is iterative and thus a bit
% safer in calculation. The allowed range is set as a fraction of the whole
% population. For gaussian distributions stdf converts to this fraction.
% As usually stdf is not varied much, we have implemented a lookup list
% here.
%
use_new_outlier = 1;
ps.shear_stdf = 4; % arbitrary guess
stdf = ps.shear_stdf;
if use_new_outlier==1
fracs = 1 - [31.8,13.4,4.5,1.3,0.3]/100;
stdfs = [1,1.5,2,2.5,3];
[dummy,ind] = min(abs(stdf-stdfs));
frac = fracs(ind);
end
disp([' Maximum allowed std within calculation intervals : ',num2str(stdf)])
disp([' Data deviating more from the median will be discarded.'])
%
% check if only one istrument is to be used
%
if ps.up_dn_looker==2
% down looker only
d.weight(d.izu,:)=nan;
elseif ps.up_dn_looker==3
% up looker only
d.weight(d.izd,:)=nan;
end
%
% Apply a weight threshold to shear data.
%
% There are two variants.
% First one is used when ps.shear_throw_out_percent is not NaN.
% This one throws out the xx% shears with the lowest correlation-derived weights.
% Second one is used when ps.shear_weightmin is not NaN.
% This one sets a minimum weight to keep the shears.
%
% First one is now default with 10%.
%
disp([' Correlation-derived weights range from ',num2str(min(d.weight(:))),' to ',num2str(max(d.weight(:)))])
if 0 % ~isnan(ps.shear_throw_out_percent)
local_weight = d.weight;
ind = find(isfinite(local_weight));
[dummy,ind2] = sort(local_weight(ind));
if length(ind2)>9
local_weight(ind(ind2(1:floor(length(ind2)/10)))) = nan;
end
elseif 0 % ~isnan(ps.shear_weightmin)
disp([' Minimum weight ',num2str(ps.shear_weightmin),' for data to be used in shear calc.'])
local_weight = double(d.weight>ps.shear_weightmin);
else
disp('> No weight criterion applied to raw shear data.')
disp('> You should set ps.shear_throw_out_percent or ps.shear_weightmin')
end
local_weight = d.weight;
disp([' Removed ',int2str(100-sum(isfinite(local_weight))/sum(isfinite(d.weight))*100),...
' % of data with lowest weights from shear calculation.'])
disp([' New weights range from ',num2str(min(local_weight(:))),' to ',num2str(max(local_weight(:)))])
%
% convert the weights to 1 and NaN
%
%local_weight = replace( local_weight, local_weight<=0, nan );
%local_weight = replace( local_weight, local_weight>0, 1 );
local_weight(find(local_weight <= 0)) = nan;
local_weight(find(local_weight > 0)) = 1;
%
% compute shear
%
% two ways are offered here
% first: central differences for the shears
% second: single differences
% the first is similar to the ways of the old calc_shear2.m
%
% central differences
if diff_type==2
local_weight = [repmat(nan,1,size(local_weight,2));diff2(local_weight)+1;repmat(nan,1,size(local_weight,2))];
ushear = [NaN*d.ru(1,:);diff2(d.ru(:,:))./diff2(d.izm);NaN*d.ru(1,:)].*local_weight;
vshear = [NaN*d.rv(1,:);diff2(d.rv(:,:))./diff2(d.izm);NaN*d.rv(1,:)].*local_weight;
wshear = [NaN*d.rw(1,:);diff2(d.rw(:,:))./diff2(d.izm);NaN*d.rw(1,:)].*local_weight;
zshear = -d.izm;
% single differences
else
ushear = diff( d.ru.*local_weight )./diff(d.izm);
vshear = diff( d.rv.*local_weight )./diff(d.izm);
wshear = diff( d.rw.*local_weight )./diff(d.izm);
zshear = -(d.izm(1:end-1,:)+d.izm(2:end,:))/2;
end
ds.ushear = ushear;
ds.vshear = vshear;
ds.wshear = wshear;
ds.zshear = zshear;
%
% set depth levels
%
z = dr.z;
%
% prepare shear solution result vectors
%
ds.usm = repmat(nan,length(z),1);
ds.vsm = ds.usm;
ds.wsm = ds.usm;
ds.usmd = ds.usm;
ds.vsmd = ds.usm;
ds.use = ds.usm;
ds.vse = ds.usm;
ds.wse = ds.usm;
ds.nn = ds.usm;
ds.z = z;
%
% loop over depth levels and calculate the average shear at that level
%
% in the case of central differences this is oversampled here
% but by sticking with the same resolution it makes the results easier
% to work with
%
for n=[1:length(z)]
i1 = find( ( abs( zshear - z(n) ) <= shear_dz/2 ) & isfinite( ushear + vshear ) );
ds.nn(n) = length(i1);
if ds.nn(n) > 2
% two ways to select outliers
% first: select all that are beyond a fixed range around the median
% second: iteratively reject the worst (largest distance from mean)
% until a fixed fraction is rejected
% the second is usually the safer calculation but is a bit slower
if 1
usmm = median( ushear(i1) );
ussd1 = std( ushear(i1) );
vsmm = median( vshear(i1) );
vssd1 = std( vshear(i1) );
wsmm = median( wshear(i1) );
wssd1 = std( wshear(i1) );
ii1 = i1( find(abs(ushear(i1)-usmm)<stdf*ussd1) );
ii2 = i1( find(abs(vshear(i1)-vsmm)<stdf*vssd1) );
ii3 = i1( find(abs(wshear(i1)-wsmm)<stdf*wssd1) );
else
[dummy,ii1] = meanoutlier(ushear(i1),frac);
[dummy,ii2] = meanoutlier(vshear(i1),frac);
[dummy,ii3] = meanoutlier(wshear(i1),frac);
ii1 = i1(ii1);
ii2 = i1(ii2);
ii3 = i1(ii3);
end
% two ways of calculating the mean and std of the selected shears
% first: if there is a rejected one in any of u,v,w shears then use it
% for non of the calculations
% second: if there is a rejected one in any of u,v,w shears then use it
% only in u,v or w calculations
% the second one is the one used by the old calc_shear2.m
% but to me this does not make sense, GK May 2015
if 1
dummy = zeros(size(ushear));
dummy(ii1) = 1;
dummy(ii2) = dummy(ii2)+1;
dummy(ii3) = dummy(ii3)+1;
ii = find(dummy==3);
if length(ii)>1
ds.usm(n) = mean(ushear(ii));
ds.usmd(n) = median(ushear(ii));
ds.use(n) = std(ushear(ii));
ds.ii(n) = length(ii);
ds.vsm(n) = mean(vshear(ii));
ds.vsmd(n) = median(vshear(ii));
ds.vse(n) = std(vshear(ii));
ds.wsm(n) = mean(wshear(ii));
ds.wsmd(n) = median(wshear(ii));
ds.wse(n) = std(wshear(ii));
% debugging plot
if 0
figure(3)
clf
subplot(3,1,1)
hist(ushear(ii),30)
hold on
ax = axis;
plot([1,1]*mean(ushear(ii)),ax(3:4),'r')
plot([1,1]*median(ushear(ii)),ax(3:4),'m')
ind = find(dr.z==z(n));
if ind<length(dr.u)
plot([1,1]*(dr.u(ind-1)-dr.u(ind+1))/20,ax(3:4),'g')
end
title(int2str(z(n)))
subplot(3,1,2)
hist(vshear(ii),30)
hold on
ax = axis;
plot([1,1]*mean(vshear(ii)),ax(3:4),'r')
plot([1,1]*median(vshear(ii)),ax(3:4),'m')
ind = find(dr.z==z(n));
if ind<length(dr.u)
plot([1,1]*(dr.v(ind-1)-dr.v(ind+1))/20,ax(3:4),'g')
end
subplot(3,1,3)
hist( zshear(ii)-z(n), 30 )
pause
end
end
else
if length(ii1)>1
ds.usm(n) = mean(ushear(ii1));
ds.usmd(n) = median(ushear(ii1));
ds.use(n) = std(ushear(ii1));
end
if length(ii2)>1
ds.vsm(n) = mean(vshear(ii2));
ds.vsmd(n) = median(vshear(ii2));
ds.vse(n) = std(vshear(ii2));
end
if length(ii3)>1
ds.wsm(n) = mean(wshear(ii3));
ds.wsmd(n) = median(wshear(ii3));
ds.wse(n) = nstd(wshear(ii3));
end
end
end
end
%
% a debugging figure
%
if 0
sfigure(3);
clf
orient tall
subplot(1,2,1)
plot(ushear,zshear,'b.','markersize',3)
hold on
plot(ds.usm,ds.z,'r')
plot(ds.usmd,ds.z,'k')
inv_shear_u = -diff(dr.u-mean(dr.u))/dz;
plot(inv_shear_u,(z(1:end-1)+z(2:end))/2,'g')
set(gca,'ydir','reverse')
subplot(1,2,2)
plot(vshear,zshear,'b.','markersize',3)
hold on
plot(ds.vsm,ds.z,'r')
plot(ds.vsmd,ds.z,'k')
inv_shear_v = -diff(dr.v-mean(dr.v))/dz;
plot(inv_shear_v,(z(1:end-1)+z(2:end))/2,'g')
set(gca,'ydir','reverse')
sfigure(2)
end
%
% integrate shear profile (from bottom up)
%
%%ds.usm = replace( ds.usm, isnan(ds.usm), 0 );
%%ds.vsm = replace( ds.vsm, isnan(ds.vsm), 0 );
%%ds.wsm = replace( ds.wsm, isnan(ds.wsm), 0 );
ds.usm(find(isnan(ds.usm))) = 0;
ds.vsm(find(isnan(ds.vsm))) = 0;
ds.wsm(find(isnan(ds.wsm))) = 0;
if length(ind)/length(ds.usm)*100>5
disp(['> Found ',num2str(length(ind)/length(ds.usm)*100),...
'% Nan in shear data. Integration result might be problematic.'])
end
if 1
ds.ur = flipud(cumsum(flipud(ds.usm)))*dz;
ds.vr = flipud(cumsum(flipud(ds.vsm)))*dz;
ds.wr = flipud(cumsum(flipud(ds.wsm)))*dz;
else
ds.ur = flipud(cumsum(flipud(ds.usmd)))*dz;
ds.vr = flipud(cumsum(flipud(ds.vsmd)))*dz;
ds.wr = flipud(cumsum(flipud(ds.wsmd)))*dz;
end
ds.ur = ds.ur-mean(ds.ur);
ds.vr = ds.vr-mean(ds.vr);
ds.wr = ds.wr-mean(ds.wr);
%
% This is a peculiar place for the single ping error estimate. But
% as it is based on the variability in the data itself, it makes sense.
% The assumption is that there should be basically zero shear in the
% vertical velocities. At least it is so small as to be not detectable
% here. Thus any variability in the vertical shear is caused by the
% errors/noise of the measurement. Together with an angular conversion
% factor this gives an error/noise value for the horizontal velocities.
%
if ~isfield(d,'zd')
dz2 = abs(mean(diff(d.z)));
else
dz2 = diff_type*abs(mean(diff(d.zd)));
end
if isfield(d,'down')
fac = 1/tan(d.down.Beam_angle*pi/180)*sqrt(2)*dz2;
else
fac = 1/tan(d.up.Beam_angle*pi/180)*sqrt(2)*dz2;
end
ds.ensemble_vel_err = ds.wse*fac;
dr.ensemble_vel_err = ds.wse*fac;
%
% store result and give text output
%
dr.u_shear_method = ds.ur;
dr.v_shear_method = ds.vr;
dr.w_shear_method = ds.wr;
uds = std(dr.u-mean(dr.u)-ds.ur);
vds = std(dr.v-mean(dr.v)-ds.vr);
uvds = sqrt(uds^2+vds^2);
disp([' Inversion average error : ',num2str(mean( dr.uerr ) ),' m/s'])
if uvds>mean(dr.uerr)*1.5
error_increase_factor = 1/mean(dr.uerr)*uvds/1.5;
warn = ('> Increasing error estimate because of elevated shear - inverse difference');
disp(warn)
disp(['> by a factor of ',num2str(error_increase_factor)])
disp(['> std of difference between regular and shear profile : ',num2str(uvds),' m/s'])
p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
dr.uerr = dr.uerr * error_increase_factor;
end
disp([' Final average error : ',num2str(mean( dr.uerr ) ),' m/s'])
%--------------------------------------------------
function x = diff2(x,k,dn)
%DIFF2 Difference function. If X is a vector [x(1) x(2) ... x(n)],
% then DIFF(X) returns a vector of central differences between
% elements [x(3)-x(1) x(4)-x(2) ... x(n)-x(n-2)]. If X is a
% matrix, the differences are calculated down each column:
% DIFF(X) = X(3:n,:) - X(1:n-2,:).
% DIFF(X,n) is the n'th difference function.
% J.N. Little 8-30-85
% Copyright (c) 1985, 1986 by the MathWorks, Inc.
if nargin < 2, k = 1; end
if nargin < 3, dn = 2; end
for i=1:k
[m,n] = size(x);
if m == 1
x = x(1+dn:n) - x(1:n-dn);
else
x = x(1+dn:m,:) - x(1:m-dn,:);
end
end