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
sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+Mp(3:Nx-2,3:Ny-2,k).*pv(3:Nx-2,3:Ny-2,k)*dz;
sumN(3:Nx-2,3:Ny-2)=sumN(3:Nx-2,3:Ny-2)+Np(3:Nx-2,3:Ny-2,k).*pv(3:Nx-2,3:Ny-2,k)*dz;
sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+Dp(3:Nx-2,3:Ny-2,k).*pv(3:Nx-2,3:Ny-2,k)*dz;
end
sumW(3:Nx-2,3:Ny-2)=0.5*meanw(3:Nx-2,3:Ny-2,Nz1).*pv(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:20

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

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


%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(Nx1:Nx2,Ny1:Ny2,1))',v)
hold on
        contour(squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
        text(10,98,'Buoyancy diffusion integral','Fontsize',14);
hold off
% xlabel('X (gridpoints)')
% ylabel('Y (gridpoints)')
 set(gca,'DataAspectRatio',[1,1.6,2])

subplot(3,1,2)
        contour(squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v)

 set(gca,'DataAspectRatio',[1,1.6,1])
        text(10,98,'Momentum diffusion integral','Fontsize',14);
hold on
        contour(squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
hold off
% xlabel('X (gridpoints)')
% ylabel('Y (gridpoints)')


subplot(3,1,3)
        contour(squeeze(sumN1(Nx1:Nx2,Ny1:Ny2,1))',v)

 set(gca,'DataAspectRatio',[1,1.6,1])
        text(10,98,'Nonlinear terms','Fontsize',14);
hold on
        contour(squeeze(sumN1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
hold off
% xlabel('X (gridpoints)')
% ylabel('Y (gridpoints)')

        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(Nx1:Nx2,Ny1:Ny2,1))',2*v)
hold on
        contour(squeeze(sumN1(Nx1:Nx2,Ny1:Ny2,1))',2*v1,'--')
hold off




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

sum=0;
sumT=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
sum=sum+0.5*meanw(i,j,Nz1)*pv(i,j,Nz1)*pv(i,j,Nz1)*dx*dy;
sumT=sum-0.5*meanw(i,j,Nz2)*pv(i,j,Nz2)*pv(i,j,Nz2)*dx*dy;
end
end

for j=Ny1:Ny2
for k=Nz1:Nz2
 sumT=sumT+0.5*meanu(Nx2,j,k)*pv(Nx2,j,k)*pv(Nx2,j,k)*dz*dy;
 sumT=sumT-0.5*meanu(Nx1,j,k)*pv(Nx1,j,k)*pv(Nx1,j,k)*dz*dy;
end
end

for i=Nx1:Nx2
for k=Nz1:Nz2
sumT=sumT+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
sumT=sumT-0.5*meanv(i,Ny1,k)*pv(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
sum=sum+sumD(i,j)*dx*dy;
end
end
'DISSIPATION BY EDDY-DIFFUSIVITY'
sum

sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
sum=sum+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