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	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