t0=input('Enter temperature : ')

% Interpolate Dw on T=T0 surface

h=zeros(Nx,Ny);
hh=zeros(Nx,Ny);
hhh=zeros(Nx,Ny);
Dh=zeros(Nx,Ny);

Nz1=4;
he=dz*(Nz1-1);

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

%figure

%pcolor(hhh.');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) ]);
%figure


we=0;
ws=0;
WE=zeros(Nx,Ny);
WS=zeros(Nx,Ny);
hh=zeros(Nx,Ny);

%for i=Nx/2:Nx
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;
WE(i,j)=meanw(i,j,Nz1);
WS(i,j)=Dh(i,j);

% if r<0.40
if h(i,j)>he
we=meanw(i,j,Nz1)*dx*dx+we;
ws=Dh(i,j)*dx*dx+ws;
% end
end
end
end

hh=h;

for k=1:3
WE1(:,:,k)=WE(:,:);
WS1(:,:,k)=WS(:,:);
hh1(:,:,k)=hh(:,:);
end


for k=1:5
WE1=smooth3(WE1);
WS1=smooth3(WS1);
hh1=smooth3(hh1);
end

hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)+hh1(3:Nx,2:Ny-1,1)+...
hh1(2:Nx-1,1:Ny-2,1)+hh1(2:Nx-1,3:Ny,1)-4*hh1(2:Nx-1,2:Ny-1,1));

%hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)-hh1(3:Nx,2:Ny-1,1)).*...
%(hh1(1:Nx-2,2:Ny-1,1)-hh1(3:Nx,2:Ny-1,1))+...
%10000*(hh1(2:Nx-1,1:Ny-2,1)-hh1(2:Nx-1,3:Ny,1)).*...
%(hh1(2:Nx-1,1:Ny-2,1)-hh1(2:Nx-1,3:Ny,1));

i=find(h<0.015);
h(i)=0;
% i=find(abs(hhh)>0.001);
% hhh(i)=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
%hhh(i,j)=NaN;
%hh(i,j)=NaN;
%WS1(i,j,1)=NaN;
%end

if h(i,j)<he
WE1(i,j,1)=0;
WS1(i,j,1)=NaN;
hh(i,j)=NaN;
hhh(i,j)=NaN;
end

end
end


%pcolor(hh(1:Nx,:).');shading flat; colorbar; axis square
% set(gca,'DataAspectRatio',[2,2,2])
% title('h','FontSize',16);
%figure

x=1:3:Nx-2;
y=1:3:Ny;
hh(:,:)=max(he,hh(:,:));
subplot('position',[0.30 0.7 0.4 0.3])
 surf(x,y,200*(-hh(1:3:Nx-2,1:3:Ny)'+he));
axis([1 Nx-4 0 Ny-2 -10.0 0])
a=0:1:63;
a=a/63;
cmap=[a' a' a'];
colormap(cmap);
set(gca,'DataAspectRatio',[20,20,7])
set(gca,'XtickLabel','|')
set(gca,'YtickLabel','|')
set(gca,'ZtickLabel','|')
xlabel('x','FontSize',20)
ylabel('y','FontSize',20)
zlabel('z','FontSize',20,'Rotation',[0])
%title('h','FontSize',19);
view(-20,10);

hhh(1:Nx,:)=-10*hhh(1:Nx,:);
% hhh(Nx/2:Nx,:)=100*hhh(Nx/2:Nx,:);
a=max(max(hhh(1:Nx,:)));

%pcolor(hhh(1:Nx,:).');caxis([0 a]);shading flat; colorbar; axis square
%set(gca,'DataAspectRatio',[2,2,2])
%title('-\nabla^2 h','FontSize',12);
%xlabel('x','Fontsize',13);
%ylabel('y','Fontsize',13);
%text(50,-10,'x');
%text(-10,60,'y');

%pcolor(squeeze(WE1(Nx/2:Nx,:,1)).');shading flat; colorbar; axis square
% set(gca,'DataAspectRatio',[2,2,2])
%title('Ek, smooth');
%figure

%pcolor(meanw(Nx/2:Nx,:,3).');shading flat; colorbar; axis square
%set(gca,'DataAspectRatio',[2,2,2])
%title('Ek');
figure

subplot('position',[0.15 0.05 0.7 0.5])
WS1(1:Nx,:,1)=-WS1(1:Nx,:,1)*16667;
a=max(max(WS1(1:Nx,:,1)));
pcolor(squeeze(WS1(1:Nx,:,1)).');caxis([0 a]);shading flat;colormap('default'); 
set(gca,'DataAspectRatio',[2,2,2])
set(gca,'Xtick',[1 100],'XtickLabel','|')
set(gca,'Ytick',[1 120],'YtickLabel','|')
title('w^*','FontSize',20);
xlabel('x','Fontsize',20);
ylabel('y','Fontsize',20,'Rotation',[0]);


we
ws