function plotpv( iter)
clear 
clear path 

global Nx Ny Nz
global lat long dz dm mdep
global delt_su su_its t_su delt
global descriptor this_path
global f deltaf Q beta r_expt r_heat H
global time rots it
global g Cp rho_bar alpha
global u v t w
global iterations


param_file_name =   ...
   input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
feval(param_file_name) ;

iterations

it = input(' Please enter iteration : ','s') 

path   = this_path 
cmdstr=['cd ' path ];
eval(cmdstr);
path=pwd

ufilename=(['U.' it ]);
vfilename=(['V.' it ]);
tfilename=(['T.' it ]);
u=rdmds(ufilename,'b');
v=rdmds(vfilename,'b');
t=rdmds(tfilename,'b');

pv=zeros(Nx,Ny,Nz);
ff=zeros(Nx,Ny,Nz);
dx=dm;
dy=dm;

for j=1:Ny
ff(:,j,:)=f+(j-1)*dy*beta;
end

pv(1:Nx-1,1:Ny-1,1:Nz-1)=(v(2:Nx,1:Ny-1,1:Nz-1)/dx -v(1:Nx-1,1:Ny-1,1:Nz-1)/dx ...
                    -u(1:Nx-1,2:Ny,1:Nz-1)/dy +u(1:Nx-1,1:Ny-1,1:Nz-1)/dy ...
		    +ff(1:Nx-1,1:Ny-1,1:Nz-1)) ...
		    .*(t(1:Nx-1,1:Ny-1,1:Nz-1)-t(1:Nx-1,1:Ny-1,2:Nz))/dz;

t0=input('Enter temperature : ')

% Interpolate PV on T=T0 surface

h=zeros(Nz,Ny);
pvh=zeros(Nz,Ny);

for i=1:Nx
for j=1:Ny
kk=find(t(i,j,:)<t0);

if kk(1)>1
h(i,j)=(kk(1)-1)*dz+dz*(t(i,j,kk(1)-1)-t0)/(t(i,j,kk(1)-1)-t(i,j,kk(1)));
pvh(i,j)=pv(i,j,kk(1)-1) ...
+(pv(i,j,kk(1))-pv(i,j,kk(1)-1))*(t(i,j,kk(1)-1)-t0)/(t(i,j,kk(1)-1)-t(i,j,kk(1)));
else
h(i,j)=0;
pvh(i,j)=0;
end

end
end

hmax=max(max(h))


%cmin=0;
%cmax=0.1;
%V=[cmin,cmax];
%contourf(h,50);caxis(V);colorbar;grid

%        title='depth of a density interface';
%        imagesc(lat,long,h);shading flat;axis image;caxis(V);colorbar('horizontal');
%        text(0,-30,descriptor);
%        text(0,140,title);

figure
pcolor(pvh.');shading flat; colorbar; axis square
title(['pv timestep ' num2str(eval(it)) ' level of t=' num2str(t0) ]);
drawnow

figure

pcolor(h.');shading flat; colorbar; axis square
title(['h timestep' num2str(eval(it)) ' level of t=' num2str(t0) ]);
drawnow

return