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

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))]);

  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