t0=input('Enter temperature : ')
% t0=27;

% Interpolate Dw on T=T0 surface

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

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

Nz1=4;
Nx1=4;
Nx2=Nx-3;
Ny1=4;
Ny2=Ny-3;

he=dz*(Nz1-1);
hmax=max(max(h))


close all


we=0;
ws=0;
WE=zeros(Nx,Ny);
WS=zeros(Nx,Ny);
hh=zeros(Nx,Ny);
WE1=zeros(Nx,Ny,3);
WS1=zeros(Nx,Ny,3);
hh1=zeros(Nx,Ny,3);

for i=Nx1:Nx2
for j=Ny1:Ny2

WE(i,j)=meanw(i,j,Nz1);
WS(i,j)=Dh(i,j);

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



end
end

we
ws

for i=Nx1:Nx2
for k=Nz1:Nz

z=dz*(k-1);
if z>he 
if z<h(i,Ny1)
we=we-meanv(i,Ny1,k)*dx*dz;
end
end

if z>he
if z<h(i,Ny2)
we=we+meanv(i,Ny2,k)*dx*dz;
end
end

end
end

for j=Ny1:Ny2
for k=Nz1:Nz

z=dz*(k-1);
if z>he
if z<h(Nx1,j)
we=we-meanu(Nx1,j,k)*dy*dz;
end
end

if z>he
if z<h(Nx2,j)
we=we+meanu(Nx2,j,k)*dy*dz;
end
end

end
end

we

'ratio', ws/we


hh=h;

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


for k=1:8
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));


i=find(h<he);
h(i)=0;
% i=find(abs(hhh)>0.001);
% hhh(i)=0;

for i=1:Nx
for j=1:Ny

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


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(:,:).');caxis([0 a]);shading flat; axis square
%set(gca,'DataAspectRatio',[2,2,2],'XtickLabel','|','YtickLabel','|')
%title('-\nabla^2 h','FontSize',15);
%xlabel('X','FontSize',15);
%ylabel('Y','Rotation',[0],'FontSize',13);

figure
WS1(1:Nx,:,1)=-WS1(1:Nx,:,1)*16667;
a=max(max(WS1(1:Nx,:,1)));
pcolor(squeeze(WS1(:,:,1)).');caxis([-0.3*a a]);shading flat; axis square;colorbar;
set(gca,'DataAspectRatio',[2,2,2],'XtickLabel','|','YtickLabel','|')
title('w^*','FontSize',15);
xlabel('X','FontSize',15);
ylabel('Y','Rotation',[0],'FontSize',13);


%WS2=squeeze(WS1(:,:,1));
%i=find(~isnan(hhh));
%WS2=WS2(i);
%WS=-WS*16667;
%hhh=hhh(i);
%sqrt(mean(WS2.*WS2)/ mean(hhh.*hhh))
%mean(WS2)/mean(hhh)
%max(WS2)/max(hhh)