global Nx Ny Nz
global dx dy lat long 
global dz depth
global delt iteration it
global descriptor 
global this_path
global g Cp rho_bar alpha
global f0 B0 Q0 beta r_heat r_expt H 
global u v w t  

param_file_name =  'datadonut';
feval(param_file_name) ;

iteration

it = input(' Please enter iteration : ') 

path   = this_path 
eval(['cd ' path '/it',int2str(it)]);
path=pwd

fprintf(1,' Reading output ... \n') ;

tfn='theta.sun.b';
ufn='u.sun.b';
vfn='v.sun.b';

% T (m/s)
t=Readmodel(tfn,[Nx Ny Nz],'float32');
t=Extract([3 Nx Ny Nz t(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
t=reshape(t,Nx,Ny,Nz);

% U (m/s)
u=Readmodel(ufn,[Nx Ny Nz],'float32');
u=Extract([3 Nx Ny Nz u(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
u=reshape(u,Nx,Ny,Nz);

% V (m/s)
v=Readmodel(vfn,[Nx Ny Nz],'float32');
v=Extract([3 Nx Ny Nz v(5:4+Nx*Ny*Nz)],[1 Nx 1 Ny 1 Nz]);
v=reshape(v,Nx,Ny,Nz);

% W (m/s)

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

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

return