timour_matlab/mscripts/DIAGP2.m

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
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*Mp(3:Nx-2,3:Ny-2,k)*dz;
sumN(3:Nx-2,3:Ny-2)=sumN(3:Nx-2,3:Ny-2)+a*Np(3:Nx-2,3:Ny-2,k)*dz;
sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+a*Dp(3:Nx-2,3:Ny-2,k)*dz;
end
sumW(3:Nx-2,3:Ny-2)=meanw(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:25

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

for k=1:10
sumN1=smooth3(sumN1);
end

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

sumD1(4:Nx,Ny,1)=0;
sumM1(4:Nx,Ny,1)=0;
sumN1(4:Nx,Ny,1)=0;
sumD1(4:Nx,Ny-1,1)=0;
sumM1(4:Nx,Ny-1,1)=0;
sumN1(4:Nx,Ny-1,1)=0;


%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(4:Nx,1:Ny,1))',v,'k')
hold on
        contour(squeeze(sumD1(4:Nx,1:Ny,1))',v1,'k--')
        text(10,107,'Buoyancy transfer','Fontsize',14);
        text(-20,90,'a)','Fontsize',18);
hold off
 xlabel('X ')
 ylabel('Y ')
 set(gca,'XtickLabel','||')
 set(gca,'YtickLabel','||')

 set(gca,'DataAspectRatio',[1,1.6,2])

subplot(3,1,2)
        contour(squeeze(sumM1(4:Nx,1:Ny,1))',v,'k')

 set(gca,'DataAspectRatio',[1,1.6,1])
        text(10,107,'Momentum transfer','Fontsize',14);
        text(-20,90,'b)','Fontsize',18);
hold on
        contour(squeeze(sumM1(4:Nx,1:Ny,1))',v1,'k--')
hold off
 xlabel('X ')
 ylabel('Y ')
 set(gca,'XtickLabel','||')
 set(gca,'YtickLabel','||')


subplot(3,1,3)
        contour(squeeze(sumN1(4:Nx,1:Ny,1))',v,'k')

 set(gca,'DataAspectRatio',[1,1.6,1])
        text(10,107,'Inertial term ','Fontsize',14);
        text(-20,90,'c)','Fontsize',18);
hold on
        contour(squeeze(sumN1(4:Nx,1:Ny,1))',v1,'k--')
hold off
 xlabel('X ')
 ylabel('Y ')
 set(gca,'XtickLabel','||')
 set(gca,'YtickLabel','||')


%        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(1:Nx,1:Ny,1))',v)
hold on
        contour(squeeze(sumN1(1:Nx,1:Ny,1))',v1,'--')
hold off


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

sum=0;
sumT=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*meanw(i,j,Nz1)*pv(i,j,Nz1)*dx*dy;
sumT=sum-a*b*meanw(i,j,Nz2)*pv(i,j,Nz2)*dx*dy;
end
end

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

 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;
end
end

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

sumT=sumT+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;
end
end
'VORTICITY GENERATION'
sum
sumT


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

sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
sum=sum+sumN(i,j)*dx*dy;
end
end
'NONLINEARITY'
sum
1	edhill	1.1	sumM=zeros(Nx,Ny);
2			sumN=zeros(Nx,Ny);
3			sumD=zeros(Nx,Ny);
4			sumW=zeros(Nx,Ny);
5
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
17
18			for k=Nz1:Nz2
19			if k==Nz1
20			a=0.5;
21			elseif k==Nz2
22			a=0.5;
23			else
24			a=1;
25			end
26
27			sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+aMp(3:Nx-2,3:Ny-2,k)dz;
28			sumN(3:Nx-2,3:Ny-2)=sumN(3:Nx-2,3:Ny-2)+aNp(3:Nx-2,3:Ny-2,k)dz;
29			sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+aDp(3:Nx-2,3:Ny-2,k)dz;
30			end
31			sumW(3:Nx-2,3:Ny-2)=meanw(3:Nx-2,3:Ny-2,Nz1).*pv(3:Nx-2,3:Ny-2,Nz1);
32
33
34			for k=1:3
35			sumM1(:,:,k)=sumM(:,:);
36			sumN1(:,:,k)=sumN(:,:);
37			sumD1(:,:,k)=sumD(:,:);
38			sumW1(:,:,k)=sumW(:,:);
39			end
40
41			for k=1:25
42
43			sumM1=smooth3(sumM1);
44			sumD1=smooth3(sumD1);
45			sumW1=smooth3(sumW1);
46			end
47
48			for k=1:10
49			sumN1=smooth3(sumN1);
50			end
51
52			v=zeros(25,1);
53			v1=zeros(25,1);
54			for i=1:25
55			% v(i)=-0.1(i-0.5)7;
56			% v1(i)=0.1(i-0.5)7;
57			v(i)=0.001(i)10;
58			v1(i)=-0.001(i)10;
59			end
60
61			sumD1(4:Nx,Ny,1)=0;
62			sumM1(4:Nx,Ny,1)=0;
63			sumN1(4:Nx,Ny,1)=0;
64			sumD1(4:Nx,Ny-1,1)=0;
65			sumM1(4:Nx,Ny-1,1)=0;
66			sumN1(4:Nx,Ny-1,1)=0;
67
68
69
70			%V=[-10 10];
71			% title='Vorticity input';
72			% imagesc(lat,long,sumW1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
73			% set(gca,'ydir','norm')
74			% text(0,110,title);
75			%figure
76			% title='Buoyancy diffusion integral';
77			% imagesc(lat,long,sumD1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
78			% set(gca,'ydir','norm')
79			% text(0,110,title);
80			%figure
81			% title='Non-linear terms integral';
82			% imagesc(lat,long,sumN1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
83			% set(gca,'ydir','norm')
84			% text(0,110,title);
85			%figure
86			% title='Momentum diffusion integral';
87			% imagesc(lat,long,sumM1(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
88			% set(gca,'ydir','norm')
89			% text(0,110,title);
90			subplot(3,1,1)
91			contour(squeeze(sumD1(4:Nx,1:Ny,1))',v,'k')
92			hold on
93			contour(squeeze(sumD1(4:Nx,1:Ny,1))',v1,'k--')
94			text(10,107,'Buoyancy transfer','Fontsize',14);
95			text(-20,90,'a)','Fontsize',18);
96			hold off
97			xlabel('X ')
98			ylabel('Y ')
99			set(gca,'XtickLabel','\|\|')
100			set(gca,'YtickLabel','\|\|')
101
102			set(gca,'DataAspectRatio',[1,1.6,2])
103
104			subplot(3,1,2)
105			contour(squeeze(sumM1(4:Nx,1:Ny,1))',v,'k')
106
107			set(gca,'DataAspectRatio',[1,1.6,1])
108			text(10,107,'Momentum transfer','Fontsize',14);
109			text(-20,90,'b)','Fontsize',18);
110			hold on
111			contour(squeeze(sumM1(4:Nx,1:Ny,1))',v1,'k--')
112			hold off
113			xlabel('X ')
114			ylabel('Y ')
115			set(gca,'XtickLabel','\|\|')
116			set(gca,'YtickLabel','\|\|')
117
118
119
120			subplot(3,1,3)
121			contour(squeeze(sumN1(4:Nx,1:Ny,1))',v,'k')
122
123			set(gca,'DataAspectRatio',[1,1.6,1])
124			text(10,107,'Inertial term ','Fontsize',14);
125			text(-20,90,'c)','Fontsize',18);
126			hold on
127			contour(squeeze(sumN1(4:Nx,1:Ny,1))',v1,'k--')
128			hold off
129			xlabel('X ')
130			ylabel('Y ')
131			set(gca,'XtickLabel','\|\|')
132			set(gca,'YtickLabel','\|\|')
133
134
135			% text(40,-25,'countour interval - 1.5 (non-dimensional units)','FontSize',7);
136			% text(40,-30,'negative values - solid line','FontSize',7);
137			% text(40,-35,'positive values - dashed line','FontSize',7);
138			figure
139			contour(squeeze(sumN1(1:Nx,1:Ny,1))',v)
140			hold on
141			contour(squeeze(sumN1(1:Nx,1:Ny,1))',v1,'--')
142			hold off
143
144
145
146
147			%-----TEST-------------------------------------------------
148
149			sum=0;
150			sumT=0;
151			for i=Nx1:Nx2
152			for j=Ny1:Ny2
153			if i==Nx1
154			a=0.5;
155			elseif i==Nx2
156			a=0.5;
157			else
158			a=1;
159			end
160
161			if j==Ny1
162			b=0.5;
163			elseif j==Ny2
164			b=0.5;
165			else
166			b=1;
167			end
168
169			sum=sum+abmeanw(i,j,Nz1)pv(i,j,Nz1)dx*dy;
170			sumT=sum-abmeanw(i,j,Nz2)pv(i,j,Nz2)dx*dy;
171			end
172			end
173
174			for j=Ny1:Ny2
175			for k=Nz1:Nz2
176			if k==Nz1
177			a=0.5;
178			elseif k==Nz2
179			a=0.5;
180			else
181			a=1;
182			end
183
184			if j==Ny1
185			b=0.5;
186			elseif j==Ny2
187			b=0.5;
188			else
189			b=1;
190			end
191
192			sumT=sumT+abmeanu(Nx2,j,k)pv(Nx2,j,k)dz*dy;
193			sumT=sumT-abmeanu(Nx1,j,k)pv(Nx1,j,k)dz*dy;
194			end
195			end
196
197			for i=Nx1:Nx2
198			for k=Nz1:Nz2
199			if k==Nz1
200			a=0.5;
201			elseif k==Nz2
202			a=0.5;
203			else
204			a=1;
205			end
206
207			if i==Nx1
208			b=0.5;
209			elseif i==Nx2
210			b=0.5;
211			else
212			b=1;
213			end
214
215			sumT=sumT+abmeanv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
216			sumT=sumT-abmeanv(i,Ny1,k)pv(i,Ny1,k)dz*dx;
217			end
218			end
219			'VORTICITY GENERATION'
220			sum
221			sumT
222
223
224			sum=0;
225			for i=Nx1:Nx2
226			for j=Ny1:Ny2
227			if i==Nx1
228			a=0.5;
229			elseif i==Nx2
230			a=0.5;
231			else
232			a=1;
233			end
234
235			if j==Ny1
236			b=0.5;
237			elseif j==Ny2
238			b=0.5;
239			else
240			b=1;
241			end
242
243			sum=sum+absumD(i,j)dxdy;
244			end
245			end
246			'DISSIPATION BY EDDY-DIFFUSIVITY'
247			sum
248
249			sum=0;
250			for i=Nx1:Nx2
251			for j=Ny1:Ny2
252			if i==Nx1
253			a=0.5;
254			elseif i==Nx2
255			a=0.5;
256			else
257			a=1;
258			end
259
260			if j==Ny1
261			b=0.5;
262			elseif j==Ny2
263			b=0.5;
264			else
265			b=1;
266			end
267
268			sum=sum+absumM(i,j)dxdy;
269			end
270			end
271			'DISSIPATION BY EDDY-VISCOSITY'
272			sum
273
274			sum=0;
275			for i=Nx1:Nx2
276			for j=Ny1:Ny2
277			sum=sum+sumN(i,j)dxdy;
278			end
279			end
280			'NONLINEARITY'
281			sum