sumM=zeros(Nx,Ny);
sumN=zeros(Nx,Ny);
sumD=zeros(Nx,Ny);
sumW=zeros(Nx,Ny);


dz=0.005;

% DEFINE BOX OF INTEGRATION
Nx1=4;
Nx2=Nx-3;
Ny1=4;
Ny2=Ny-3;
Nz1=5;
Nz2=Nz-4;


for k=Nz1:Nz2
if k==Nz1
a=0.5;
elseif k==Nz2
a=0.5;
else
a=1;
end

sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+a*Mp(3:Nx-2,3:Ny-2,k)*dz;
sumN(3:Nx-2,3:Ny-2)=sumN(3:Nx-2,3:Ny-2)+a*Np(3:Nx-2,3:Ny-2,k)*dz;
sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+a*Dp(3:Nx-2,3:Ny-2,k)*dz;
end
sumW(3:Nx-2,3:Ny-2)=meanw(3:Nx-2,3:Ny-2,Nz1).*pv(3:Nx-2,3:Ny-2,Nz1);


for k=1:3
sumM1(:,:,k)=sumM(:,:);
sumN1(:,:,k)=sumN(:,:);
sumD1(:,:,k)=sumD(:,:);
sumW1(:,:,k)=sumW(:,:);
end

for k=1:25

sumM1=smooth3(sumM1);
sumD1=smooth3(sumD1);
sumW1=smooth3(sumW1);
end

for k=1:10
sumN1=smooth3(sumN1);
end

v=zeros(25,1);
v1=zeros(25,1);
for i=1:25
% v(i)=-0.1*(i-0.5)*7;
% v1(i)=0.1*(i-0.5)*7;
v(i)=0.001*(i)*10;
v1(i)=-0.001*(i)*10;
end

sumD1(4:Nx,Ny,1)=0;
sumM1(4:Nx,Ny,1)=0;
sumN1(4:Nx,Ny,1)=0;
sumD1(4:Nx,Ny-1,1)=0;
sumM1(4:Nx,Ny-1,1)=0;
sumN1(4:Nx,Ny-1,1)=0;



%V=[-10 10];
%	title='Vorticity input';
%        imagesc(lat,long,sumW1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
%        set(gca,'ydir','norm')
%        text(0,110,title);
%figure
%	title='Buoyancy diffusion integral';
%        imagesc(lat,long,sumD1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
%        set(gca,'ydir','norm')
%        text(0,110,title);
%figure
%	title='Non-linear terms integral';
%        imagesc(lat,long,sumN1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
%        set(gca,'ydir','norm')
%        text(0,110,title);
%figure
%	title='Momentum diffusion integral';
%        imagesc(lat,long,sumM1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
%        set(gca,'ydir','norm')
%        text(0,110,title);
subplot(3,1,1)
        contour(squeeze(sumD1(4:Nx,1:Ny,1))',v,'k')
hold on
        contour(squeeze(sumD1(4:Nx,1:Ny,1))',v1,'k--')
        text(10,107,'Buoyancy transfer','Fontsize',14);
        text(-20,90,'a)','Fontsize',18);
hold off
 xlabel('X ')
 ylabel('Y ')
 set(gca,'XtickLabel','||')
 set(gca,'YtickLabel','||')

 set(gca,'DataAspectRatio',[1,1.6,2])

subplot(3,1,2)
        contour(squeeze(sumM1(4:Nx,1:Ny,1))',v,'k')

 set(gca,'DataAspectRatio',[1,1.6,1])
        text(10,107,'Momentum transfer','Fontsize',14);
        text(-20,90,'b)','Fontsize',18);
hold on
        contour(squeeze(sumM1(4:Nx,1:Ny,1))',v1,'k--')
hold off
 xlabel('X ')
 ylabel('Y ')
 set(gca,'XtickLabel','||')
 set(gca,'YtickLabel','||')



subplot(3,1,3)
        contour(squeeze(sumN1(4:Nx,1:Ny,1))',v,'k')

 set(gca,'DataAspectRatio',[1,1.6,1])
        text(10,107,'Inertial term ','Fontsize',14);
        text(-20,90,'c)','Fontsize',18);
hold on
        contour(squeeze(sumN1(4:Nx,1:Ny,1))',v1,'k--')
hold off
 xlabel('X ')
 ylabel('Y ')
 set(gca,'XtickLabel','||')
 set(gca,'YtickLabel','||')


%        text(40,-25,'countour interval - 1.5 (non-dimensional units)','FontSize',7);
%        text(40,-30,'negative values - solid line','FontSize',7);
%        text(40,-35,'positive values - dashed line','FontSize',7);
figure
        contour(squeeze(sumN1(1:Nx,1:Ny,1))',v)
hold on
        contour(squeeze(sumN1(1:Nx,1:Ny,1))',v1,'--')
hold off




%-----TEST-------------------------------------------------

sum=0;
sumT=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
if i==Nx1
a=0.5;
elseif i==Nx2
a=0.5;
else
a=1;
end

if j==Ny1
b=0.5;
elseif j==Ny2
b=0.5;
else
b=1;
end

sum=sum+a*b*meanw(i,j,Nz1)*pv(i,j,Nz1)*dx*dy;
sumT=sum-a*b*meanw(i,j,Nz2)*pv(i,j,Nz2)*dx*dy;
end
end

for j=Ny1:Ny2
for k=Nz1:Nz2
if k==Nz1
a=0.5;
elseif k==Nz2
a=0.5;
else
a=1;
end

if j==Ny1
b=0.5;
elseif j==Ny2
b=0.5;
else
b=1;
end

 sumT=sumT+a*b*meanu(Nx2,j,k)*pv(Nx2,j,k)*dz*dy;
 sumT=sumT-a*b*meanu(Nx1,j,k)*pv(Nx1,j,k)*dz*dy;
end
end

for i=Nx1:Nx2
for k=Nz1:Nz2
if k==Nz1
a=0.5;
elseif k==Nz2
a=0.5;
else
a=1;
end

if i==Nx1
b=0.5;
elseif i==Nx2
b=0.5;
else
b=1;
end

sumT=sumT+a*b*meanv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
sumT=sumT-a*b*meanv(i,Ny1,k)*pv(i,Ny1,k)*dz*dx;
end
end
'VORTICITY GENERATION'
sum
sumT


sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
if i==Nx1
a=0.5;
elseif i==Nx2
a=0.5;
else
a=1;
end

if j==Ny1
b=0.5;
elseif j==Ny2
b=0.5;
else
b=1;
end

sum=sum+a*b*sumD(i,j)*dx*dy;
end
end
'DISSIPATION BY EDDY-DIFFUSIVITY'
sum

sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
if i==Nx1
a=0.5;
elseif i==Nx2
a=0.5;
else
a=1;
end

if j==Ny1
b=0.5;
elseif j==Ny2
b=0.5;
else
b=1;
end

sum=sum+a*b*sumM(i,j)*dx*dy;
end
end
'DISSIPATION BY EDDY-VISCOSITY'
sum

sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
sum=sum+sumN(i,j)*dx*dy;
end
end
'NONLINEARITY'
sum