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

% compute the E-W asymmetry of the thermocline and its depth

iterations

itstart = input(' Please enter start iteration : ','s') 
itend = input(' Please enter end iteration : ','s') 


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;
meantheta(:,:,1)=20;
meantheta(:,:,2)=20;
meantheta(:,:,3)=20;


%cmin=20;
%cmax=24,
%V=[cmin cmax];
%figure
%caxis('manual')
%
%subplot(1,2,1)
%        title='average temperature at z=2';
%        imagesc(lat,long,squeeze(meantheta(:,:,2))');shading flat;axis image;caxis(V);colorbar('horizontal');
%        text(0,-30,descriptor);text(50,-30,num2str(nitend-nitstart+1));text(80,-30,'files');
%        text(0,140,title);
%        text(0,-20,'from timestep');text(50,-20,num2str(itstart));
%	text(0,-10,'to'); text(50,-10,num2str(itend));
%
%subplot(1,2,2)
%        contour(flipud(squeeze(meantheta(:,Ny/2,:))'),50);caxis(V);colorbar;grid
%        text(0,-2,'average temperature at mid-tank')
%


%
NNN=Ny/2
%

tx=zeros(Nx,1);

for i=1:Nx,

for k=1:Nz,
for j=NNN-10:NNN+10,
tx(i)=tx(i)+meantheta(i,j,k);
end
end
tx(i)=tx(i)/Nz;
tx(i)=tx(i)/21;
end

i=find(tx<20);
tx(i)=20;

tmax=max(tx);

lmax=find(tx==tmax);


tmin=min(tx);

te=tmin+(tmax-tmin)/2.71;

tsign=zeros(Nx-1,1);

for i=1:Nx-1
tsign(i)=(tx(i)-te)*(tx(i+1)-te);
end

ii=find(tsign<0);

xw=ii(1)+0.5;
xe=ii(2)+0.5;

x0=(lmax-xw)/(xe-xw)

txe=zeros(Nx,1);
txe(:)=te;

plot(tx)
hold on
plot(txe,'--')
hold off

% Depth
%
%tz=zeros(Nz,1);
%for k=1:Nz,
%for i=1:Nx,
%tz(k)=tz(k)+meantheta(i,NNN,k);
%end
%tz(k)=tz(k)/Nx;
%end
%
%tzmax=max(tz);
%tzmin=min(tz);
%
%kk=zeros(Nz-5,1);
%tza=zeros(Nz-5,1);
%for k=1:Nz-5
%kk(k)=k;
%tza(k)=tz(k)-tzmin;
%end
%
%logp1=polyfit(kk,log(tza),1);
%logpred1=2.71.^polyval(logp1,kk);
%figure
%semilogy(kk,logpred1,'-',kk,tza,'+');
%grid on
%
%h0=-1/logp1(1)