Dw=zeros(Nx,Ny,Nz);

dz=-abs(dz);

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 : ')
t0=0.5*(max(max(max(meantheta)))+min(min(min(meantheta(1:Nx-1,1:Ny-1,1:Nz-1)))));

% Interpolate Dw on T=T0 surface

h=zeros(Nx,Ny);
Dh=zeros(Nx,Ny);
dz=abs(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