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