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='data11';
feval(param_file_name) ;

itstart='0000200000';
itend='0000240000';
iii=5;
istart=142;
iend=230;
t0=25;


iterations

sizeit=size(iterations);
for i=1:sizeit(1)
iter(i)=eval(iterations(i,1:10));
end

nitstart=find(iter==eval(itstart));
nitend=find(iter==eval(itend));

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

sumtheta=zeros(Nx,Ny,Nz);
counter=0;

for i=nitstart:nitend
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');
sumtheta=sumtheta+t;
counter=counter+1;
end
meantheta=sumtheta/counter;

W1=smooth3(meantheta);

for i=1:iii,
W1=smooth3(W1);
end

hh=zeros(Nx,1);
for i=1:Nx,
kk=find(W1(i,Ny/2,:)<t0);
if kk(1)>1
hh(i)=(kk(1)-1)+(W1(i,Ny/2,kk(1)-1)-t0)/(W1(i,Ny/2,kk(1)-1)-W1(i,Ny/2,kk(1)));
else
hh(i)=0;
end
end
hmax=max(hh);
I=find(hh==hmax);
x0=(I-istart)/(iend-istart)


        contour(flipud(squeeze(W1(:,Ny/2,3:30))'),[25,25],'b')
	hold on

itstart='0000240000';
itend='0000300000';
nitstart=find(iter==eval(itstart));
nitend=find(iter==eval(itend));

sumtheta=zeros(Nx,Ny,Nz);
counter=0;
for i=nitstart:nitend
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');
sumtheta=sumtheta+t;
counter=counter+1;
end
meantheta=sumtheta/counter;
W2=smooth3(meantheta);
for i=1:iii,
W2=smooth3(W2);
end
hh=zeros(Nx,1);
for i=1:Nx,
kk=find(W2(i,Ny/2,:)<t0);
if kk(1)>1
hh(i)=(kk(1)-1)+(W2(i,Ny/2,kk(1)-1)-t0)/(W2(i,Ny/2,kk(1)-1)-W2(i,Ny/2,kk(1)));
else
hh(i)=0;
end
end
hmax=max(hh);
I=find(hh==hmax);
x0=(I-istart)/(iend-istart)
        contour(flipud(squeeze(W2(:,Ny/2,3:30))'),[25,25],'b')
	hold on

itstart='0000300000';
itend='0000350000';
nitstart=find(iter==eval(itstart));
nitend=find(iter==eval(itend));

sumtheta=zeros(Nx,Ny,Nz);
counter=0;
for i=nitstart:nitend
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');
sumtheta=sumtheta+t;
counter=counter+1;
end
meantheta=sumtheta/counter;
W3=smooth3(meantheta);
for i=1:iii,
W3=smooth3(W3);
end
hh=zeros(Nx,1);
for i=1:Nx,
kk=find(W3(i,Ny/2,:)<t0);
if kk(1)>1
hh(i)=(kk(1)-1)+(W3(i,Ny/2,kk(1)-1)-t0)/(W3(i,Ny/2,kk(1)-1)-W3(i,Ny/2,kk(1)));
else
hh(i)=0;
end
end
hmax=max(hh);
I=find(hh==hmax);
x0=(I-istart)/(iend-istart)
        contour(flipud(squeeze(W3(:,Ny/2,3:30))'),[25,25],'b')
        hold on

itstart='0000350000';
itend='0000400000';
nitstart=find(iter==eval(itstart));
nitend=find(iter==eval(itend));

sumtheta=zeros(Nx,Ny,Nz);
counter=0;
for i=nitstart:nitend
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');
sumtheta=sumtheta+t;
counter=counter+1;
end
meantheta=sumtheta/counter;
W4=smooth3(meantheta);
for i=1:iii,
W4=smooth3(W4);
end
hh=zeros(Nx,1);
for i=1:Nx,
kk=find(W4(i,Ny/2,:)<t0);
if kk(1)>1
hh(i)=(kk(1)-1)+(W4(i,Ny/2,kk(1)-1)-t0)/(W4(i,Ny/2,kk(1)-1)-W4(i,Ny/2,kk(1)));
else
hh(i)=0;
end
end
hmax=max(hh);
I=find(hh==hmax);
x0=(I-istart)/(iend-istart)
	contour(flipud(squeeze(W4(:,Ny/2,3:30))'),[25,25],'b')
        hold on

itstart='0000400000';
itend='0000500000';
nitstart=find(iter==eval(itstart));
nitend=find(iter==eval(itend));

sumtheta=zeros(Nx,Ny,Nz);
counter=0;
for i=nitstart:nitend
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');
sumtheta=sumtheta+t;
counter=counter+1;
end
meantheta=sumtheta/counter;
W5=smooth3(meantheta);
for i=1:iii,
W5=smooth3(W5);
end
hh=zeros(Nx,1);
for i=1:Nx,
kk=find(W5(i,Ny/2,:)<t0);
if kk(1)>1
hh(i)=(kk(1)-1)+(W5(i,Ny/2,kk(1)-1)-t0)/(W5(i,Ny/2,kk(1)-1)-W5(i,Ny/2,kk(1)));
else
hh(i)=0;
end
end
hmax=max(hh);
I=find(hh==hmax);
x0=(I-istart)/(iend-istart)
	contour(flipud(squeeze(W5(:,Ny/2,3:30))'),[25,25],'b')
        hold off