timour_matlab/mscripts/DIAG3.m

% Enstrophy 'DIAG3'
sumM=zeros(Nx,Ny);
sumD=zeros(Nx,Ny);
sumW=zeros(Nx,Ny);

dz=0.005;

% DEFINE BOX OF INTEGRATION
Nx1=140;
Nx2=Nx-4;
Ny1=5;
Ny2=Ny-4;
Nz1=3;
Nz2=Nz/2-2;

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

%  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
 sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+Mp1(3:Nx-2,3:Ny-2,k).*pv1(3:Nx-2,3:Ny-2,k)*dz;
 sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+Dp1(3:Nx-2,3:Ny-2,k).*pv1(3:Nx-2,3:Ny-2,k)*dz;
 end
 sumW(3:Nx-2,3:Ny-2)=0.5*meanw(3:Nx-2,3:Ny-2,Nz1).*pv1(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: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.1*(i)*7;
v1(i)=0.1*(i)*7;
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

subplot(2,1,1)
        contour(squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',50)
        text(10,Ny2,'Buoyancy diffusion integral','Fontsize',15);
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 diffusion 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(squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',50)

 set(gca,'DataAspectRatio',[1,1,1])
        text(10,Ny2,'Momentum diffusion integral','Fontsize',15);
 xlabel('X (gridpoints)')
 ylabel('Y (gridpoints)')
        text(40,-25,'countour interval - 0.7 (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+0.5*meanw(i,j,Nz1)*pv(i,j,Nz1)*pv(i,j,Nz1)*dx*dy;
sumT=sum-0.5*meanw(i,j,Nz2)*pv(i,j,Nz2)*pv(i,j,Nz2)*dx*dy;
end
end

for j=Ny1:Ny2
for k=Nz1:Nz2
 sumT=sumT+0.5*meanu(Nx2,j,k)*pv(Nx2,j,k)*pv(Nx2,j,k)*dz*dy;
 sumT=sumT-0.5*meanu(Nx1,j,k)*pv(Nx1,j,k)*pv(Nx1,j,k)*dz*dy;
end
end

for i=Nx1:Nx2
for k=Nz1:Nz2
sumT=sumT+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
sumY=sumY+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
sumT=sumT-0.5*meanv(i,Ny1,k)*pv(i,Ny1,k)*pv(i,Ny1,k)*dz*dx;
end
end
'VORTICITY GENERATION'
sum
sumY
sumT


sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
sum=sum+sumD(i,j)*dx*dy;
end
end
'DISSIPATION BY EDDY-DIFFUSIVITY'
sum

sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
sum=sum+sumM(i,j)*dx*dy;
end
end
'DISSIPATION BY EDDY-VISCOSITY'
sum
1	edhill	1.1	% Enstrophy 'DIAG3'
2			sumM=zeros(Nx,Ny);
3			sumD=zeros(Nx,Ny);
4			sumW=zeros(Nx,Ny);
5
6			dz=0.005;
7
8			% DEFINE BOX OF INTEGRATION
9			Nx1=140;
10			Nx2=Nx-4;
11			Ny1=5;
12			Ny2=Ny-4;
13			Nz1=3;
14			Nz2=Nz/2-2;
15
16			pv1=pv;
17			Mp1=Mp;
18			Dp1=Dp;
19
20			% pv1=smooth3(pv);
21			% Mp1=smooth3(Mp);
22			% Dp1=smooth3(Dp);
23
24			% for kk=1:4
25			% pv=smooth3(pv);
26			% Mp=smooth3(Mp);
27			% Dp=smooth3(Dp);
28			% end
29
30			for k=Nz1:Nz2
31			sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+Mp1(3:Nx-2,3:Ny-2,k).pv1(3:Nx-2,3:Ny-2,k)dz;
32			sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+Dp1(3:Nx-2,3:Ny-2,k).pv1(3:Nx-2,3:Ny-2,k)dz;
33			end
34			sumW(3:Nx-2,3:Ny-2)=0.5meanw(3:Nx-2,3:Ny-2,Nz1).pv1(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1);
35
36			for k=1:3
37			sumM1(:,:,k)=sumM(:,:);
38			sumD1(:,:,k)=sumD(:,:);
39			sumW1(:,:,k)=sumW(:,:);
40			end
41
42			for k=1:25
43
44			sumM1=smooth3(sumM1);
45			sumD1=smooth3(sumD1);
46			sumW1=smooth3(sumW1);
47			end
48
49			V=[0.25*min(min(sumW(Nx1+2:Nx2,Ny1:Ny2-4))) 0];
50
51			% title='Vorticity input';
52			% imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
53			% set(gca,'ydir','norm')
54			% text(0,110,title);
55			%figure
56			v=zeros(10,1);
57			v1=zeros(10,1);
58			for i=1:10
59			%v(i)=-0.1(i-0.5)7;
60			%v1(i)=0.1(i-0.5)7;
61			v(i)=-0.1(i)7;
62			v1(i)=0.1(i)7;
63			end
64
65			% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v)
66			% hold on
67			% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
68			% hold off
69			% text(20,0,'countour interval is 0.7 in the non-dimensional units');
70			% text(20,95,'Vorticity input','Fontsize',16);
71
72			%figure
73
74			%V=[-0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))];
75			% title='Buoyancy diffusion integral';
76			% imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
77			% set(gca,'ydir','norm')
78			% text(0,110,title);
79
80
81			% figure
82
83			subplot(2,1,1)
84			contour(squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',50)
85			text(10,Ny2,'Buoyancy diffusion integral','Fontsize',15);
86			hold off
87			xlabel('X (gridpoints)')
88			ylabel('Y (gridpoints)')
89			set(gca,'DataAspectRatio',[2,2,2])
90
91			%figure
92			%V=[-0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))];
93			% title='Momentum diffusion integral';
94			% imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
95			% set(gca,'ydir','norm')
96			% text(0,110,title);
97
98			subplot(2,1,2)
99			contour(squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',50)
100
101			set(gca,'DataAspectRatio',[1,1,1])
102			text(10,Ny2,'Momentum diffusion integral','Fontsize',15);
103			xlabel('X (gridpoints)')
104			ylabel('Y (gridpoints)')
105			text(40,-25,'countour interval - 0.7 (non-dimensional units)','FontSize',7);
106			text(40,-30,'negative values - solid line','FontSize',7);
107			text(40,-35,'positive values - dashed line','FontSize',7);
108
109
110
111			%-----TEST-------------------------------------------------
112
113			sum=0;
114			sumT=0;
115			sumY=0;
116			for i=Nx1:Nx2
117			for j=Ny1:Ny2
118			sum=sum+0.5meanw(i,j,Nz1)pv(i,j,Nz1)pv(i,j,Nz1)dx*dy;
119			sumT=sum-0.5meanw(i,j,Nz2)pv(i,j,Nz2)pv(i,j,Nz2)dx*dy;
120			end
121			end
122
123			for j=Ny1:Ny2
124			for k=Nz1:Nz2
125			sumT=sumT+0.5meanu(Nx2,j,k)pv(Nx2,j,k)pv(Nx2,j,k)dz*dy;
126			sumT=sumT-0.5meanu(Nx1,j,k)pv(Nx1,j,k)pv(Nx1,j,k)dz*dy;
127			end
128			end
129
130			for i=Nx1:Nx2
131			for k=Nz1:Nz2
132			sumT=sumT+0.5meanv(i,Ny2,k)pv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
133			sumY=sumY+0.5meanv(i,Ny2,k)pv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
134			sumT=sumT-0.5meanv(i,Ny1,k)pv(i,Ny1,k)pv(i,Ny1,k)dz*dx;
135			end
136			end
137			'VORTICITY GENERATION'
138			sum
139			sumY
140			sumT
141
142
143			sum=0;
144			for i=Nx1:Nx2
145			for j=Ny1:Ny2
146			sum=sum+sumD(i,j)dxdy;
147			end
148			end
149			'DISSIPATION BY EDDY-DIFFUSIVITY'
150			sum
151
152			sum=0;
153			for i=Nx1:Nx2
154			for j=Ny1:Ny2
155			sum=sum+sumM(i,j)dxdy;
156			end
157			end
158			'DISSIPATION BY EDDY-VISCOSITY'
159			sum