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