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