% PV diagnostics 'DIAG5'
sumM=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=42;
Nz1=5;
Nz2=Nz-4;

pv1=pv;
Mp1=Mp;
Dp1=Dp;

%  pv1=smooth3(pv);
%  Mp1=smooth3(Mp);
%  Dp1=smooth3(Dp);

%  for kk=1:4
%  pv=smooth3(pv);
%  Mp=smooth3(Mp);
%  Dp=smooth3(Dp);
%  end

 for k=Nz1:Nz2
 sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+Mp1(3:Nx-2,3:Ny-2,k)*dz;
 sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+Dp1(3:Nx-2,3:Ny-2,k)*dz;
 end
 sumW(3:Nx-2,3:Ny-2)=meanw(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1);

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

for k=1:25

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

V=[0.25*min(min(sumW(Nx1+2:Nx2,Ny1:Ny2-4))) 0];

%	title='Vorticity input';
%        imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
%        set(gca,'ydir','norm')
%        text(0,110,title);
%figure
v=zeros(10,1);
v1=zeros(10,1);
for i=1:10
%v(i)=-0.1*(i-0.5)*7;
%v1(i)=0.1*(i-0.5)*7;
v(i)=-0.001*(i)*7;
v1(i)=0.001*(i)*7;
end

% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v)
% hold on
% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
% hold off 
%        text(20,0,'countour interval is 0.7 in the non-dimensional units');
%        text(20,95,'Vorticity input','Fontsize',16);

%figure

%V=[-0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))];
%	title='Buoyancy diffusion integral';
%        imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
%        set(gca,'ydir','norm')
%        text(0,110,title);


% figure

subplot(2,1,1)
	contour(squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',50)
        text(10,Ny2,'Buoyancy diffusion integral','Fontsize',15);
hold off
 xlabel('X (gridpoints)')
 ylabel('Y (gridpoints)')
 set(gca,'DataAspectRatio',[2,2,2])

%figure
%V=[-0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))];
%	title='Momentum diffusion integral';
%        imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
%        set(gca,'ydir','norm')
%        text(0,110,title);

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

 set(gca,'DataAspectRatio',[1,1,1])
        text(10,Ny2,'Momentum diffusion integral','Fontsize',15);
 xlabel('X (gridpoints)')
 ylabel('Y (gridpoints)')
        text(40,-25,'countour interval - 0.7 (non-dimensional units)','FontSize',7);
        text(40,-30,'negative values - solid line','FontSize',7);
        text(40,-35,'positive values - dashed line','FontSize',7);



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

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

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

for i=Nx1:Nx2
for k=Nz1:Nz2
sumT=sumT+meanv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
sumY=sumY+meanv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
sumT=sumT-meanv(i,Ny1,k)*pv(i,Ny1,k)*dz*dx;
end
end
'VORTICITY GENERATION'
sum
sumY
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