timour_matlab/mscripts/meanpsb1.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) ;


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

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

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

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

sum=zeros(Nx,Ny);
for k=Nz:-1:1,
sum=sum+dz*0.0002*(meant(:,:,k)-20);
end


cmax=max(max(sum));
cmin=min(min(sum));
cc=cmax-cmin;
c1=cmin+0.25*cc;
c2=cmin+0.75*cc;
V=[cmin cmax];
figure
caxis('manual')

        title='mean density anomaly';
        imagesc(lat,long,sum(142:230,18:106)');shading flat;axis image;caxis(V);colorbar('vertical');
        set(gca,'ydir','norm')
        text(50,130,descriptor);text(50,-30,num2str(nitend-nitstart+1));text(80,-30,'files');
        text(0,130,title);
        text(0,-10,'from timestep');text(30,-10,num2str(itstart));
        text(70,-10,'to'); text(80,-10,num2str(itend));

hold on

contour(sum',[c1 c1],'b')

hold on

contour(sum',[c2 c2],'b')

hold off

NNN=Ny/2;
ii=find(sum(:,NNN)>c1);
iw1=ii(1)-1+(c1-sum(ii(1)-1,NNN))/(sum(ii(1),NNN)-sum(ii(1)-1,NNN))
% iii=size(ii);
% ie1=ii(iii(1))+0.5

ii=find(sum(:,NNN)>c2);
iw2=ii(1)-1+(c2-sum(ii(1)-1,NNN))/(sum(ii(1),NNN)-sum(ii(1)-1,NNN))
% iii=size(ii);
% ie2=ii(iii(1))+0.5
xi=iw2-iw1

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
19			iterations
20
21			itstart = input(' Please enter start iteration : ','s')
22			itend = input(' Please enter end iteration : ','s')
23
24
25			sizeit=size(iterations);
26			for i=1:sizeit(1)
27			iter(i)=eval(iterations(i,1:10));
28			end
29			nitstart=find(iter==eval(itstart))
30			nitend=find(iter==eval(itend))
31
32			path = this_path
33			cmdstr=['cd ' path ];
34			eval(cmdstr);
35			path=pwd
36
37			sum=zeros(Nx,Ny,Nz);
38			counter=0;
39
40			for i=nitstart:nitend
41			tfilename=(['T.' iterations((i),1:10) ]) ;
42			t=rdmds(tfilename,'b');
43			sum=sum+t;
44			counter=counter+1;
45			end
46			meant=sum/counter;
47
48			i=find(meant<20);
49			meant(i)=20;
50
51			sum=zeros(Nx,Ny);
52			for k=Nz:-1:1,
53			sum=sum+dz0.0002(meant(:,:,k)-20);
54			end
55
56
57			cmax=max(max(sum));
58			cmin=min(min(sum));
59			cc=cmax-cmin;
60			c1=cmin+0.25*cc;
61			c2=cmin+0.75*cc;
62			V=[cmin cmax];
63			figure
64			caxis('manual')
65
66			title='mean density anomaly';
67			imagesc(lat,long,sum(142:230,18:106)');shading flat;axis image;caxis(V);colorbar('vertical');
68			set(gca,'ydir','norm')
69			text(50,130,descriptor);text(50,-30,num2str(nitend-nitstart+1));text(80,-30,'files');
70			text(0,130,title);
71			text(0,-10,'from timestep');text(30,-10,num2str(itstart));
72			text(70,-10,'to'); text(80,-10,num2str(itend));
73
74			hold on
75
76			contour(sum',[c1 c1],'b')
77
78			hold on
79
80			contour(sum',[c2 c2],'b')
81
82			hold off
83
84			NNN=Ny/2;
85			ii=find(sum(:,NNN)>c1);
86			iw1=ii(1)-1+(c1-sum(ii(1)-1,NNN))/(sum(ii(1),NNN)-sum(ii(1)-1,NNN))
87			% iii=size(ii);
88			% ie1=ii(iii(1))+0.5
89
90			ii=find(sum(:,NNN)>c2);
91			iw2=ii(1)-1+(c2-sum(ii(1)-1,NNN))/(sum(ii(1),NNN)-sum(ii(1)-1,NNN))
92			% iii=size(ii);
93			% ie2=ii(iii(1))+0.5
94			xi=iw2-iw1
95