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

mufilename  =  (['mudgeu.' it ]);
mvfilename  =  (['mudgev.' it ]);
ufilename   =  (['U.' it ]);
vfilename   =  (['V.' it ]);
tfilename   =  (['T.' it ]);
udfilename  =  (['Udsk.' it ]);
vdfilename  =  (['Vdsk.' it ]);
usfilename  =  (['Uspin.' it ]);
vsfilename  =  (['Vspin.' it ]);
mudgeu      =  rdmeta(mufilename,'b');
mudgev      =  rdmeta(mvfilename,'b');
udisk       =  rdmeta(udfilename,'b');
vdisk       =  rdmeta(vdfilename,'b');
uspin       =  rdmeta(usfilename,'b');
vspin       =  rdmeta(vsfilename,'b');
u           =  rdmeta(ufilename,'b');
v           =  rdmeta(vfilename,'b');
t           =  rdmeta(tfilename,'b');

w=zeros(Nx,Ny,Nz+1);
dx=dm;
dy=dm;

for k=Nz:-1:1,
 w(1:Nx-1,1:Ny-1,k)= w(1:Nx-1,1:Ny-1,k+1)...
                    -u(2:Nx,1:Ny-1,k)*dz(k)/dx +u(1:Nx-1,1:Ny-1,k)*dz(k)/dx ...
                    -v(1:Nx-1,2:Ny,k)*dz(k)/dy +v(1:Nx-1,1:Ny-1,k)*dz(k)/dy ;
end
                udminusu=-udisk-u(:,:,1);
                udminusus=-udisk-uspin(:,:,1);
                vdminusv=vdisk-v(:,:,1);
                vdminusvs=vdisk-vspin(:,:,1);
figure
subplot(4,4,1);pcolor(udisk);axis square;colorbar;shading flat;
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('udisk');
subplot(4,4,2);pcolor(vdisk);axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('vdisk') 
subplot(4,4,3);pcolor(sqrt(vdisk.*vdisk+udisk.*udisk));axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('diskspeed')
subplot(4,4,4);pcolor(vspin);axis square;colorbar;shading flat                          
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('vspin')    
subplot(4,4,5);pcolor(uspin);axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('uspin')                                 
subplot(4,4,6);pcolor(sqrt(vdisk.*vdisk+udisk.*udisk));axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('spinspeed') 
subplot(4,4,7);pcolor(mudgeu);axis square;colorbar;shading flat                         
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('mudgeu')   
subplot(4,4,8);pcolor(mudgev);axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('mudgev')                                 
subplot(4,4,9);pcolor(u(:,:,1));axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('u')     
subplot(4,4,10);pcolor(v(:,:,1));axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('v')                                        
subplot(4,4,11);pcolor(sqrt(v(:,:,1).*v(:,:,1)+u(:,:,1).*u(:,:,1)));axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('speed')
subplot(4,4,12);pcolor(udminusu);axis square;colorbar;shading flat                       
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('udminu')
subplot(4,4,13);pcolor(udminusus);axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('udminusus')
subplot(4,4,14);pcolor(vdminusv);axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('vdminusv') 
subplot(4,4,15);pcolor(vdminusvs);axis square;colorbar;shading flat
set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]);title('vdminusvs')
subplot(4,4,16);set(gca,'XTickLabel',[]);set(gca,'YTickLabel',[]); title(descriptor)