getdpth.m

%======================================================================
%                    G E T D P T H . M 
%                    doc: Fri Nov 19 12:17:32 2021
%                    dlm: Fri Nov 19 12:24:52 2021
%                    (c) 2021 A.M. Thurnherr
%                    uE-Info: 19 0 NIL 0 0 72 0 2 4 NIL ofnI
%======================================================================

function [d,p]=getdpth(d,p)
% function [d,p]=getdpth(d,p)
% LADCP-2 processing software v 5.0
%
%  -- make depth from raw data 
% Matin Visbeck and Gerd Krahmann April-2000
% changed December 2002

% CHANGES BY ANT:
%	Nov 19, 2021: - added error message

disp('GETDEPTH: Integrating depth from vertical velocity')

% set default start, deepst  and end depth
p=setdefv(p,'zpar',[0 NaN 0]);
p=setdefv(p,'cut',15);
p=setdefv(p,'dzbelow',[2 -1]*medianan(abs(diff(d.zd))));
p.zpar=abs(p.zpar);

% get time difference for w-values
dt=diff(d.time_jul)*24*3600;
dt=mean([0,dt;dt,0]);

d.izmflag = d.rw*0;

p=setdefv(p,'ctddepth',0);

figure(4)
orient tall


% two sweeps to properly remove bottom values

d=getmeanw(d,p);

for n=1:2

  % set depth below which to delete all data 
  dzbelow=p.dzbelow(n);

  disp([' starting run ',int2str(n),' to find bottom depth'])

  % integrate vertical velocity to obtain depth
  %  assume zero depth at beginning and end of cast
  %  first sweep to get sound speed profile


  if n==1
   zzd=cumsum(d.wm.*dt);
   [zzdmax,zzdimax]=max(zzd);
   disp([' maxdepth form down-trace is ',num2str(zzdmax)])
   zzu=fliplr(cumsum(fliplr(-d.wm.*dt)));
   [zzumax,zzuimax]=max(zzu);
   disp([' maxdepth from up-trace   is ',num2str(zzumax)])
   % find mean maximum depth
   zzimax=round((zzdimax+zzuimax)/2);
   zzmax=(zzd(zzimax)+zzu(zzimax))/2;
   % merge up and down trace depth by forcing them to agree
   zz=[zzd(1:zzimax)/zzd(zzimax)*zzmax,zzu(zzimax+1:end)/zzu(zzimax)*zzmax];
   zz0=zz;
   if existf(d,'ctdprof_ss')
    disp(' take soundspeed from CTD profile ')
    zctd=d.ctdprof_z;
    zctd(1)=-1e5;
    zctd(end)=1e5;
    ss=interp1q(zctd,d.ctdprof_ss,zz')';
   else
    disp(' make soundspeed based on pressure and ADCP temp')
    pp=press(abs(zz));
    ss=sounds(pp,d.temp(1,:),34.5);
   end
   % sound speed correction
   if d.soundc==0
    sc=ss./d.sv(1,:);
    d.wm=d.wm.*sc;
    if existf(d,'hsurf')
     d.hsurf=d.hsurf.*sc;
    end
    if existf(d,'hbot')
     d.hbot=d.hbot.*sc;
    end
   end
  end

  % to surface corrections
  if ~isfinite(p.zpar(1)), p.zpar(1)=10; end
  if ~isfinite(p.zpar(3)), p.zpar(3)=10; end
  zsc=zz(1)-p.zpar(1);
  zec=zz(end)-p.zpar(3);

  [dum,ibot]=max(zz);
  % for very shallow stations turn of surface detection
  if dum<100, p.surfdist=0; disp(' shallow station no surface detection '),end
  ii=1:ibot;
  iok=ii(find(zz(ii)<200 & zz(ii)>30 ));
  if length(iok)>2
   iok2=1:iok(end);
   subplot(321)
   plot(iok2,-zz(iok2),'.r',iok2,iok2*0,'-k'),
   hold on
   title('start depth (.red) (.blue) Surface distance')
   xlabel('time in ensembles')
   ylabel('depth in meter')
   ax=axis;
   ax(4)=max(ax(4),20);
   axis(ax);
  end

  ii=ibot:length(zz);
  iok=ii(find(zz(ii)<200 & zz(ii)>30 ));
  if length(iok)>2
   iok2=iok(1):length(zz);
   subplot(322)
   plot(iok2,-zz(iok2),'.r',iok2,iok2*0,'-k'),
   hold on
   title('end depth')
   xlabel('time in ensembles')
   ylabel('depth in meter')
   ax=axis;
   ax(4)=max(ax(4),20);
   axis(ax);
  end

  % surface distance to find start depth
  if existf(d,'hsurf') & p.surfdist
   if sum(isfinite(d.hsurf))>10
    [zmax,ibot]=max(zz);
   % start depth
    ii=find(isfinite(d.hsurf(1:ibot)));
    iok=ii(find(zz(ii)<200 & zz(ii)>40 & abs(d.hsurf(ii)-zz(ii))<30));
    if ~isempty(iok)
     % use surface return for start depth
      start_depth_corr=median(d.hsurf(iok)-zz(iok));
    % temporary plot of surface detection
      iok2=1:iok(end);
      subplot(321)
      if n==1
       plot(iok,d.hsurf(iok)-zz(iok),'.',iok2,-zz(iok2),'.r'),
       hold on, plot(iok2,iok2*0+start_depth_corr,'--b') 
      else
       plot(iok,d.hsurf(iok)-zz(iok),'.',iok2,-zz(iok2),'.r'),
      end
      if std(d.hsurf(iok)-zz(iok))<20
       zsc=-start_depth_corr;
       disp([' start depth correction from surface return: ',int2str(zsc)])
      end
      ax=axis;
      ax(4)=max(ax(4),20);
      axis(ax);
    end

   % surface distance to find end depth
    ii=find(isfinite(d.hsurf(ibot:end)))+ibot-1;
    iok=ii(find(zz(ii)<200 & zz(ii)>40 & abs(d.hsurf(ii)-zz(ii))<30));
    if ~isempty(iok)
     % use surface return for end depth
      end_depth_corr=median(d.hsurf(iok)-zz(iok));
    % temporary plot of surface detection
      subplot(322)
      iok2=iok(1):length(zz);
      if n==1
       plot(iok,d.hsurf(iok)-zz(iok),'.',iok2,-zz(iok2),'.r'),
       hold on, plot(iok2,iok2*0+end_depth_corr,'--b') 
      else
       plot(iok,d.hsurf(iok)-zz(iok),'.',iok2,-zz(iok2),'.r'),
      end
      if std(d.hsurf(iok)-zz(iok))<20
       zec=-end_depth_corr;
       disp([' end depth correction from surface return: ',int2str(zec)])
      end
      ax=axis;
      ax(4)=max(ax(4),20);
      axis(ax);
    end
   end 
  end 

  disp([' last depth from int W is :',num2str(zz0(end))])
  disp([' should be                :',num2str(p.zpar(3))])
  % correct for start and end depth
  d.z_uncorr = -zz0;
  zz1=linspace(zsc,zec,length(zz));
  if length(zz1)>1
    zz=zz-zz1;
  else
    zz=zz-(zz(1)-p.zpar(1));
  end

  disp([' maximum depth from int W is :',num2str(max(zz))])
  disp([' should be                   :',num2str(p.zpar(2))])
  % correct for maximum depth
  if isfinite(p.zpar(2))
    zz=zz/max(zz)*p.zpar(2);
  end

  if p.ctddepth==0 
   % save results only if CTD-depth was not availabel
   d.z=-zz;
   p.ladcpdepth=1;
   disp(' use LADCP depth from integrated W')
  else
   p.ladcpdepth=0;
  end

  if p.ctddepth==1 & n>1
   d.z_ladcp=-zz;
   zz=-d.z;
   dz=d.z_ladcp-d.z;
   ii=find(isfinite(dz));
   p.ladcpr_CTD_depth_std=[mean(dz), std(dz)];
   disp(' use CTD time series depth ')
   disp([' LADCP minus CTD depth mean: ',num2str(p.ladcpr_CTD_depth_std(1)),...
          '  std: ',num2str(p.ladcpr_CTD_depth_std(2))]);
   if ~isfinite(p.ladcpr_CTD_depth_std(1))
	error('non-numeric result, try reprocessing with p.getdepth = 2');
   end
  end
  [p.maxdepth,ibottom]=max(-d.z);

  % bottom depth
  p.zbottom=NaN;
  subplot(312)
  iok=find((max(zz)-zz)<200);
  iok2=iok(1):iok(end);
  plot(iok2,d.z(iok2),'.r'),
  if sum(isfinite(d.hbot))>10
    % look for bottom only close to deepest CTD depth
    iok=find((max(zz)-zz)<200 & d.hbot>0 & abs(d.wm)>0);
    if ~isempty(iok)
      % fit polynom to bottom depth time series
      c=polyfit(iok,d.hbot(iok)-d.z(iok),1);
      if n>1
      % use deepest point to set bottom depth
       zbottomerr= polyval(c,iok)-(d.hbot(iok)-d.z(iok)) ;
       iok=iok(find(abs(zbottomerr)<1.5*std(zbottomerr) | abs(zbottomerr)<50 ));
       c=polyfit(iok,d.hbot(iok)-d.z(iok),2);
       zbottomerr= polyval(c,iok)-(d.hbot(iok)-d.z(iok)) ;
       p.zbottom=polyval(c,ibottom);
      else
       p.zbottom=medianan(d.hbot(iok)-d.z(iok));
       zbottomerr= p.zbottom-(d.hbot(iok)-d.z(iok)) ;
      end
      p.zbottomerror = medianan(abs(zbottomerr)); 
    % temporary plot of bottom detection
     iok2=iok(1):iok(end);
      subplot(312)
      plot(iok,-d.hbot(iok)+d.z(iok),'.',iok2,d.z(iok2),'.r'),
      hold on, plot(iok2,iok2*0-p.zbottom,'--k') 
      hold on, plot(iok2,-polyval(c,iok2),'-b')
    % remove outlier
      ii=find(abs(zbottomerr)>2*std(zbottomerr) | abs(zbottomerr)>100 );
      if n==1
       hold on, plot(iok(ii),-d.hbot(iok(ii))+d.z(iok(ii)),'*g'), 
      end
      title('bottom (--k)  bottom distance (.b)')
      xlabel('time in ensembles')
      ylabel('depth in meter')

      d.hbot(iok(ii))=NaN;
      d.bvel(iok(ii),:)=NaN;
    else
      p.zbottomerror = nan;
    end
   % check if bottom is shallower that maxctd-depth an
    if ((p.zbottom-p.maxdepth<-20 & isfinite(p.zbottom)) |...
       	p.zbottomerror > 20 )
      disp('  no bottom found')
      disp(['   given maximum profile depth : ',int2str(p.maxdepth)])
      disp(['   extracted bottom depth      : ',int2str(p.zbottom)])
      disp(['        bottom depth error     : ',int2str(p.zbottomerror)])
      p.zbottom=NaN;
    else
      disp(['  bottom found at ',int2str(p.zbottom),' +/- ',...
                                 int2str(p.zbottomerror),' m'])
      if (p.zbottom<p.maxdepth)
        disp('  extracted bottom within 10m above given maximum profile depth')
      end
    end
  end
  pause(0.1)

  [izm1,izm]=meshgrid([fliplr(d.zu),-d.zd],d.z);
  d.izm=izm'+izm1';

  % flag all data below bottom as bad
  if ~isnan(p.zbottom)
    ii = find(d.izm<-p.zbottom-dzbelow);
    d.izmflag(ii)=NaN; 
  end

  if length(d.zu)>0
   % flag all data close to the surface as bad
   ii = find(d.izm>-d.zu(1));
   d.izmflag(ii)=NaN; 
  end
  

end 


% set velocities deeper than bottom to NaN
bad = find( isnan(d.izmflag) & isfinite(d.ru) );
if ~isempty(bad)
  disp([' removing ',int2str(length(bad)),...
	' values  below recognized bottom'])
end

%d.ru = d.ru+d.izmflag;
%d.rv = d.rv+d.izmflag;
%d.rw = d.rw+d.izmflag;
d.weight = d.weight+d.izmflag;

% compute pressure from depth
d.p=press(abs(d.z));

if d.z(1)<-50, 
 warn=[' first LADCP depth is ',int2str(d.z(1))];
 disp(warn)
 p.warn(size(p.warn,1)+1,1:length(warn))=warn;
end

if d.z(end)<-50, 
 warn=[' last LADCP depth is ',int2str(d.z(end))];
 disp(warn)
 p.warn(size(p.warn,1)+1,1:length(warn))=warn;
end

% get sound speed time series if not available already
if ~existf(d,'ss')
 if existf(d,'ctd_ss')
    disp(' take soundspeed from CTD time series')
  d.ss=d.ctd_ss;
 elseif existf(d,'ctdprof_ss')
    disp(' take soundspeed from CTD profile ')
    zctd=d.ctdprof_z;
    zctd(1)=-1e5;
    zctd(end)=1e5;
    d.ss=interp1q(zctd,d.ctdprof_ss,d.z')';
 else
   if existf(d,'ctd_temp')
    disp(' make soundspeed based on pressure and CTD temp')
    d.ss=sounds(d.p,d.ctd_temp,34.5);
   else
    disp(' make soundspeed based on pressure and ADCP temp')
    d.ss=sounds(d.p,d.temp(1,:),34.5);
   end
 end
end

% correct velocity for sound speed
if d.soundc==0
    disp(' correct velocities for sound speed ')
    sc=meshgrid(d.ss./d.sv(1,:),d.izd);
    d.ru(d.izd,:)=d.ru(d.izd,:).*sc;
    d.rv(d.izd,:)=d.rv(d.izd,:).*sc;
    d.rw(d.izd,:)=d.rw(d.izd,:).*sc;
    if length(d.zd)~=length(d.ru(:,1))
     sc=meshgrid(d.ss./d.sv(2,:),d.izu);
     d.ru(d.izu,:)=d.ru(d.izu,:).*sc;
     d.rv(d.izu,:)=d.rv(d.izu,:).*sc;
     d.rw(d.izu,:)=d.rw(d.izu,:).*sc;
    end
    d.soundc=1;
else
    disp(' will not correct for sound speed twice')
end

% cut raw data to only include profile
i1=find(d.z(1:ibottom)>-p.cut);
i2=find(d.z(ibottom:end)>-p.cut)+ibottom-1;
if length(i1)==0
 i1=1;
end
if length(i2)==0
 i2=length(d.z);
end
ii=d.z*0;
ic=i1(end):i2(1);
ii(ic)=1;
if (sum(ii)~=length(ii)) & p.cut>0 & existf(d,'cutindx')~=1
 disp(' remove data at begining and end of cast')
 d=cutstruct(d,ii);
 subplot(321), ax=axis; plot([1 1]*ic(1),ax(3:4),'--k')
 subplot(322), ax=axis; plot([1 1]*ic(end),ax(3:4),'--k')
 %p.time_start=gregoria(d.time_jul(1));
 %p.time_end=gregoria(d.time_jul(end));
 p.zpar(1)=max([0,-d.z(1)]);
 p.zpar(3)=max([0,-d.z(end)]);
end

streamer([p.name,'  Figure 4']);
pause(0.01)

%-----------------------------------------------------------------
function d=getmeanw(d,p)
% function d=getmeanw(d,p)
  [d.wm1,dum]=medianan(d.rw(p.wizr,:)+d.izmflag(p.wizr,:),1);
  ii=find(sum(isfinite(dum))<2);
  d.wm1(ii)=NaN;
  d.wm=d.wm1;
  ii=2:(length(d.wm1)-1);
  d.wm(ii)=meannan([d.wm1(ii-1);d.wm1(ii);d.wm1(ii);d.wm1(ii+1)]);
  % try to replace bad W-data by propagating them from the ends
  for dummy=1:3
   ii=find(~isfinite(d.wm));
   if length(ii)==0, break, end
   dat=[];
   for nn=1:dummy
    i1=ii-nn;
    i=find(i1<1);
    i1(i)=1;
    i2=ii+nn;
    i=find(i2>length(d.wm));
    i2(i)=length(d.wm);
    dat=[dat;d.wm(i1);d.wm(i2)];
   end
   d.wm(ii)=meannan(dat);
  end
  % replace missing w-profiles by zero
  ii=find(~isfinite(d.wm));
  d.wm(ii)=0;


%==============================================================
function a=cutstruct(a,ii)
% reduce array size in structure
lz=length(ii);
iok=find(ii==1);
a.cutindx=[iok(1) iok(end)];
if isstruct(a)
  fnames = fieldnames(a);
  for n=1:size(fnames,1)
   dummy = getfield(a,fnames{n});
   [ly,lx]=size(dummy);
   if ly==lz
    a=setfield(a,fnames{n},dummy(iok,:));
   elseif lx==lz
    a=setfield(a,fnames{n},dummy(:,iok));
   end
  end
end