timour_matlab/mscripts/asym1.m

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;


h=zeros(Nx,Ny);
hh=zeros(Nx,1);

t0=input('Enter temperature : ')
% t0=20.5;
meantheta(:,:,Nz)=0.;

for i=1:Nx,
for j=1:Ny,
kk=find(meantheta(i,j,:)<t0);
if kk(1)>1
h(i,j)=(kk(1)-1)*dz+dz*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
else
h(i,j)=0;
end
end
end

% hmax=max(max(h))

% [I,J]=find(h==hmax)
% NNN=J


%
NNN=Ny/2
%
hh=h(:,NNN);
hmax=max(hh);
I=find(hh==hmax);

ii=find(h(:,NNN)>0);
ii1=size(ii);
istart=ii(1)
iend=ii(ii1(1))
x0=(I-istart)/(iend-istart)


plot(hh);
1	edhill	1.1	clear path
2
3			global Nx Ny Nz
4			global lat long dz dm mdep
5			global delt_su su_its t_su delt
6			global descriptor this_path
7			global f deltaf Q beta r_expt r_heat H
8			global time rots it
9			global g Cp rho_bar alpha
10			global u v t w
11			global iterations
12
13
14			param_file_name = ...
15			input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
16			feval(param_file_name) ;
17
18			% compute the E-W asymmetry of the thermocline and its depth
19
20			iterations
21
22			itstart = input(' Please enter start iteration : ','s')
23			itend = input(' Please enter end iteration : ','s')
24
25
26			sizeit=size(iterations);
27			for i=1:sizeit(1)
28			iter(i)=eval(iterations(i,1:10));
29			end
30			nitstart=find(iter==eval(itstart))
31			nitend=find(iter==eval(itend))
32
33			path = this_path
34			cmdstr=['cd ' path ];
35			eval(cmdstr);
36			path=pwd
37
38			sumtheta=zeros(Nx,Ny,Nz);
39			counter=0;
40
41			for i=nitstart:nitend
42			tfilename=(['T.' iterations((i),1:10) ]) ;
43			t=rdmds(tfilename,'b');
44			sumtheta=sumtheta+t;
45			counter=counter+1;
46			end
47
48			meantheta=sumtheta/counter;
49
50
51
52
53			h=zeros(Nx,Ny);
54			hh=zeros(Nx,1);
55
56			t0=input('Enter temperature : ')
57			% t0=20.5;
58			meantheta(:,:,Nz)=0.;
59
60			for i=1:Nx,
61			for j=1:Ny,
62			kk=find(meantheta(i,j,:)<t0);
63			if kk(1)>1
64			h(i,j)=(kk(1)-1)dz+dz(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
65			else
66			h(i,j)=0;
67			end
68			end
69			end
70
71			% hmax=max(max(h))
72
73			% [I,J]=find(h==hmax)
74			% NNN=J
75
76
77			%
78			NNN=Ny/2
79			%
80			hh=h(:,NNN);
81			hmax=max(hh);
82			I=find(hh==hmax);
83
84			ii=find(h(:,NNN)>0);
85			ii1=size(ii);
86			istart=ii(1)
87			iend=ii(ii1(1))
88			x0=(I-istart)/(iend-istart)
89
90
91			plot(hh);