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