timour_matlab/mscripts/DIAGP.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
tfilename=(['T.' iterations((i),1:10) ]) ;
t=rdmds(tfilename,'b');

ufilename=(['U.' iterations((i),1:10) ]) ;
u=rdmds(ufilename,'b');

vfilename=(['V.' iterations((i),1:10) ]) ;
v=rdmds(vfilename,'b');

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');

vfilename=(['V.' iterations((i),1:10) ]) ;
v=rdmds(vfilename,'b');

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);
Nu=zeros(Nx,Ny,Nz);
Nv=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);

Nu(2:Nx-1,2:Ny-1,2:Nz-1)=-(meanu(3:Nx,2:Ny-1,2:Nz-1).*meanu(3:Nx,2:Ny-1,2:Nz-1) ...
        -meanu(1:Nx-2,2:Ny-1,2:Nz-1).*meanu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
        -(meanv(2:Nx-1,3:Ny,2:Nz-1).*meanu(2:Nx-1,3:Ny,2:Nz-1) ...
        -meanv(2:Nx-1,1:Ny-2,2:Nz-1).*meanu(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
        -(meanw(2:Nx-1,2:Ny-1,3:Nz).*meanu(2:Nx-1,2:Ny-1,3:Nz) ...
        -meanw(2:Nx-1,2:Ny-1,1:Nz-2).*meanu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);

Nv(2:Nx-1,2:Ny-1,2:Nz-1)=-(meanv(3:Nx,2:Ny-1,2:Nz-1).*meanu(3:Nx,2:Ny-1,2:Nz-1) ...
        -meanv(1:Nx-2,2:Ny-1,2:Nz-1).*meanu(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ...
        -(meanv(2:Nx-1,3:Ny,2:Nz-1).*meanv(2:Nx-1,3:Ny,2:Nz-1) ...
        -meanv(2:Nx-1,1:Ny-2,2:Nz-1).*meanv(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
        -(meanw(2:Nx-1,2:Ny-1,3:Nz).*meanv(2:Nx-1,2:Ny-1,3:Nz) ...
        -meanw(2:Nx-1,2:Ny-1,1:Nz-2).*meanv(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)=(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);

Mp=zeros(Nx,Ny,Nz);
Np=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).* ...
ff(3:Nx-2,3:Ny-2,3:Nz-2);

Np(3:Nx-2,3:Ny-2,3:Nz-2)=((Nv(4:Nx-1,3:Ny-2,3:Nz-2)-Nv(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx)-  ...
(Nu(3:Nx-2,4:Ny-1,3:Nz-2)-Nu(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) ...
        -(Nv(3:Nx-2,3:Ny-2,4:Nz-1)-Nv(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) ...
+(Nu(3:Nx-2,3:Ny-2,4:Nz-1)-Nu(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);

sumM=zeros(Nx,Ny);
sumN=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;
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;
sumN(3:Nx-2,3:Ny-2)=sumN(3:Nx-2,3:Ny-2)+Np(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=[-10 10];
        title='Vorticity input';
        imagesc(lat,long,sumW(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,sumD(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,sumN(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,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;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

sum=0;
for i=Nx1:Nx2
for j=Ny1:Ny2
sum=sum+sumN(i,j)*dx*dy;
end
end
'NONLINEARITY'
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			tfilename=(['T.' iterations((i),1:10) ]) ;
44			t=rdmds(tfilename,'b');
45
46			ufilename=(['U.' iterations((i),1:10) ]) ;
47			u=rdmds(ufilename,'b');
48
49			vfilename=(['V.' iterations((i),1:10) ]) ;
50			v=rdmds(vfilename,'b');
51
52			wfilename=(['W.' iterations((i),1:10) ]) ;
53			w=rdmds(wfilename,'b');
54
55			sumtheta=sumtheta+t;
56			sumu=sumu+u;
57			sumv=sumv+v;
58			sumw=sumw+w;
59			counter=counter+1;
60			end
61
62			meantheta=sumtheta/counter;
63			meanu=sumu/counter;
64			meanv=sumv/counter;
65			meanw=sumw/counter;
66
67			sumut=zeros(Nx,Ny,Nz);
68			sumvt=zeros(Nx,Ny,Nz);
69			sumwt=zeros(Nx,Ny,Nz);
70
71			sumuu=zeros(Nx,Ny,Nz);
72			sumvu=zeros(Nx,Ny,Nz);
73			sumwu=zeros(Nx,Ny,Nz);
74			sumvv=zeros(Nx,Ny,Nz);
75			sumwv=zeros(Nx,Ny,Nz);
76
77			meanu=smooth3(meanu);
78			meanv=smooth3(meanv);
79			meanw=smooth3(meanw);
80			meantheta=smooth3(meantheta);
81
82			for i=nitstart:nitend
83			i
84			tfilename=(['T.' iterations((i),1:10) ]) ;
85			t=rdmds(tfilename,'b');
86
87			ufilename=(['U.' iterations((i),1:10) ]) ;
88			u=rdmds(ufilename,'b');
89
90			vfilename=(['V.' iterations((i),1:10) ]) ;
91			v=rdmds(vfilename,'b');
92
93			wfilename=(['W.' iterations((i),1:10) ]) ;
94			w=rdmds(wfilename,'b');
95
96			t=t-meantheta;
97			u=u-meanu;
98			v=v-meanv;
99			w=w-meanw;
100
101			sumut=sumut+u.*t;
102			sumvt=sumvt+v.*t;
103			sumwt=sumwt+w.*t;
104
105			sumuu=sumuu+u.*u;
106			sumvu=sumvu+v.*u;
107			sumwu=sumwu+w.*u;
108			sumvv=sumvv+v.*v;
109			sumwv=sumwv+w.*v;
110			end
111
112			sumut=sumut/counter;
113			sumvt=sumvt/counter;
114			sumwt=sumwt/counter;
115
116			sumuu=sumuu/counter;
117			sumvu=sumvu/counter;
118			sumwu=sumwu/counter;
119			sumvv=sumvv/counter;
120			sumwv=sumwv/counter;
121
122			Mu=zeros(Nx,Ny,Nz);
123			Mv=zeros(Nx,Ny,Nz);
124			Nu=zeros(Nx,Ny,Nz);
125			Nv=zeros(Nx,Ny,Nz);
126			D=zeros(Nx,Ny,Nz);
127
128			dx=dm;
129			dy=dm;
130			dz=-dz;
131
132			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) ...
133			-(sumvu(2:Nx-1,3:Ny,2:Nz-1)-sumvu(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
134			-(sumwu(2:Nx-1,2:Ny-1,3:Nz)-sumwu(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
135
136			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) ...
137			-(sumvv(2:Nx-1,3:Ny,2:Nz-1)-sumvv(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
138			-(sumwv(2:Nx-1,2:Ny-1,3:Nz)-sumwv(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
139
140			Nu(2:Nx-1,2:Ny-1,2:Nz-1)=-(meanu(3:Nx,2:Ny-1,2:Nz-1).*meanu(3:Nx,2:Ny-1,2:Nz-1) ...
141			-meanu(1:Nx-2,2:Ny-1,2:Nz-1).meanu(1:Nx-2,2:Ny-1,2:Nz-1))/(2dx) ...
142			-(meanv(2:Nx-1,3:Ny,2:Nz-1).*meanu(2:Nx-1,3:Ny,2:Nz-1) ...
143			-meanv(2:Nx-1,1:Ny-2,2:Nz-1).meanu(2:Nx-1,1:Ny-2,2:Nz-1))/(2dy) ...
144			-(meanw(2:Nx-1,2:Ny-1,3:Nz).*meanu(2:Nx-1,2:Ny-1,3:Nz) ...
145			-meanw(2:Nx-1,2:Ny-1,1:Nz-2).meanu(2:Nx-1,2:Ny-1,1:Nz-2))/(2dz);
146
147			Nv(2:Nx-1,2:Ny-1,2:Nz-1)=-(meanv(3:Nx,2:Ny-1,2:Nz-1).*meanu(3:Nx,2:Ny-1,2:Nz-1) ...
148			-meanv(1:Nx-2,2:Ny-1,2:Nz-1).meanu(1:Nx-2,2:Ny-1,2:Nz-1))/(2dx) ...
149			-(meanv(2:Nx-1,3:Ny,2:Nz-1).*meanv(2:Nx-1,3:Ny,2:Nz-1) ...
150			-meanv(2:Nx-1,1:Ny-2,2:Nz-1).meanv(2:Nx-1,1:Ny-2,2:Nz-1))/(2dy) ...
151			-(meanw(2:Nx-1,2:Ny-1,3:Nz).*meanv(2:Nx-1,2:Ny-1,3:Nz) ...
152			-meanw(2:Nx-1,2:Ny-1,1:Nz-2).meanv(2:Nx-1,2:Ny-1,1:Nz-2))/(2dz);
153
154			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) ...
155			-(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ...
156			-(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz);
157
158			pv=zeros(Nx,Ny,Nz);
159			ff=zeros(Nx,Ny,Nz);
160
161
162			for j=1:Ny
163			ff(:,j,:)=f+(j-1)dybeta;
164			end
165
166			pv(2:Nx-1,2:Ny-1,2:Nz-1)=(ff(2:Nx-1,2:Ny-1,2:Nz-1)) ...
167			.(meantheta(2:Nx-1,2:Ny-1,3:Nz)-meantheta(2:Nx-1,2:Ny-1,1:Nz-2))/(2dz);
168
169			Mp=zeros(Nx,Ny,Nz);
170			Np=zeros(Nx,Ny,Nz);
171			Dp=zeros(Nx,Ny,Nz);
172
173			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)- ...
174			(Mu(3:Nx-2,4:Ny-1,3:Nz-2)-Mu(3:Nx-2,2:Ny-3,3:Nz-2))/(2dy)). ...
175			(meantheta(3:Nx-2,3:Ny-2,4:Nz-1)-meantheta(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz) ...
176			-(Mv(3:Nx-2,3:Ny-2,4:Nz-1)-Mv(3:Nx-2,3:Ny-2,2:Nz-3))/(2dz). ...
177			(meantheta(4:Nx-1,3:Ny-2,3:Nz-2)-meantheta(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx) ...
178			+(Mu(3:Nx-2,3:Ny-2,4:Nz-1)-Mu(3:Nx-2,3:Ny-2,2:Nz-3))/(2dz). ...
179			(meantheta(3:Nx-2,4:Ny-1,3:Nz-2)-meantheta(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy);
180
181			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). ...
182			ff(3:Nx-2,3:Ny-2,3:Nz-2);
183
184			Np(3:Nx-2,3:Ny-2,3:Nz-2)=((Nv(4:Nx-1,3:Ny-2,3:Nz-2)-Nv(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx)- ...
185			(Nu(3:Nx-2,4:Ny-1,3:Nz-2)-Nu(3:Nx-2,2:Ny-3,3:Nz-2))/(2dy)). ...
186			(meantheta(3:Nx-2,3:Ny-2,4:Nz-1)-meantheta(3:Nx-2,3:Ny-2,2:Nz-3))/(2*dz) ...
187			-(Nv(3:Nx-2,3:Ny-2,4:Nz-1)-Nv(3:Nx-2,3:Ny-2,2:Nz-3))/(2dz). ...
188			(meantheta(4:Nx-1,3:Ny-2,3:Nz-2)-meantheta(2:Nx-3,3:Ny-2,3:Nz-2))/(2*dx) ...
189			+(Nu(3:Nx-2,3:Ny-2,4:Nz-1)-Nu(3:Nx-2,3:Ny-2,2:Nz-3))/(2dz). ...
190			(meantheta(3:Nx-2,4:Ny-1,3:Nz-2)-meantheta(3:Nx-2,2:Ny-3,3:Nz-2))/(2*dy);
191
192			sumM=zeros(Nx,Ny);
193			sumN=zeros(Nx,Ny);
194			sumD=zeros(Nx,Ny);
195			sumW=zeros(Nx,Ny);
196
197			dz=-dz;
198
199			% DEFINE BOX OF INTEGRATION
200			Nx1=3;
201			Nx2=Nx-2;
202			Ny1=3;
203			Ny2=Ny-2;
204			Nz1=5;
205			Nz2=Nz-4;
206
207
208			for k=Nz1:Nz2
209			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;
210			sumN(3:Nx-2,3:Ny-2)=sumN(3:Nx-2,3:Ny-2)+Np(3:Nx-2,3:Ny-2,k).pv(3:Nx-2,3:Ny-2,k)dz;
211			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;
212			end
213			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);
214
215			%contourf(sumM',20);colorbar;
216			%figure
217			%contourf(sumD',20);colorbar;
218
219			V=[-10 10];
220			title='Vorticity input';
221			imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
222			set(gca,'ydir','norm')
223			text(0,110,title);
224			figure
225			title='Buoyancy diffusion integral';
226			imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
227			set(gca,'ydir','norm')
228			text(0,110,title);
229			figure
230			title='Non-linear terms integral';
231			imagesc(lat,long,sumN(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
232			set(gca,'ydir','norm')
233			text(0,110,title);
234			figure
235			title='Momentum diffusion integral';
236			imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;axis image;colorbar('vertical');
237			set(gca,'ydir','norm')
238			text(0,110,title);
239
240
241			%-----TEST-------------------------------------------------
242
243			sum=0;
244			sumT=0;
245			for i=Nx1:Nx2
246			for j=Ny1:Ny2
247			sum=sum+0.5meanw(i,j,Nz1)pv(i,j,Nz1)pv(i,j,Nz1)dx*dy;
248			sumT=sum-0.5meanw(i,j,Nz2)pv(i,j,Nz2)pv(i,j,Nz2)dx*dy;
249			end
250			end
251
252			for j=Ny1:Ny2
253			for k=Nz1:Nz2
254			sumT=sumT+0.5meanu(Nx2,j,k)pv(Nx2,j,k)pv(Nx2,j,k)dz*dy;
255			sumT=sumT-0.5meanu(Nx1,j,k)pv(Nx1,j,k)pv(Nx1,j,k)dz*dy;
256			end
257			end
258
259			for i=Nx1:Nx2
260			for k=Nz1:Nz2
261			sumT=sumT+0.5meanv(i,Ny2,k)pv(i,Ny2,k)pv(i,Ny2,k)dz*dx;
262			sumT=sumT-0.5meanv(i,Ny1,k)pv(i,Ny1,k)pv(i,Ny1,k)dz*dx;
263			end
264			end
265			'VORTICITY GENERATION'
266			sum
267			sumT
268
269
270			sum=0;
271			for i=Nx1:Nx2
272			for j=Ny1:Ny2
273			sum=sum+sumD(i,j)dxdy;
274			end
275			end
276			'DISSIPATION BY EDDY-DIFFUSIVITY'
277			sum
278
279			sum=0;
280			for i=Nx1:Nx2
281			for j=Ny1:Ny2
282			sum=sum+sumM(i,j)dxdy;
283			end
284			end
285			'DISSIPATION BY EDDY-VISCOSITY'
286			sum
287
288			sum=0;
289			for i=Nx1:Nx2
290			for j=Ny1:Ny2
291			sum=sum+sumN(i,j)dxdy;
292			end
293			end
294			'NONLINEARITY'
295			sum