function [field,x,y]=peigui_depthslice(filein,depth,time,var)

%       USAGE: [field,x,y]=peigui_depthslice(filein,depth,time,var)
%
% peigui_DEPTHSLICE  Returns a matrix containing a horizontal slice (field),
%   of the variable var and the  position (x,y) at a 
%   specified depth and time step from an HOPS output (PE) file. 
%   The HOPS model is a sigtma model, what mean that to pass form the
%    the output file to depth levels is necesary an interpolation, this
%    have been optimized using the matlab vectors.
%
% 		filein -- the .nc file
% 		time -- time of teh simulatio
%		var -- string with the variable name
%		depth -- depth where the slice is taken
%
% Pedro Velez (UIB-IEO), July 1999


if strcmp('pbar',strrep(var,' ',''))
[field,x,y]=peigui_2dslice(filein,time,var);
field=field/1.e+12;
disp(sprintf('>>>>Pbar is a 2d variable \n'));
return;
end

if strcmp('qbar',strrep(var,' ',''))
[field,x,y]=peigui_2dslice(filein,time,var);
disp(sprintf('>>>>qbar is a 2d variable \n'));
return;
end

if strcmp('|Ut|',strrep(var,' ',''))
[u,x,y]=hopsgui_depthsliceu(filein,depth,time,'vtot');
[v,x,y]=hopsgui_depthslicev(filein,depth,time,'vtot');
field=sqrt(u.^2+v.^2);
return;
end

if strcmp('|Uc|',strrep(var,' ',''))
[u,x,y]=peigui_depthsliceu(filein,depth,time,'vclin');
[v,x,y]=peigui_depthslicev(filein,depth,time,'vclin');
field=sqrt(u.^2+v.^2);
return;
end

if strcmp('|Ug|',strrep(var,' ',''))
[u,x,y]=peigui_depthsliceu(filein,depth,time,'vgeo');
[v,x,y]=peigui_depthslicev(filein,depth,time,'vgeo');
field=sqrt(u.^2+v.^2);
return;
end

if strcmp('ut',strrep(var,' ',''))
[field,x,y]=peigui_depthsliceu(filein,depth,time,'vtot');
return;
end

if strcmp('vt',strrep(var,' ',''))
[field,x,y]=peigui_depthslicev(filein,depth,time,'vtot');
return;
end

if strcmp('ug',strrep(var,' ',''))
[field,x,y]=peigui_depthsliceu(filein,depth,time,'vgeo');
return;
end

if strcmp('vg',strrep(var,' ',''))
[field,x,y]=peigui_depthslicev(filein,depth,time,'vgeo');
return;
end


if strcmp('uc',strrep(var,' ',''))
[field,x,y]=peigui_depthsliceu(filein,depth,time,'vclin');
return;
end

if strcmp('vc',strrep(var,' ',''))
[field,x,y]=peigui_depthslicev(filein,depth,time,'vclin');
return;
end



if  strcmp('potden',strrep(var,' ',''))
[temp,x,y]=peigui_depthslice(filein,depth,time,'temp');
[salt,x,y]=peigui_depthslice(filein,depth,time,'salt');
field=potden80(salt,potemp80(salt,temp,ones(size(temp))*-depth),ones(size(temp))*-depth);
return;
end


vargrid='tgrid3';
if strcmp('wvsv',strrep(var,' ',''))
	vargrid='wvgrid3';
elseif strcmp('wvst',strrep(var,' ',''))
	vargrid='wtgrid3';
elseif strcmp('wvzv',strrep(var,' ',''))	
	vargrid='vgrid3';
end

x=getnc(filein,vargrid,[-1 -1 1 1],[-1 -1 1 1],1,1,2);
y=getnc(filein,vargrid,[-1 -1 1 2],[-1 -1 1 2],1,1,2);
zgrid3=getnc(filein,vargrid,[-1 -1 -1 3],[-1 -1 -1 3],1,1,2);
uz=getnc(filein,var,[time -1 -1 -1],[time -1 -1 -1],1,1,2);

%uz(find(uz==0))=NaN;
	
[ny,nx,nz]=size(zgrid3);
uz=reshape(uz,nx*ny,nz)';
zlev=reshape(zgrid3,nx*ny,nz)';

% find all the values of zlev that are less than the requested depth
% zind consists of zeros and ones.  One if the cell depth is less
% than (below) the requested z level, zero if cell level is greater.

zind=zlev < depth;

% find the indices IK of the cells that are just above the required level 

dz=diff(zind);
dz(nz,:)=zeros(size(dz(nz-1,:)));
ik=find(dz==1);

% find the indices ISUM where the requested depth is between two data values
isum=find(sum(dz));
field=zeros(nx*ny,1);

% Do linear interpolation
field(isum)=uz(ik)+(uz(ik+1)-uz(ik)).*(depth-zlev(ik))./(zlev(ik+1)-zlev(ik));

% put back in matrix form
field=reshape(field,ny,nx);

field(find(field==0))=NaN;

land=getnc(filein,'landt',[-1 -1],[-1 -1],1,1,1);
field(find(land==0))=NaN;




