t1=input('Enter temperature 1: ')
%t2=input('Enter temperature 2: ')

%t1=25;
t2=25;


% Interpolate Dw on T=T0 surface

h1=zeros(Nx,Ny);
h2=zeros(Nx,Ny);
hh=zeros(Nx,Ny);
hhh=zeros(Nx,Ny);
Dh1=zeros(Nx,Ny);
Dh2=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,:)<t1);

if kk(1)>1
h1(i,j)=(kk(1)-1)*dz+dz*(meantheta(i,j,kk(1)-1)-t1)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));

Dh1(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)-t1)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
else
h1(i,j)=0;
Dh1(i,j)=0;
end

end
end

for i=1:Nx
for j=1:Ny
kk=find(meantheta(i,j,:)<t2);

if kk(1)>1
h2(i,j)=(kk(1)-1)*dz+dz*(meantheta(i,j,kk(1)-1)-t2)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));

Dh2(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)-t2)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
else
h2(i,j)=0;
Dh2(i,j)=0;
end

end
end

hmax=max(max(h2))


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=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 h1(i,j)>he
we=meanw(i,j,Nz1)*dx*dx+we;
ws=Dh(i,j)*dx*dx+ws;
% end
end
end
end

hh=h1;

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 h1(i,j)<(he+dz)
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(:,:));
surf(x,y,200*(-hh(1:3:Nx-2,1:3:Ny)'+he));
axis([0 Nx-2 0 Ny -10.0 0])
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)
%title('h','FontSize',19);
view(-20,15);
figure

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

subplot(2,1,1)
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);


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

subplot(2,1,2)
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; colorbar; axis square;
set(gca,'DataAspectRatio',[2,2,2])
 title('w_1^*','FontSize',12);
xlabel('x','Fontsize',13);
ylabel('y','Fontsize',13);


we
ws