timour_matlab/mscripts/DIAGW.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

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

for i=nitstart:nitend
i,0
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');

ufilename=(['U.' iterations((i),1:10) ]) ;
u=rdmds(ufilename,'b');
u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));

vfilename=(['V.' iterations((i),1:10) ]) ;
v=rdmds(vfilename,'b');
v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));

wfilename=(['W.' iterations((i),1:10) ]) ;
w=rdmds(wfilename,'b');

sumtheta=sumtheta+t;
sumu=sumu+u;
sumv=sumv+v;
sumw=sumw+w;
counter=counter+1;
end

meantheta=sumtheta/counter;
meanu=sumu/counter;
meanv=sumv/counter;
meanw=sumw/counter;

sumut=zeros(Nx,Ny,Nz);
sumvt=zeros(Nx,Ny,Nz);
sumwt=zeros(Nx,Ny,Nz);

meanu=smooth3(meanu);
meanv=smooth3(meanv);
meanw=smooth3(meanw);
meantheta=smooth3(meantheta);

for i=nitstart:nitend
i
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');

ufilename=(['U.' iterations((i),1:10) ]) ;
u=rdmds(ufilename,'b');
u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));

vfilename=(['V.' iterations((i),1:10) ]) ;
v=rdmds(vfilename,'b');
v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));

wfilename=(['W.' iterations((i),1:10) ]) ;
w=rdmds(wfilename,'b');

t=t-meantheta;
u=u-meanu;
v=v-meanv;
w=w-meanw;

sumut=sumut+u.*t;
sumvt=sumvt+v.*t;
sumwt=sumwt+w.*t;

end

sumut=sumut/counter;
sumvt=sumvt/counter;
sumwt=sumwt/counter;


D=zeros(Nx,Ny,Nz);

dx=dm;
dy=dm;
dz=-dz;

D(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumut(3:Nx,2:Ny-1,2:Nz-1)-sumut(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
        -(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
        -(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);

Dw=zeros(Nx,Ny,Nz);
ff=zeros(Nx,Ny,Nz);


for j=1:Ny
ff(:,j,:)=f+(j-1)*dy*beta;
end

Dw(:,:,2:Nz-1)=D(:,:,2:Nz-1)*2*dz./(meantheta(:,:,3:Nz)-meantheta(:,:,1:Nz-2));

i=find(meantheta(:,:,3:Nz)==meantheta(:,:,1:Nz-2));
Dw(i)=0;


t0=input('Enter temperature : ')

% Interpolate Dw on T=T0 surface

h=zeros(Nx,Ny);
Dh=zeros(Nx,Ny);
dz=-dz

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

Dh(i,j)=Dw(i,j,kk(1)-1) ...
+(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ...
*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
else
h(i,j)=0;
Dh(i,j)=0;
end

end
end

hmax=max(max(h))

figure

pcolor(Dh.');shading flat; colorbar; axis square
 title('mean Wstar');
% title(['mean Wstar from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]);

figure

pcolor(h.');shading flat; colorbar; axis square
%pcolor(flipud(h.'));shading flat; colorbar; axis square
%title(['mean H from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]);


we=0;
ws=0;
for i=1:Nx
for j=1:Ny

r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx;

%if r<0.44

we=meanw(i,j,3)*dx*dx+we;
ws=Dh(i,j)*dx*dx+ws;
%end
end
end

we
ws


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			% iterations
19
20			itstart = input(' Please enter start iteration : ','s')
21			itend = input(' Please enter end iteration : ','s')
22
23
24			sizeit=size(iterations);
25			for i=1:sizeit(1)
26			iter(i)=eval(iterations(i,1:10));
27			end
28			nitstart=find(iter==eval(itstart))
29			nitend=find(iter==eval(itend))
30
31			path = this_path
32			cmdstr=['cd ' path ];
33			eval(cmdstr);
34			path=pwd
35
36			sumtheta=zeros(Nx,Ny,Nz);
37			sumu=zeros(Nx,Ny,Nz);
38			sumv=zeros(Nx,Ny,Nz);
39			sumw=zeros(Nx,Ny,Nz);
40			counter=0;
41
42			for i=nitstart:nitend
43			i,0
44			tfilename=(['T.' iterations((i),1:10) ]) ;
45			t=rdmds(tfilename,'b');
46
47			ufilename=(['U.' iterations((i),1:10) ]) ;
48			u=rdmds(ufilename,'b');
49			u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
50
51			vfilename=(['V.' iterations((i),1:10) ]) ;
52			v=rdmds(vfilename,'b');
53			v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
54
55			wfilename=(['W.' iterations((i),1:10) ]) ;
56			w=rdmds(wfilename,'b');
57
58			sumtheta=sumtheta+t;
59			sumu=sumu+u;
60			sumv=sumv+v;
61			sumw=sumw+w;
62			counter=counter+1;
63			end
64
65			meantheta=sumtheta/counter;
66			meanu=sumu/counter;
67			meanv=sumv/counter;
68			meanw=sumw/counter;
69
70			sumut=zeros(Nx,Ny,Nz);
71			sumvt=zeros(Nx,Ny,Nz);
72			sumwt=zeros(Nx,Ny,Nz);
73
74			meanu=smooth3(meanu);
75			meanv=smooth3(meanv);
76			meanw=smooth3(meanw);
77			meantheta=smooth3(meantheta);
78
79			for i=nitstart:nitend
80			i
81			tfilename=(['T.' iterations((i),1:10) ]) ;
82			t=rdmds(tfilename,'b');
83
84			ufilename=(['U.' iterations((i),1:10) ]) ;
85			u=rdmds(ufilename,'b');
86			u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
87
88			vfilename=(['V.' iterations((i),1:10) ]) ;
89			v=rdmds(vfilename,'b');
90			v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
91
92			wfilename=(['W.' iterations((i),1:10) ]) ;
93			w=rdmds(wfilename,'b');
94
95			t=t-meantheta;
96			u=u-meanu;
97			v=v-meanv;
98			w=w-meanw;
99
100			sumut=sumut+u.*t;
101			sumvt=sumvt+v.*t;
102			sumwt=sumwt+w.*t;
103
104			end
105
106			sumut=sumut/counter;
107			sumvt=sumvt/counter;
108			sumwt=sumwt/counter;
109
110
111			D=zeros(Nx,Ny,Nz);
112
113			dx=dm;
114			dy=dm;
115			dz=-dz;
116
117			D(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumut(3:Nx,2:Ny-1,2:Nz-1)-sumut(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
118			-(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
119			-(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
120
121			Dw=zeros(Nx,Ny,Nz);
122			ff=zeros(Nx,Ny,Nz);
123
124
125			for j=1:Ny
126			ff(:,j,:)=f+(j-1)dybeta;
127			end
128
129			Dw(:,:,2:Nz-1)=D(:,:,2:Nz-1)2dz./(meantheta(:,:,3:Nz)-meantheta(:,:,1:Nz-2));
130
131			i=find(meantheta(:,:,3:Nz)==meantheta(:,:,1:Nz-2));
132			Dw(i)=0;
133
134
135			t0=input('Enter temperature : ')
136
137			% Interpolate Dw on T=T0 surface
138
139			h=zeros(Nx,Ny);
140			Dh=zeros(Nx,Ny);
141			dz=-dz
142
143			meantheta(:,:,Nz)=0.;
144			for i=1:Nx
145			for j=1:Ny
146			kk=find(meantheta(i,j,:)<t0);
147
148			if kk(1)>1
149			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)));
150
151			Dh(i,j)=Dw(i,j,kk(1)-1) ...
152			+(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ...
153			*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
154			else
155			h(i,j)=0;
156			Dh(i,j)=0;
157			end
158
159			end
160			end
161
162			hmax=max(max(h))
163
164			figure
165
166			pcolor(Dh.');shading flat; colorbar; axis square
167			title('mean Wstar');
168			% title(['mean Wstar from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]);
169
170			figure
171
172			pcolor(h.');shading flat; colorbar; axis square
173			%pcolor(flipud(h.'));shading flat; colorbar; axis square
174			%title(['mean H from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]);
175
176
177			we=0;
178			ws=0;
179			for i=1:Nx
180			for j=1:Ny
181
182			r=sqrt((i-3Nx/4)(i-3Nx/4)+(j-Ny/2)(j-Ny/2))*dx;
183
184			%if r<0.44
185
186			we=meanw(i,j,3)dxdx+we;
187			ws=Dh(i,j)dxdx+ws;
188			%end
189			end
190			end
191
192			we
193			ws
194
195
196
197
198
199