% PV diagnostics 'DIAGtest'

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;
% for single gyre
% Ny2=Ny-3;
% for double-gyre
Ny2=45;

Nz1=5;
Nz2=Nz-4;

%Nx1=30;
%Nx2=32;
%Ny1=30;
%Ny2=32;
%Nz1=10;
%Nz2=12;

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

%  for kk=1:4
%  pv1=smooth3(pv1);
%  Mp1=smooth3(Mp1);
%  Dp1=smooth3(Dp1);
%  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);

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

for k=1:15

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:Ny-2;
subplot(2,1,1)
	contourf(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny-2,1))',10);colorbar;
%	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','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';
%        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)
	contourf(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny-2,1))',10);colorbar;
%	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-------------------------------------------------

dz

bx=zeros(Nx,Ny,Nz);
by=zeros(Nx,Ny,Nz);
bz=zeros(Nx,Ny,Nz);
zetax=zeros(Nx,Ny,Nz);
zetay=zeros(Nx,Ny,Nz);
zetaz=zeros(Nx,Ny,Nz);

bx(2:Nx-1,2:Ny-1,2:Nz-1)=(meantheta(3:Nx,2:Ny-1,2:Nz-1)-meantheta(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx);
by(2:Nx-1,2:Ny-1,2:Nz-1)=(meantheta(2:Nx-1,3:Ny,2:Nz-1)-meantheta(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy);
bz(2:Nx-1,2:Ny-1,2:Nz-1)=-(meantheta(2:Nx-1,2:Ny-1,3:Nz)-meantheta(2:Nx-1,2:Ny-1,1:Nz-2))/(2*abs(dz));

zetax(2:Nx-1,2:Ny-1,2:Nz-1)=(meanv(2:Nx-1,2:Ny-1,3:Nz)-meanv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*abs(dz));
zetay(2:Nx-1,2:Ny-1,2:Nz-1)=-(meanu(2:Nx-1,2:Ny-1,3:Nz)-meanu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*abs(dz));
zetaz(2:Nx-1,2:Ny-1,2:Nz-1)=(meanv(3:Nx,2:Ny-1,2:Nz-1)/(2*dx)-meanv(1:Nx-2,2:Ny-1,2:Nz-1)/(2*dx) ...
                    -meanu(2:Nx-1,3:Ny,2:Nz-1)/(2*dy)+meanu(2:Nx-1,1:Ny-2,2:Nz-1)/(2*dy) ...
                    +ff(2:Nx-1,2:Ny-1,2:Nz-1));


sum=0;
sumT=0;
sumV=0;
sumY=0;
summ=0;
sumd=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


summ=summ+a*b*(Mu(i,j,Nz1)*by(i,j,Nz1)-Mv(i,j,Nz1)*bx(i,j,Nz1))*dx*dy;
summ=summ-a*b*(Mu(i,j,Nz2)*by(i,j,Nz2)-Mv(i,j,Nz2)*bx(i,j,Nz2))*dx*dy;
sumd=sumd+a*b*zetaz(i,j,Nz1)*D(i,j,Nz1)*dx*dy;
sumd=sumd-a*b*zetaz(i,j,Nz2)*D(i,j,Nz2)*dx*dy;
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;
sumV=sumV+a*b*meanw(i,j,Nz1)*dx*dy;
sumV=sumV-a*b*meanw(i,j,Nz2)*dx*dy;
end
end

'xy'
summ
sumV
sumd

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



summ=summ+a*b*Mv(Nx2,j,k)*bz(Nx2,j,k)*dz*dy;
summ=summ-a*b*Mv(Nx1,j,k)*bz(Nx1,j,k)*dz*dy;
sumd=sumd+a*b*zetax(Nx2,j,k)*D(Nx2,j,k)*dz*dy;
sumd=sumd-a*b*zetax(Nx1,j,k)*D(Nx1,j,k)*dz*dy;
 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;
 sumV=sumV+a*b*meanu(Nx2,j,k)*dz*dy;
 sumV=sumV-a*b*meanu(Nx1,j,k)*dz*dy;
end
end

'xy+yz'
summ
sumd

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


summ=summ-a*b*Mu(i,Ny2,k)*bz(i,Ny2,k)*dz*dx;
summ=summ+a*b*Mu(i,Ny1,k)*bz(i,Ny1,k)*dz*dx;
sumd=sumd+a*b*zetay(i,Ny2,k)*D(i,Ny2,k)*dz*dx;
sumd=sumd-a*b*zetay(i,Ny1,k)*D(i,Ny1,k)*dz*dx;
sumT=sumT+a*b*meanv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
sumY=sumY+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;
sumV=sumV+a*b*meanv(i,Ny2,k)*dz*dx;
sumV=sumV-a*b*meanv(i,Ny1,k)*dz*dx;
end
end
'TOTAL'
sum
sumY
sumT
sumV
summ
sumd


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

'IMBALANCE'
(sumT-(summ+sumd))/sumd