timour_matlab/mscripts/DIAGtest.m

% 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
1	edhill	1.1	% PV diagnostics 'DIAGtest'
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			% for single gyre
14			% Ny2=Ny-3;
15			% for double-gyre
16			Ny2=45;
17
18			Nz1=5;
19			Nz2=Nz-4;
20
21			%Nx1=30;
22			%Nx2=32;
23			%Ny1=30;
24			%Ny2=32;
25			%Nz1=10;
26			%Nz2=12;
27
28			pv1=pv;
29			Mp1=Mp;
30			Dp1=Dp;
31
32			% for kk=1:4
33			% pv1=smooth3(pv1);
34			% Mp1=smooth3(Mp1);
35			% Dp1=smooth3(Dp1);
36			% end
37
38
39			for k=Nz1:Nz2
40
41			if k==Nz1
42			a=0.5;
43			elseif k==Nz2
44			a=0.5;
45			else
46			a=1;
47			end
48
49			sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+aMp1(3:Nx-2,3:Ny-2,k)dz;
50			sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+aDp1(3:Nx-2,3:Ny-2,k)dz;
51			end
52			sumW(3:Nx-2,3:Ny-2)=meanw(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1);
53
54			for k=1:3
55			sumM1(:,:,k)=sumM(:,:);
56			sumD1(:,:,k)=sumD(:,:);
57			sumW1(:,:,k)=sumW(:,:);
58			end
59
60			for k=1:15
61
62			sumM1=smooth3(sumM1);
63			sumD1=smooth3(sumD1);
64			sumW1=smooth3(sumW1);
65			end
66
67			V=[0.25*min(min(sumW(Nx1+2:Nx2,Ny1:Ny2-4))) 0];
68
69			% title='Vorticity input';
70			% imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
71			% set(gca,'ydir','norm')
72			% text(0,110,title);
73			%figure
74			v=zeros(10,1);
75			v1=zeros(10,1);
76			for i=1:10
77			%v(i)=-0.1(i-0.5)7;
78			%v1(i)=0.1(i-0.5)7;
79			v(i)=-0.01*(i);
80			v1(i)=0.01*(i);
81			end
82
83			% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v)
84			% hold on
85			% contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
86			% hold off
87			% text(20,0,'countour interval is 0.7 in the non-dimensional units');
88			% text(20,95,'Vorticity input','Fontsize',16);
89
90			%figure
91
92			%V=[-0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))];
93			% title='Buoyancy diffusion integral';
94			% imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
95			% set(gca,'ydir','norm')
96			% text(0,110,title);
97
98
99			% figure
100			x=Nx1:Nx2;
101			y=Ny1:Ny-2;
102			subplot(2,1,1)
103			contourf(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny-2,1))',10);colorbar;
104			% contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v)
105			%hold on
106			% contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
107			text(10,Ny2+5,'Buoyancy eddy-transfer','Fontsize',13);
108			%hold off
109			xlabel('X (gridpoints)')
110			ylabel('Y (gridpoints)')
111			set(gca,'DataAspectRatio',[2,2,2])
112
113			%figure
114			%V=[-0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))];
115			% title='Momentum eddy-transfer';
116			% imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
117			% set(gca,'ydir','norm')
118			% text(0,110,title);
119
120			subplot(2,1,2)
121			contourf(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny-2,1))',10);colorbar;
122			% contour(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v)
123
124			set(gca,'DataAspectRatio',[1,1,1])
125			text(10,Ny2+5,'Momentum eddy-transfer integral','Fontsize',13);
126			%hold on
127			% contour(x,y,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
128			%hold off
129			xlabel('X (gridpoints)')
130			ylabel('Y (gridpoints)')
131			% text(40,-20,'countour interval - 0.01 (non-dimensional units)','FontSize',7);
132			% text(40,-30,'negative values - solid line','FontSize',7);
133			% text(40,-35,'positive values - dashed line','FontSize',7);
134
135
136
137			%-----TEST-------------------------------------------------
138
139			dz
140
141			bx=zeros(Nx,Ny,Nz);
142			by=zeros(Nx,Ny,Nz);
143			bz=zeros(Nx,Ny,Nz);
144			zetax=zeros(Nx,Ny,Nz);
145			zetay=zeros(Nx,Ny,Nz);
146			zetaz=zeros(Nx,Ny,Nz);
147
148			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);
149			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);
150			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));
151
152			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));
153			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));
154			zetaz(2:Nx-1,2:Ny-1,2:Nz-1)=(meanv(3:Nx,2:Ny-1,2:Nz-1)/(2dx)-meanv(1:Nx-2,2:Ny-1,2:Nz-1)/(2dx) ...
155			-meanu(2:Nx-1,3:Ny,2:Nz-1)/(2dy)+meanu(2:Nx-1,1:Ny-2,2:Nz-1)/(2dy) ...
156			+ff(2:Nx-1,2:Ny-1,2:Nz-1));
157
158
159			sum=0;
160			sumT=0;
161			sumV=0;
162			sumY=0;
163			summ=0;
164			sumd=0;
165			for i=Nx1:Nx2
166			for j=Ny1:Ny2
167
168			if i==Nx1
169			a=0.5;
170			elseif i==Nx2
171			a=0.5;
172			else
173			a=1;
174			end
175
176			if j==Ny1
177			b=0.5;
178			elseif j==Ny2
179			b=0.5;
180			else
181			b=1;
182			end
183
184
185			summ=summ+ab(Mu(i,j,Nz1)by(i,j,Nz1)-Mv(i,j,Nz1)bx(i,j,Nz1))dxdy;
186			summ=summ-ab(Mu(i,j,Nz2)by(i,j,Nz2)-Mv(i,j,Nz2)bx(i,j,Nz2))dxdy;
187			sumd=sumd+abzetaz(i,j,Nz1)D(i,j,Nz1)dx*dy;
188			sumd=sumd-abzetaz(i,j,Nz2)D(i,j,Nz2)dx*dy;
189			sum=sum+abmeanw(i,j,Nz1)pv(i,j,Nz1)dx*dy;
190			sumT=sum-abmeanw(i,j,Nz2)pv(i,j,Nz2)dx*dy;
191			sumV=sumV+abmeanw(i,j,Nz1)dxdy;
192			sumV=sumV-abmeanw(i,j,Nz2)dxdy;
193			end
194			end
195
196			'xy'
197			summ
198			sumV
199			sumd
200
201			for j=Ny1:Ny2
202			for k=Nz1:Nz2
203
204			if k==Nz1
205			a=0.5;
206			elseif k==Nz2
207			a=0.5;
208			else
209			a=1;
210			end
211
212			if j==Ny1
213			b=0.5;
214			elseif j==Ny2
215			b=0.5;
216			else
217			b=1;
218			end
219
220
221
222			summ=summ+abMv(Nx2,j,k)bz(Nx2,j,k)dz*dy;
223			summ=summ-abMv(Nx1,j,k)bz(Nx1,j,k)dz*dy;
224			sumd=sumd+abzetax(Nx2,j,k)D(Nx2,j,k)dz*dy;
225			sumd=sumd-abzetax(Nx1,j,k)D(Nx1,j,k)dz*dy;
226			sumT=sumT+abmeanu(Nx2,j,k)pv(Nx2,j,k)dz*dy;
227			sumT=sumT-abmeanu(Nx1,j,k)pv(Nx1,j,k)dz*dy;
228			sumV=sumV+abmeanu(Nx2,j,k)dzdy;
229			sumV=sumV-abmeanu(Nx1,j,k)dzdy;
230			end
231			end
232
233			'xy+yz'
234			summ
235			sumd
236
237			for i=Nx1:Nx2
238			for k=Nz1:Nz2
239
240			if k==Nz1
241			a=0.5;
242			elseif k==Nz2
243			a=0.5;
244			else
245			a=1;
246			end
247
248			if i==Nx1
249			b=0.5;
250			elseif i==Nx2
251			b=0.5;
252			else
253			b=1;
254			end
255
256
257			summ=summ-abMu(i,Ny2,k)bz(i,Ny2,k)dz*dx;
258			summ=summ+abMu(i,Ny1,k)bz(i,Ny1,k)dz*dx;
259			sumd=sumd+abzetay(i,Ny2,k)D(i,Ny2,k)dz*dx;
260			sumd=sumd-abzetay(i,Ny1,k)D(i,Ny1,k)dz*dx;
261			sumT=sumT+abmeanv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
262			sumY=sumY+abmeanv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
263			sumT=sumT-abmeanv(i,Ny1,k)pv(i,Ny1,k)dz*dx;
264			sumV=sumV+abmeanv(i,Ny2,k)dzdx;
265			sumV=sumV-abmeanv(i,Ny1,k)dzdx;
266			end
267			end
268			'TOTAL'
269			sum
270			sumY
271			sumT
272			sumV
273			summ
274			sumd
275
276
277			sum=0;
278			for i=Nx1:Nx2
279			for j=Ny1:Ny2
280
281			if i==Nx1
282			a=0.5;
283			elseif i==Nx2
284			a=0.5;
285			else
286			a=1;
287			end
288
289			if j==Ny1
290			b=0.5;
291			elseif j==Ny2
292			b=0.5;
293			else
294			b=1;
295			end
296
297			sum=sum+absumD(i,j)dxdy;
298			end
299			end
300			'DISSIPATION BY EDDY-DIFFUSIVITY'
301			sum
302
303			sum=0;
304			for i=Nx1:Nx2
305			for j=Ny1:Ny2
306
307			if i==Nx1
308			a=0.5;
309			elseif i==Nx2
310			a=0.5;
311			else
312			a=1;
313			end
314
315			if j==Ny1
316			b=0.5;
317			elseif j==Ny2
318			b=0.5;
319			else
320			b=1;
321			end
322
323			sum=sum+absumM(i,j)dxdy;
324			end
325			end
326			'DISSIPATION BY EDDY-VISCOSITY'
327			sum
328
329			'IMBALANCE'
330			(sumT-(summ+sumd))/sumd