timour_matlab/mscripts/DIAG2.m

clear path 

global Nx Ny Nz
global lat long dz dm mdep
global delt_su su_its t_su delt
global descriptor this_path
global f deltaf Q beta r_expt r_heat H
global time rots it
global g Cp rho_bar alpha
global u v t w 
global iterations


param_file_name =   ...
   input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
feval(param_file_name) ;

% iterations

itstart = input(' Please enter start iteration : ','s') 
itend = input(' Please enter end iteration : ','s') 


sizeit=size(iterations);
for i=1:sizeit(1)
iter(i)=eval(iterations(i,1:10));
end
nitstart=find(iter==eval(itstart))
nitend=find(iter==eval(itend))

path   = this_path 
cmdstr=['cd ' path ];
eval(cmdstr);
path=pwd

sumtheta=zeros(Nx,Ny,Nz);
sumu=zeros(Nx,Ny,Nz);
sumv=zeros(Nx,Ny,Nz);
sumw=zeros(Nx,Ny,Nz);
counter=0;

for i=nitstart:nitend
i,0
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');

ufilename=(['U.' iterations((i),1:10) ]) ;
u=rdmds(ufilename,'b');
u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));

vfilename=(['V.' iterations((i),1:10) ]) ;
v=rdmds(vfilename,'b');
v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));

wfilename=(['W.' iterations((i),1:10) ]) ;
w=rdmds(wfilename,'b');

sumtheta=sumtheta+t;
sumu=sumu+u;
sumv=sumv+v;
sumw=sumw+w;
counter=counter+1;
end

meantheta=sumtheta/counter;
meanu=sumu/counter;
meanv=sumv/counter;
meanw=sumw/counter;

sumut=zeros(Nx,Ny,Nz);
sumvt=zeros(Nx,Ny,Nz);
sumwt=zeros(Nx,Ny,Nz);

sumuu=zeros(Nx,Ny,Nz);
sumvu=zeros(Nx,Ny,Nz);
sumwu=zeros(Nx,Ny,Nz);
sumvv=zeros(Nx,Ny,Nz);
sumwv=zeros(Nx,Ny,Nz);

%meanu=smooth3(meanu);
%meanv=smooth3(meanv);
%meanw=smooth3(meanw);
%meantheta=smooth3(meantheta);

for i=nitstart:nitend
i
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');

ufilename=(['U.' iterations((i),1:10) ]) ;
u=rdmds(ufilename,'b');
u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));

vfilename=(['V.' iterations((i),1:10) ]) ;
v=rdmds(vfilename,'b');
v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));

wfilename=(['W.' iterations((i),1:10) ]) ;
w=rdmds(wfilename,'b');

t=t-meantheta;
u=u-meanu;
v=v-meanv;
w=w-meanw;

sumut=sumut+u.*t;
sumvt=sumvt+v.*t;
sumwt=sumwt+w.*t;

sumuu=sumuu+u.*u;
sumvu=sumvu+v.*u;
sumwu=sumwu+w.*u;
sumvv=sumvv+v.*v;
sumwv=sumwv+w.*v;
end

sumut=sumut/counter;
sumvt=sumvt/counter;
sumwt=sumwt/counter;

sumuu=sumuu/counter;
sumvu=sumvu/counter;
sumwu=sumwu/counter;
sumvv=sumvv/counter;
sumwv=sumwv/counter;

Mu=zeros(Nx,Ny,Nz);
Mv=zeros(Nx,Ny,Nz);
D=zeros(Nx,Ny,Nz);

dx=dm;
dy=dm;
dz=-dz;

Mu(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumuu(3:Nx,2:Ny-1,2:Nz-1)-sumuu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
        -(sumvu(2:Nx-1,3:Ny,2:Nz-1)-sumvu(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
        -(sumwu(2:Nx-1,2:Ny-1,3:Nz)-sumwu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);

Mv(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumvu(3:Nx,2:Ny-1,2:Nz-1)-sumvu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
        -(sumvv(2:Nx-1,3:Ny,2:Nz-1)-sumvv(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
        -(sumwv(2:Nx-1,2:Ny-1,3:Nz)-sumwv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);

D(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumut(3:Nx,2:Ny-1,2:Nz-1)-sumut(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
        -(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
        -(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);

pv=zeros(Nx,Ny,Nz);
ff=zeros(Nx,Ny,Nz);


for j=1:Ny
ff(:,j,:)=f+(j-1)*dy*beta;
end

pv(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)) ...
                    .*(meantheta(2:Nx-1,2:Ny-1,3:Nz)-meantheta(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
                -(meanv(2:Nx-1,2:Ny-1,3:Nz)-meanv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
                .*(meantheta(3:Nx,2:Ny-1,2:Nz-1)-meantheta(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx)+ ...
                (meanu(2:Nx-1,2:Ny-1,3:Nz)-meanu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
                .*(meantheta(2:Nx-1,3:Ny,2:Nz-1)-meantheta(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy);

Mp=zeros(Nx,Ny,Nz);
Dp=zeros(Nx,Ny,Nz);

Mp(3:Nx-2,3:Ny-2,3:Nz-2)=((Mv(4:Nx-1,3:Ny-2,3:Nz-2)-Mv(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx)-  ...
(Mu(3:Nx-2,4:Ny-1,3:Nz-2)-Mu(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy)).* ...
(meantheta(3:Nx-2,3:Ny-2,4:Nz-1)-meantheta(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz) ...
        -(Mv(3:Nx-2,3:Ny-2,4:Nz-1)-Mv(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz).* ...
        (meantheta(4:Nx-1,3:Ny-2,3:Nz-2)-meantheta(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx) ...
+(Mu(3:Nx-2,3:Ny-2,4:Nz-1)-Mu(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz).* ...
(meantheta(3:Nx-2,4:Ny-1,3:Nz-2)-meantheta(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy);

Dp(3:Nx-2,3:Ny-2,3:Nz-2)=(D(3:Nx-2,3:Ny-2,4:Nz-1)-D(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz).* ...
(meanv(4:Nx-1,3:Ny-2,3:Nz-2)/(2*dx)-meanv(2:Nx-3,3:Ny-2,3:Nz-2)/(2*dx) ...
-meanu(3:Nx-2,4:Ny-1,3:Nz-2)/(2*dy)+meanu(3:Nx-2,2:Ny-3,3:Nz-2)/(2*dy)+ff(3:Nx-2,3:Ny-2,3:Nz-2)) ...
        -(D(4:Nx-1,3:Ny-2,3:Nz-2)-D(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx).*  ...
        (meanv(3:Nx-2,3:Ny-2,4:Nz-1)-meanv(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz)+  ...
(D(3:Nx-2,4:Ny-1,3:Nz-2)-D(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy).* ...
(meanu(3:Nx-2,3:Ny-2,4:Nz-1)-meanu(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz);


sumM=zeros(Nx,Ny);
sumD=zeros(Nx,Ny);
sumW=zeros(Nx,Ny);

dz=-dz;

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


for k=Nz1:Nz2
sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+Mp(3:Nx-2,3:Ny-2,k).*pv(3:Nx-2,3:Ny-2,k)*dz;
sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+Dp(3:Nx-2,3:Ny-2,k).*pv(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).*pv(3:Nx-2,3:Ny-2,Nz1).*pv(3:Nx-2,3:Ny-2,Nz1);

%contourf(sumM',20);colorbar;
%figure
%contourf(sumD',20);colorbar;

V=[0.25*min(min(sumW(Nx1:Nx2,Ny1:Ny2-4))) 0];

        title='Vorticity input';
        imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
        set(gca,'ydir','norm')
        text(0,110,title);
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 eddy-transfer integral';
        imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
        set(gca,'ydir','norm')
        text(0,110,title);
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 integral';
        imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
        set(gca,'ydir','norm')
        text(0,110,title);


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

sum=0;
sumT=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;
sumT=sumT-0.5*meanv(i,Ny1,k)*pv(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
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	clear path
2
3			global Nx Ny Nz
4			global lat long dz dm mdep
5			global delt_su su_its t_su delt
6			global descriptor this_path
7			global f deltaf Q beta r_expt r_heat H
8			global time rots it
9			global g Cp rho_bar alpha
10			global u v t w
11			global iterations
12
13
14			param_file_name = ...
15			input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
16			feval(param_file_name) ;
17
18			% iterations
19
20			itstart = input(' Please enter start iteration : ','s')
21			itend = input(' Please enter end iteration : ','s')
22
23
24			sizeit=size(iterations);
25			for i=1:sizeit(1)
26			iter(i)=eval(iterations(i,1:10));
27			end
28			nitstart=find(iter==eval(itstart))
29			nitend=find(iter==eval(itend))
30
31			path = this_path
32			cmdstr=['cd ' path ];
33			eval(cmdstr);
34			path=pwd
35
36			sumtheta=zeros(Nx,Ny,Nz);
37			sumu=zeros(Nx,Ny,Nz);
38			sumv=zeros(Nx,Ny,Nz);
39			sumw=zeros(Nx,Ny,Nz);
40			counter=0;
41
42			for i=nitstart:nitend
43			i,0
44			tfilename=(['T.' iterations((i),1:10) ]) ;
45			t=rdmds(tfilename,'b');
46
47			ufilename=(['U.' iterations((i),1:10) ]) ;
48			u=rdmds(ufilename,'b');
49			u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
50
51			vfilename=(['V.' iterations((i),1:10) ]) ;
52			v=rdmds(vfilename,'b');
53			v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
54
55			wfilename=(['W.' iterations((i),1:10) ]) ;
56			w=rdmds(wfilename,'b');
57
58			sumtheta=sumtheta+t;
59			sumu=sumu+u;
60			sumv=sumv+v;
61			sumw=sumw+w;
62			counter=counter+1;
63			end
64
65			meantheta=sumtheta/counter;
66			meanu=sumu/counter;
67			meanv=sumv/counter;
68			meanw=sumw/counter;
69
70			sumut=zeros(Nx,Ny,Nz);
71			sumvt=zeros(Nx,Ny,Nz);
72			sumwt=zeros(Nx,Ny,Nz);
73
74			sumuu=zeros(Nx,Ny,Nz);
75			sumvu=zeros(Nx,Ny,Nz);
76			sumwu=zeros(Nx,Ny,Nz);
77			sumvv=zeros(Nx,Ny,Nz);
78			sumwv=zeros(Nx,Ny,Nz);
79
80			%meanu=smooth3(meanu);
81			%meanv=smooth3(meanv);
82			%meanw=smooth3(meanw);
83			%meantheta=smooth3(meantheta);
84
85			for i=nitstart:nitend
86			i
87			tfilename=(['T.' iterations((i),1:10) ]) ;
88			t=rdmds(tfilename,'b');
89
90			ufilename=(['U.' iterations((i),1:10) ]) ;
91			u=rdmds(ufilename,'b');
92			u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:));
93
94			vfilename=(['V.' iterations((i),1:10) ]) ;
95			v=rdmds(vfilename,'b');
96			v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:));
97
98			wfilename=(['W.' iterations((i),1:10) ]) ;
99			w=rdmds(wfilename,'b');
100
101			t=t-meantheta;
102			u=u-meanu;
103			v=v-meanv;
104			w=w-meanw;
105
106			sumut=sumut+u.*t;
107			sumvt=sumvt+v.*t;
108			sumwt=sumwt+w.*t;
109
110			sumuu=sumuu+u.*u;
111			sumvu=sumvu+v.*u;
112			sumwu=sumwu+w.*u;
113			sumvv=sumvv+v.*v;
114			sumwv=sumwv+w.*v;
115			end
116
117			sumut=sumut/counter;
118			sumvt=sumvt/counter;
119			sumwt=sumwt/counter;
120
121			sumuu=sumuu/counter;
122			sumvu=sumvu/counter;
123			sumwu=sumwu/counter;
124			sumvv=sumvv/counter;
125			sumwv=sumwv/counter;
126
127			Mu=zeros(Nx,Ny,Nz);
128			Mv=zeros(Nx,Ny,Nz);
129			D=zeros(Nx,Ny,Nz);
130
131			dx=dm;
132			dy=dm;
133			dz=-dz;
134
135			Mu(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumuu(3:Nx,2:Ny-1,2:Nz-1)-sumuu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
136			-(sumvu(2:Nx-1,3:Ny,2:Nz-1)-sumvu(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
137			-(sumwu(2:Nx-1,2:Ny-1,3:Nz)-sumwu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
138
139			Mv(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumvu(3:Nx,2:Ny-1,2:Nz-1)-sumvu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
140			-(sumvv(2:Nx-1,3:Ny,2:Nz-1)-sumvv(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
141			-(sumwv(2:Nx-1,2:Ny-1,3:Nz)-sumwv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
142
143			D(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumut(3:Nx,2:Ny-1,2:Nz-1)-sumut(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
144			-(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
145			-(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
146
147			pv=zeros(Nx,Ny,Nz);
148			ff=zeros(Nx,Ny,Nz);
149
150
151			for j=1:Ny
152			ff(:,j,:)=f+(j-1)dybeta;
153			end
154
155			pv(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) ...
156			-meanu(2:Nx-1,3:Ny,2:Nz-1)/(2dy)+meanu(2:Nx-1,1:Ny-2,2:Nz-1)/(2dy) ...
157			+ff(2:Nx-1,2:Ny-1,2:Nz-1)) ...
158			.(meantheta(2:Nx-1,2:Ny-1,3:Nz)-meantheta(2:Nx-1,2:Ny-1,1:Nz-2))/(2dz) ...
159			-(meanv(2:Nx-1,2:Ny-1,3:Nz)-meanv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
160			.(meantheta(3:Nx,2:Ny-1,2:Nz-1)-meantheta(1:Nx-2,2:Ny-1,2:Nz-1))/(2dx)+ ...
161			(meanu(2:Nx-1,2:Ny-1,3:Nz)-meanu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz) ...
162			.(meantheta(2:Nx-1,3:Ny,2:Nz-1)-meantheta(2:Nx-1,1:Ny-2,2:Nz-1))/(2dy);
163
164			Mp=zeros(Nx,Ny,Nz);
165			Dp=zeros(Nx,Ny,Nz);
166
167			Mp(3:Nx-2,3:Ny-2,3:Nz-2)=((Mv(4:Nx-1,3:Ny-2,3:Nz-2)-Mv(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx)- ...
168			(Mu(3:Nx-2,4:Ny-1,3:Nz-2)-Mu(3:Nx-2,2:Ny-3,3:Nz-2))/(2dy)). ...
169			(meantheta(3:Nx-2,3:Ny-2,4:Nz-1)-meantheta(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz) ...
170			-(Mv(3:Nx-2,3:Ny-2,4:Nz-1)-Mv(3:Nx-2,3:Ny-2,2:Nz-3))/(2dz). ...
171			(meantheta(4:Nx-1,3:Ny-2,3:Nz-2)-meantheta(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx) ...
172			+(Mu(3:Nx-2,3:Ny-2,4:Nz-1)-Mu(3:Nx-2,3:Ny-2,2:Nz-3))/(2dz). ...
173			(meantheta(3:Nx-2,4:Ny-1,3:Nz-2)-meantheta(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy);
174
175			Dp(3:Nx-2,3:Ny-2,3:Nz-2)=(D(3:Nx-2,3:Ny-2,4:Nz-1)-D(3:Nx-2,3:Ny-2,2:Nz-3))/(2dz). ...
176			(meanv(4:Nx-1,3:Ny-2,3:Nz-2)/(2dx)-meanv(2:Nx-3,3:Ny-2,3:Nz-2)/(2dx) ...
177			-meanu(3:Nx-2,4:Ny-1,3:Nz-2)/(2dy)+meanu(3:Nx-2,2:Ny-3,3:Nz-2)/(2dy)+ff(3:Nx-2,3:Ny-2,3:Nz-2)) ...
178			-(D(4:Nx-1,3:Ny-2,3:Nz-2)-D(2:Nx-3,3:Ny-2,3:Nz-2))/(2dx). ...
179			(meanv(3:Nx-2,3:Ny-2,4:Nz-1)-meanv(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz)+ ...
180			(D(3:Nx-2,4:Ny-1,3:Nz-2)-D(3:Nx-2,2:Ny-3,3:Nz-2))/(2dy). ...
181			(meanu(3:Nx-2,3:Ny-2,4:Nz-1)-meanu(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz);
182
183
184			sumM=zeros(Nx,Ny);
185			sumD=zeros(Nx,Ny);
186			sumW=zeros(Nx,Ny);
187
188			dz=-dz;
189
190			% DEFINE BOX OF INTEGRATION
191			Nx1=3;
192			Nx2=Nx-2;
193			Ny1=3;
194			%Ny2=Ny-2;
195			Ny2=Ny/2-2;
196			Nz1=5;
197			Nz2=Nz-4;
198
199
200			for k=Nz1:Nz2
201			sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+Mp(3:Nx-2,3:Ny-2,k).pv(3:Nx-2,3:Ny-2,k)dz;
202			sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+Dp(3:Nx-2,3:Ny-2,k).pv(3:Nx-2,3:Ny-2,k)dz;
203			end
204			sumW(3:Nx-2,3:Ny-2)=0.5meanw(3:Nx-2,3:Ny-2,Nz1).pv(3:Nx-2,3:Ny-2,Nz1).*pv(3:Nx-2,3:Ny-2,Nz1);
205
206			%contourf(sumM',20);colorbar;
207			%figure
208			%contourf(sumD',20);colorbar;
209
210			V=[0.25*min(min(sumW(Nx1:Nx2,Ny1:Ny2-4))) 0];
211
212			title='Vorticity input';
213			imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
214			set(gca,'ydir','norm')
215			text(0,110,title);
216			figure
217			V=[-0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))]
218			title='Buoyancy eddy-transfer integral';
219			imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
220			set(gca,'ydir','norm')
221			text(0,110,title);
222			figure
223			V=[-0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))]
224			title='Momentum eddy-transfer integral';
225			imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2-4)');shading flat;caxis(V);axis image;colorbar('vertical');
226			set(gca,'ydir','norm')
227			text(0,110,title);
228
229
230			%-----TEST-------------------------------------------------
231
232			sum=0;
233			sumT=0;
234			for i=Nx1:Nx2
235			for j=Ny1:Ny2
236			sum=sum+0.5meanw(i,j,Nz1)pv(i,j,Nz1)pv(i,j,Nz1)dx*dy;
237			sumT=sum-0.5meanw(i,j,Nz2)pv(i,j,Nz2)pv(i,j,Nz2)dx*dy;
238			end
239			end
240
241			for j=Ny1:Ny2
242			for k=Nz1:Nz2
243			sumT=sumT+0.5meanu(Nx2,j,k)pv(Nx2,j,k)pv(Nx2,j,k)dz*dy;
244			sumT=sumT-0.5meanu(Nx1,j,k)pv(Nx1,j,k)pv(Nx1,j,k)dz*dy;
245			end
246			end
247
248			for i=Nx1:Nx2
249			for k=Nz1:Nz2
250			sumT=sumT+0.5meanv(i,Ny2,k)pv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
251			sumT=sumT-0.5meanv(i,Ny1,k)pv(i,Ny1,k)pv(i,Ny1,k)dz*dx;
252			end
253			end
254			'VORTICITY GENERATION'
255			sum
256			sumT
257
258
259			sum=0;
260			for i=Nx1:Nx2
261			for j=Ny1:Ny2
262			sum=sum+sumD(i,j)dxdy;
263			end
264			end
265			'DISSIPATION BY EDDY-DIFFUSIVITY'
266			sum
267
268			sum=0;
269			for i=Nx1:Nx2
270			for j=Ny1:Ny2
271			sum=sum+sumM(i,j)dxdy;
272			end
273			end
274			'DISSIPATION BY EDDY-VISCOSITY'
275			sum