| 1 | 
 | 
| 2 | 
Dw=zeros(Nx,Ny,Nz); | 
| 3 | 
 | 
| 4 | 
dz=-abs(dz); | 
| 5 | 
 | 
| 6 | 
Dw(:,:,2:Nz-1)=D(:,:,2:Nz-1)*2*dz./(meantheta(:,:,3:Nz)-meantheta(:,:,1:Nz-2)); | 
| 7 | 
i=find(meantheta(:,:,3:Nz)==meantheta(:,:,1:Nz-2)); | 
| 8 | 
Dw(i)=0; | 
| 9 | 
 | 
| 10 | 
%t0=input('Enter temperature : ') | 
| 11 | 
t0=0.5*(max(max(max(meantheta)))+min(min(min(meantheta(1:Nx-1,1:Ny-1,1:Nz-1))))); | 
| 12 | 
 | 
| 13 | 
% Interpolate Dw on T=T0 surface | 
| 14 | 
 | 
| 15 | 
h=zeros(Nx,Ny); | 
| 16 | 
Dh=zeros(Nx,Ny); | 
| 17 | 
dz=abs(dz); | 
| 18 | 
 | 
| 19 | 
meantheta(:,:,Nz)=0.; | 
| 20 | 
for i=1:Nx | 
| 21 | 
for j=1:Ny | 
| 22 | 
kk=find(meantheta(i,j,:)<t0); | 
| 23 | 
 | 
| 24 | 
if kk(1)>1 | 
| 25 | 
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))); | 
| 26 | 
Dh(i,j)=Dw(i,j,kk(1)-1) ... | 
| 27 | 
+(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ... | 
| 28 | 
*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1))); | 
| 29 | 
else | 
| 30 | 
h(i,j)=0; | 
| 31 | 
Dh(i,j)=0; | 
| 32 | 
end | 
| 33 | 
 | 
| 34 | 
end | 
| 35 | 
end | 
| 36 | 
 | 
| 37 | 
hmax=max(max(h)) | 
| 38 | 
 | 
| 39 | 
figure | 
| 40 | 
 | 
| 41 | 
pcolor(Dh.');shading flat; colorbar; axis square | 
| 42 | 
% title('mean Wstar'); | 
| 43 | 
% title(['mean Wstar from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]); | 
| 44 | 
 | 
| 45 | 
figure | 
| 46 | 
 | 
| 47 | 
pcolor(h.');shading flat; colorbar; axis square | 
| 48 | 
 |