timour_matlab/mscripts/getdonut.m

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
1	edhill	1.1	global Nx Ny Nz
2			global dx dy lat long
3			global dz depth
4			global delt iteration it
5			global descriptor
6			global this_path
7			global g Cp rho_bar alpha
8			global f0 B0 Q0 beta r_heat r_expt H
9			global u v w t
10
11			param_file_name = 'datadonut';
12			feval(param_file_name) ;
13
14			iteration
15
16			it = input(' Please enter iteration : ')
17
18			path = this_path
19			eval(['cd ' path '/it',int2str(it)]);
20			path=pwd
21
22			fprintf(1,' Reading output ... \n') ;
23
24			tfn='theta.sun.b';
25			ufn='u.sun.b';
26			vfn='v.sun.b';
27
28			% T (m/s)
29			t=Readmodel(tfn,[Nx Ny Nz],'float32');
30			t=Extract([3 Nx Ny Nz t(5:4+NxNyNz)],[1 Nx 1 Ny 1 Nz]);
31			t=reshape(t,Nx,Ny,Nz);
32
33			% U (m/s)
34			u=Readmodel(ufn,[Nx Ny Nz],'float32');
35			u=Extract([3 Nx Ny Nz u(5:4+NxNyNz)],[1 Nx 1 Ny 1 Nz]);
36			u=reshape(u,Nx,Ny,Nz);
37
38			% V (m/s)
39			v=Readmodel(vfn,[Nx Ny Nz],'float32');
40			v=Extract([3 Nx Ny Nz v(5:4+NxNyNz)],[1 Nx 1 Ny 1 Nz]);
41			v=reshape(v,Nx,Ny,Nz);
42
43			% W (m/s)
44
45			w=zeros(Nx,Ny,Nz+1);
46
47			for k=Nz:-1:1,
48			w(1:Nx-1,1:Ny-1,k)= w(1:Nx-1,1:Ny-1,k+1)...
49			-u(2:Nx,1:Ny-1,k)dz(k)/dx +u(1:Nx-1,1:Ny-1,k)dz(k)/dx ...
50			-v(1:Nx-1,2:Ny,k)dz(k)/dy +v(1:Nx-1,1:Ny-1,k)dz(k)/dy ;
51			end
52
53			return