% PV diagnostics 'DIAG4'

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=Ny-3;
Nz1=5;
Nz2=Nz-4;
pv1=pv;
Mp1=Mp;
Dp1=Dp;
z=zeros(Nz,1);
for k=1:Nz,
z(k)=(k-1)*dz;
end

for i=1:Nx,
for j=1:Ny;
Mp1(i,j,find(z > h(i,j)))=0.*Mp1(i,j,find(z > h(i,j)));
Dp1(i,j,find(z > h(i,j)))=0.*Dp1(i,j,find(z > h(i,j)));
end
end

ii=find((Nz1-1)*dz > h);

%  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
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*Mp1(3:Nx-2,3:Ny-2,k)*dz;
 sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+a*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);
 sumW(ii)=0.*sumW(ii);

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.01*(i);
v1(i)=0.01*(i);
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
x=Nx1:Nx2;
y=Ny1:Ny2;
subplot(2,1,1)
	contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v)
hold on
	contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
        text(10,Ny2+5,'Buoyancy eddy-transfer integral','Fontsize',13);
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 eddy-transfer 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(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v)

 set(gca,'DataAspectRatio',[1,1,1])
        text(10,Ny2+5,'Momentum eddy-transfer integral','Fontsize',13);
hold on
	contour(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
hold off
 xlabel('X (gridpoints)')
 ylabel('Y (gridpoints)')
        text(40,-20,'countour interval - 0.01 (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;
sumWW=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;
sumWW=sumWW+a*b*sumW(i,j)*dx*dy;
end
end
'DISSIPATION BY EDDY-DIFFUSIVITY'
sum
sumWW

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