timour_matlab/mscripts/DIAG4.m

% 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
1	edhill	1.1	% PV diagnostics 'DIAG4'
2
3			sumM=zeros(Nx,Ny);
4			sumD=zeros(Nx,Ny);
5			sumW=zeros(Nx,Ny);
6
7			dz=0.005;
8
9			% DEFINE BOX OF INTEGRATION
10			Nx1=4;
11			Nx2=Nx-3;
12			Ny1=4;
13			Ny2=Ny-3;
14			Nz1=5;
15			Nz2=Nz-4;
16			pv1=pv;
17			Mp1=Mp;
18			Dp1=Dp;
19			z=zeros(Nz,1);
20			for k=1:Nz,
21			z(k)=(k-1)*dz;
22			end
23
24			for i=1:Nx,
25			for j=1:Ny;
26			Mp1(i,j,find(z > h(i,j)))=0.*Mp1(i,j,find(z > h(i,j)));
27			Dp1(i,j,find(z > h(i,j)))=0.*Dp1(i,j,find(z > h(i,j)));
28			end
29			end
30
31			ii=find((Nz1-1)*dz > h);
32
33			% pv1=smooth3(pv);
34			% Mp1=smooth3(Mp);
35			% Dp1=smooth3(Dp);
36
37			% for kk=1:4
38			% pv=smooth3(pv);
39			% Mp=smooth3(Mp);
40			% Dp=smooth3(Dp);
41			% end
42
43			for k=Nz1:Nz2
44			if k==Nz1
45			a=0.5;
46			elseif k==Nz2
47			a=0.5;
48			else
49			a=1;
50			end
51
52			sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+aMp1(3:Nx-2,3:Ny-2,k)dz;
53			sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+aDp1(3:Nx-2,3:Ny-2,k)dz;
54			end
55			sumW(3:Nx-2,3:Ny-2)=meanw(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1);
56			sumW(ii)=0.*sumW(ii);
57
58			for k=1:3
59			sumM1(:,:,k)=sumM(:,:);
60			sumD1(:,:,k)=sumD(:,:);
61			sumW1(:,:,k)=sumW(:,:);
62			end
63
64			for k=1:25
65
66			sumM1=smooth3(sumM1);
67			sumD1=smooth3(sumD1);
68			%sumW1=smooth3(sumW1);
69			end
70
71			V=[0.25*min(min(sumW(Nx1+2:Nx2,Ny1:Ny2-4))) 0];
72
73			% title='Vorticity input';
74			% imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
75			% set(gca,'ydir','norm')
76			% text(0,110,title);
77			%figure
78			v=zeros(10,1);
79			v1=zeros(10,1);
80			for i=1:10
81			%v(i)=-0.1(i-0.5)7;
82			%v1(i)=0.1(i-0.5)7;
83			v(i)=-0.01*(i);
84			v1(i)=0.01*(i);
85			end
86
87			% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v)
88			% hold on
89			% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
90			% hold off
91			% text(20,0,'countour interval is 0.7 in the non-dimensional units');
92			% text(20,95,'Vorticity input','Fontsize',16);
93
94			%figure
95
96			%V=[-0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))];
97			% title='Buoyancy diffusion integral';
98			% imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
99			% set(gca,'ydir','norm')
100			% text(0,110,title);
101
102
103			% figure
104			x=Nx1:Nx2;
105			y=Ny1:Ny2;
106			subplot(2,1,1)
107			contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v)
108			hold on
109			contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
110			text(10,Ny2+5,'Buoyancy eddy-transfer integral','Fontsize',13);
111			hold off
112			xlabel('X (gridpoints)')
113			ylabel('Y (gridpoints)')
114			set(gca,'DataAspectRatio',[2,2,2])
115
116			%figure
117			%V=[-0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))];
118			% title='Momentum eddy-transfer integral';
119			% imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
120			% set(gca,'ydir','norm')
121			% text(0,110,title);
122
123			subplot(2,1,2)
124			contour(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v)
125
126			set(gca,'DataAspectRatio',[1,1,1])
127			text(10,Ny2+5,'Momentum eddy-transfer integral','Fontsize',13);
128			hold on
129			contour(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
130			hold off
131			xlabel('X (gridpoints)')
132			ylabel('Y (gridpoints)')
133			text(40,-20,'countour interval - 0.01 (non-dimensional units)','FontSize',7);
134			% text(40,-30,'negative values - solid line','FontSize',7);
135			% text(40,-35,'positive values - dashed line','FontSize',7);
136
137
138
139			%-----TEST-------------------------------------------------
140
141			sum=0;
142			sumT=0;
143			sumY=0;
144			for i=Nx1:Nx2
145			for j=Ny1:Ny2
146			sum=sum+meanw(i,j,Nz1)pv(i,j,Nz1)dx*dy;
147			sumT=sum-meanw(i,j,Nz2)pv(i,j,Nz2)dx*dy;
148			end
149			end
150
151			for j=Ny1:Ny2
152			for k=Nz1:Nz2
153			sumT=sumT+meanu(Nx2,j,k)pv(Nx2,j,k)dz*dy;
154			sumT=sumT-meanu(Nx1,j,k)pv(Nx1,j,k)dz*dy;
155			end
156			end
157
158			for i=Nx1:Nx2
159			for k=Nz1:Nz2
160			sumT=sumT+meanv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
161			sumY=sumY+meanv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
162			sumT=sumT-meanv(i,Ny1,k)pv(i,Ny1,k)dz*dx;
163			end
164			end
165			'VORTICITY GENERATION'
166			sum
167			sumY
168			sumT
169
170
171			sum=0;
172			sumWW=0;
173			for i=Nx1:Nx2
174			for j=Ny1:Ny2
175			if i==Nx1
176			a=0.5;
177			elseif i==Nx2
178			a=0.5;
179			else
180			a=1;
181			end
182
183			if j==Ny1
184			b=0.5;
185			elseif j==Ny2
186			b=0.5;
187			else
188			b=1;
189			end
190
191			sum=sum+absumD(i,j)dxdy;
192			sumWW=sumWW+absumW(i,j)dxdy;
193			end
194			end
195			'DISSIPATION BY EDDY-DIFFUSIVITY'
196			sum
197			sumWW
198
199			sum=0;
200			for i=Nx1:Nx2
201			for j=Ny1:Ny2
202			if i==Nx1
203			a=0.5;
204			elseif i==Nx2
205			a=0.5;
206			else
207			a=1;
208			end
209
210			if j==Ny1
211			b=0.5;
212			elseif j==Ny2
213			b=0.5;
214			else
215			b=1;
216			end
217
218			sum=sum+absumM(i,j)dxdy;
219			end
220			end
221			'DISSIPATION BY EDDY-VISCOSITY'
222			sum