clear all; 
close all;
clc;

 
    
 
      
load constsandscalingLarsen.mat

load LarsenHWM2002.mat
alphawind1=307;
rejimissaxeli=strcat('Uwind-LarsenHWM2002','-Da','-z0=120-Heff=8');
disp(kmin)
disp(kmax)
%rejimissaxeli='U0=200-Da-alpha315-z0=120-Heff=8';
   alphawind=alphawind1*pi/180;     

for j=1:Nx    
for k=1:Nz
  
recombcoef(j,k)=0;
N(j,k)=exp(-(((zmnorm-z(k))/sigmaHi)^2)-0*(((ymnorm-y(j))/sigmaHi)^2));
%exp(-(((zmnorm-z(k))/sigmaHi)^2)-0*(((ymnorm-y(j))/sigmaHi)^2));
Ngauss(j,k)=exp(-(((zmnorm-z(k))/sigmaHi)^2)-0*(((ymnorm-y(j))/sigmaHi)^2));
%exp(-(((zmnorm-z(k))/sigmaHi)^2)-0*(((ymnorm-y(j))/sigmaHi)^2));
     
%kappaspnormmsise(j,k)=kappaspnormmsise(jmin,500);
  %VnX(k)=Uwind(k)*cos(alphawind);
   %VnY(k)=Uwind(k)*sin(alphawind);
end
end










       

%%%%%%%%%%%%%
%%%%%%%%%%%%%

Nfirstlab=1;

Nlastlab=2;
deltak=round((kmax-5-kmin)/(Nlastlab-Nfirstlab+1));
            Jmin(Nfirstlab)=jmin;
            Jmax(Nfirstlab)=jmax;
            Kmin(Nfirstlab)=kmin;           
            Jmin(Nlastlab)=jmin;
            Jmax(Nlastlab)=jmax;
           
            Kmax(Nlastlab)=kmax;

      for i=Nfirstlab:Nlastlab-1
                  
            Jmax(i)=Jmax(Nfirstlab);
            Kmax(i)=Kmin(Nfirstlab)+deltak*i ;            
            
      end

   for i=Nfirstlab+1:Nlastlab
                   
            Jmin(i)=Jmin(Nfirstlab);
            Kmin(i)=Kmax(i-1)-1; 
            
            
    end

%parpool(Nlastlab)

spmd(2)


for t=3:12000




for j=Jmin(labindex):Jmax(labindex)+1
   for h=Kmin(labindex):Kmax(labindex)+1

       


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%ion drift velocitties  pressure  not included    
      
       ViXnop(j,h)=(VnX(h)*(kappaspnormmsise(j,h)*kappaspnormmsise(j,h)+COSI*COSI)- VnY(h)*kappaspnormmsise(j,h)*SINI-VnZ(h)*SINI*COSI)/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
       %ViXnop=(-VnZ(j,h)*SINI*COSI-VnX(j,h)*SINI*kappaspnormmsise(j,h)+VnX(j,h)*(kappaspnormmsise(j,h)*kappaspnormmsise(j,h)+COSI*COSI))/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
       %ViXnop=(VnX(j,h)*SINI+VnX(j,h)*kappaspnormmsise(j,h)+VnZ(j,h)*COSI)*kappaspnormmsise(j,h)/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
           
       ViZnop(j,h)=(-VnX(h)*SINI*COSI-VnY(h)*kappaspnormmsise(j,h)*COSI+VnZ(h)*(kappaspnormmsise(j,h)*kappaspnormmsise(j,h)+SINI*SINI))/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
      
      

   end
end
 
  
for j=Jmin(labindex):Jmax(labindex)
   for h=Kmin(labindex):Kmax(labindex)+1   
      ViXj2(j,h)=0.5*(ViXnop(j,h)+ViXnop(j+1,h));
      
      if ViXj2(j,h)>0
          tetaXj2(j,h)=1;
          rXj2(j,h)=(N(j,h)-N(j-1,h))/(N(j+1,h)-N(j,h));
                else
          rXj2(j,h)=(N(j+2,h)-N(j+1,h))/(N(j+1,h)-N(j,h));
          tetaXj2(j,h)=-1;
      end
      psiXj2(j,h)=0;
   end
end
for j=Jmin(labindex):Jmax(labindex)+1
   for h=Kmin(labindex):Kmax(labindex) 
       ViZh2(j,h)=0.5*(ViZnop(j,h)+ViZnop(j,h+1));
        if ViZh2(j,h)>0
          tetaZh2(j,h)=1;
          rZh2(j,h)=(N(j,h)-N(j,h-1))/(N(j,h+1)-N(j,h));
        else
           rZh2(j,h)=(N(j,h+2)-N(j,h+1))/(N(j,h+1)-N(j,h));
          tetaZh2(j,h)=-1;
        end
      psiZh2(j,h)=0;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% fluxes



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% psiflux(y/z)
psiflux=0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% flows(y/z)
for j=Jmin(labindex):Jmax(labindex)-1   
   for k=Kmin(labindex):Kmax(labindex) 
      
       FX2t2(j,k)=0.5*ViXj2(j,k)*((1+tetaXj2(j,k))*N(j,k)+(1-tetaXj2(j,k))*N(j+1,k))+0.5*abs(ViXj2(j,k))*(1-abs(ViXj2(j,k)*0.5*dt/dx))*psiXj2(j,k)*(N(j+1,k)-N(j,k));
       
       %{
       if ViXj2(j,k)>0
         FX2t2(j,k)=(N(j,k)+psiflux*(N(j+1,k)-N(j,k)))*ViXj2(j,k);
      else 
         FX2t2(j,k)=(N(j+1,k)+psiflux*(N(j,k)-N(j+1,k)))*ViXj2(j,k);
      end
       %}
   end
end

for j=Jmin(labindex):Jmax(labindex)   
   for k=Kmin(labindex):Kmax(labindex)-1 
      FZ2t2(j,k)=0.5*ViZh2(j,k)*((1+tetaZh2(j,k))*N(j,k)+(1-tetaZh2(j,k))*N(j,k+1))+0.5*abs(ViZh2(j,k))*(1-abs(ViZh2(j,k)*0.5*dt/(z(k)-z(k-1))))*psiZh2(j,k)*(N(j,k+1)-N(j,k));

       %{
      if  ViZh2(j,k)>0
         FZ2t2(j,k)=(N(j,k)+psiflux*(N(j,k+1)-N(j,k)))*ViZh2(j,k);
      else
         FZ2t2(j,k)=(N(j,k+1)+psiflux*(N(j,k)-N(j,k+1)))*ViZh2(j,k); 
      end
       %}
   end
end 
 %NdiffSTS1winiswinanegativetotal=0;
 %NdiffSTS1winiswinapositivetotal=0;
for j=Jmin(labindex)+1:Jmax(labindex)-1
for k=Kmin(labindex)+1:Kmax(labindex)-1

NdiffSTS1winiswina(j,k)=N(j,k)-((FX2t2(j,k)-FX2t2(j-1,k))*(0.5*dt/dx)+ (FZ2t2(j,k)-FZ2t2(j,k-1))*((0.5*dt)/(z(k)-z(k-1))));





if NdiffSTS1winiswina(j,k)<0
    

   NdiffSTS1winiswina(j,k)=0;  %        
         
end
 
end

end


%%%Groundfloor

if labindex==Nfirstlab
%{
NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
   
%}
NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
   

if NdiffSTS1winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))<0
        NdiffSTS1winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=0;    
end

end
%%%%rooof
if labindex==Nlastlab
 %{   
NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1winiswina(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1winiswina(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
%}

NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1winiswina(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1winiswina(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');



if NdiffSTS1winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)<0
         NdiffSTS1winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=0;    
end

end
%%%%walls
NdiffSTS1winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');


if NdiffSTS1winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)<0
         NdiffSTS1winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end

if NdiffSTS1winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)<0
 NdiffSTS1winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end









if labindex<Nlastlab


 for j=Jmin(labindex)-3:Jmax(labindex)+3
 for k=1:4 %Kmin(labindex)+1:Kmin(labindex)+4
a1sentup(j,k)=NdiffSTS1winiswina(j,Kmax(labindex)-k);
end
end


labSend(a1sentup,labindex+1,str2num(strcat(int2str(labindex),int2str(1),int2str(t))))
 end




   if labindex>Nfirstlab

for j=Jmin(labindex)-3:Jmax(labindex)+3
for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
a1sentdown(j,k)=NdiffSTS1winiswina(j,Kmin(labindex)+k);
end
end


labSend(a1sentdown,labindex-1,str2num(strcat(int2str(labindex),int2str(1),int2str(t))))
    end


 
  
  if labindex<Nlastlab

a1receivedfromupstairs=labReceive(labindex+1,str2num(strcat(int2str(labindex+1),int2str(1),int2str(t))));

for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS1winiswina(j,Kmax(labindex)-1+k)=a1receivedfromupstairs(j,k);

          end
end

  end

    if labindex>Nfirstlab

a1receivedfromdownstairs=labReceive(labindex-1,str2num(strcat(int2str(labindex-1),int2str(1),int2str(t))));

for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS1winiswina(j,Kmin(labindex)-k+1)=a1receivedfromdownstairs(j,k);

          end
end
    end





for j=Jmin(labindex)-3:Jmax(labindex)+3
   for k=Kmin(labindex)-3:Kmax(labindex)+3
 NdiffSTS1wina(j,k)=NdiffSTS1winiswina(j,k);
   end
end


 

for i=2:NSTS
   
for j=Jmin(labindex)+1:Jmax(labindex)-1
   for k=Kmin(labindex)+1:Kmax(labindex)-1
       NdiffSTS1(j,k)=((1-betta(j,k)*0.5*dtdiff(i)-gamma(j,k)*0.5*dtdiff(i))/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*NdiffSTS1winiswina(j,k)+...
                           +(bettax(j,k)*0.5*dtdiff(i)/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*(NdiffSTS1wina(j+1,k)+NdiffSTS1wina(j-1,k))+...
                           +(bettaz(j,k)*0.5*dtdiff(i)/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*(NdiffSTS1wina(j,k+1)+NdiffSTS1wina(j,k-1))+...
                           +(gamma(j,k)*0.5*dtdiff(i)/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*(NdiffSTS1wina(j+1,k+1)+NdiffSTS1wina(j-1,k-1));                        
         
  

if NdiffSTS1winiswina(j,k)<0
       NdiffSTS1(j,k)=0;  
         
         
end
     
   end
end
%%%Groundfloor

if labindex==Nfirstlab
    %{
NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
%}
NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');


if NdiffSTS1(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))<0
        NdiffSTS1(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=0;    
end

end
%%%%rooof

if labindex==Nlastlab
    %{
NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

%}
NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS1(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS1(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if NdiffSTS1(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)<0
         NdiffSTS1(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=0;    
end
end
%%%%walls
NdiffSTS1(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS1(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS1(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS1(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if NdiffSTS1(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)<0
         NdiffSTS1(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end

if NdiffSTS1(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)<0
 NdiffSTS1(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end











if labindex<Nlastlab


 for j=Jmin(labindex)-3:Jmax(labindex)+3
 for k=1:4 %Kmin(labindex)+1:Kmin(labindex)+4
a3sentup(j,k)=NdiffSTS1(j,Kmax(labindex)-k);
end
end


labSend(a3sentup,labindex+1,str2num(strcat(int2str(labindex),int2str(3),int2str(t))))
 end




   if labindex>Nfirstlab

for j=Jmin(labindex)-3:Jmax(labindex)+3
for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
a3sentdown(j,k)=NdiffSTS1(j,Kmin(labindex)+k);
end
end


labSend(a3sentdown,labindex-1,str2num(strcat(int2str(labindex),int2str(3),int2str(t))))
    end


 
  
  if labindex<Nlastlab

a3receivedfromupstairs=labReceive(labindex+1,str2num(strcat(int2str(labindex+1),int2str(3),int2str(t))));

for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS1(j,Kmax(labindex)-1+k)=a3receivedfromupstairs(j,k);

          end
end

  end

    if labindex>Nfirstlab

a3receivedfromdownstairs=labReceive(labindex-1,str2num(strcat(int2str(labindex-1),int2str(3),int2str(t))));


for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS1(j,Kmin(labindex)-k+1)=a3receivedfromdownstairs(j,k);

          end
end
    end


    
    for j=Jmin(labindex)-3:Jmax(labindex)+3
   for k=Kmin(labindex)-3:Kmax(labindex)+3
    NdiffSTS1winiswina(j,k)=NdiffSTS1wina(j,k);
    NdiffSTS1wina(j,k)=NdiffSTS1(j,k);
   end
    end
end






for j=Jmin(labindex)+1:Jmax(labindex)-1
   for k=Kmin(labindex)+1:Kmax(labindex)-1
                 NdiffSTS2winiswina(j,k)=NdiffSTS1(j,k)+dt*NdiffSTS1(j,k)*(dzDznorm(j,k)*Hnorm/(2*H)-recombcoef(k))-(dzDxnorm(j,k)*(NdiffSTS1(j+1,k)-NdiffSTS1(j-1,k))*(dt/dx)-(dzDznorm(j,k)+Dznorm(j,k)*Hnorm/(2*H))*(NdiffSTS1(j,k+1)-NdiffSTS1(j,k-1))*(dt/(z(k+1)-z(k-1))));
   
   
 

if NdiffSTS2winiswina(j,k)<0

       NdiffSTS2winiswina(j,k)=0;  
         
         
end

    end

end


 %%%Groundfloor

if labindex==Nfirstlab
 %{   
NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
 %} 
NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');


if NdiffSTS2winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))<0
        NdiffSTS2winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=0;    
end

end
%%%%rooof

if labindex==Nlastlab
    
%{    
NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2winiswina(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2winiswina(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
%}

NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2winiswina(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2winiswina(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2winiswina(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');


if NdiffSTS2winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)<0
         NdiffSTS2winiswina(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=0;    
end

end
%%%%walls
NdiffSTS2winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2winiswina(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2winiswina(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if NdiffSTS2winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)<0
         NdiffSTS2winiswina(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end

if NdiffSTS2winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)<0
 NdiffSTS2winiswina(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end



 







if labindex<Nlastlab


 for j=Jmin(labindex)-3:Jmax(labindex)+3
 for k=1:4 %Kmin(labindex)+1:Kmin(labindex)+4
a4sentup(j,k)=NdiffSTS2winiswina(j,Kmax(labindex)-k);
end
end


labSend(a4sentup,labindex+1,str2num(strcat(int2str(labindex),int2str(4),int2str(t))))
 end




   if labindex>Nfirstlab

for j=Jmin(labindex)-3:Jmax(labindex)+3
for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
a4sentdown(j,k)=NdiffSTS2winiswina(j,Kmin(labindex)+k);
end
end


labSend(a4sentdown,labindex-1,str2num(strcat(int2str(labindex),int2str(4),int2str(t))))
    end


 
  
  if labindex<Nlastlab

a4receivedfromupstairs=labReceive(labindex+1,str2num(strcat(int2str(labindex+1),int2str(4),int2str(t))));

for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS2winiswina(j,Kmax(labindex)-1+k)=a4receivedfromupstairs(j,k);

          end
end

  end

    if labindex>Nfirstlab

a4receivedfromdownstairs=labReceive(labindex-1,str2num(strcat(int2str(labindex-1),int2str(4),int2str(t))));


for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS2winiswina(j,Kmin(labindex)-k+1)=a4receivedfromdownstairs(j,k);

          end
end
    end


for j=Jmin(labindex)-3:Jmax(labindex)+3
   for k=Kmin(labindex)-3:Kmax(labindex)+3
NdiffSTS2wina(j,k)=NdiffSTS2winiswina(j,k);
   end
end



      
for i=2:NSTS
  
for j=Jmin(labindex)+1:Jmax(labindex)-1
   for k=Kmin(labindex)+1:Kmax(labindex)-1
       NdiffSTS2(j,k)=((1-betta(j,k)*0.5*dtdiff(i)-gamma(j,k)*0.5*dtdiff(i))/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*NdiffSTS2winiswina(j,k)+...
                           +(bettax(j,k)*0.5*dtdiff(i)/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*(NdiffSTS2wina(j+1,k)+NdiffSTS2wina(j-1,k))+...
                          +(bettaz(j,k)*0.5*dtdiff(i)/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*(NdiffSTS2wina(j,k+1)+NdiffSTS2wina(j,k-1))+...
                           +(gamma(j,k)*0.5*dtdiff(i)/(1+betta(j,k)*0.5*dtdiff(i)+gamma(j,k)*0.5*dtdiff(i)))*(NdiffSTS2wina(j+1,k+1)+NdiffSTS2wina(j-1,k-1));                        
                 
       

if NdiffSTS2(j,k)<0
        NdiffSTS2(j,k)=0;  
         
         
end

   end
end
%%%Groundfloor

if labindex==Nfirstlab
 %{   
NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
 %}
NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if NdiffSTS2(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))<0
        NdiffSTS2(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=0;    
end

end
%%%%rooof

if labindex==Nlastlab
 %{   
NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
%}
NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),NdiffSTS2(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
NdiffSTS2(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if NdiffSTS2(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)<0
         NdiffSTS2(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=0;    
end

end
%%%%walls
NdiffSTS2(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),NdiffSTS2(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
NdiffSTS2(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),NdiffSTS2(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if NdiffSTS2(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)<0
         NdiffSTS2(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end

if NdiffSTS2(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)<0
 NdiffSTS2(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end











if labindex<Nlastlab


 for j=Jmin(labindex)-3:Jmax(labindex)+3
 for k=1:4 %Kmin(labindex)+1:Kmin(labindex)+4
a6sentup(j,k)=NdiffSTS2(j,Kmax(labindex)-k);
end
end


labSend(a6sentup,labindex+1,str2num(strcat(int2str(labindex),int2str(6),int2str(t))))
 end




   if labindex>Nfirstlab

for j=Jmin(labindex)-3:Jmax(labindex)+3
for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
a6sentdown(j,k)=NdiffSTS2(j,Kmin(labindex)+k);
end
end


labSend(a6sentdown,labindex-1,str2num(strcat(int2str(labindex),int2str(6),int2str(t))))
    end


 
  
  if labindex<Nlastlab

a6receivedfromupstairs=labReceive(labindex+1,str2num(strcat(int2str(labindex+1),int2str(6),int2str(t))));

for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS2(j,Kmax(labindex)-1+k)=a6receivedfromupstairs(j,k);

          end
end

  end

    if labindex>Nfirstlab

a6receivedfromdownstairs=labReceive(labindex-1,str2num(strcat(int2str(labindex-1),int2str(6),int2str(t))));


for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
NdiffSTS2(j,Kmin(labindex)-k+1)=a6receivedfromdownstairs(j,k);

          end
end
    end
  
 for j=Jmin(labindex)-3:Jmax(labindex)+3
   for k=Kmin(labindex)-3:Kmax(labindex)+3
    NdiffSTS2winiswina(j,k)=NdiffSTS2wina(j,k);
    NdiffSTS2wina(j,k)=NdiffSTS2(j,k);
   end
    end

end







for j=Jmin(labindex):Jmax(labindex)+1
   for h=Kmin(labindex):Kmax(labindex)+1

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%ion drift velocitties  pressure  not included    
   
       ViXnop(j,h)=(VnX(h)*(kappaspnormmsise(j,h)*kappaspnormmsise(j,h)+COSI*COSI)-VnY(h)*kappaspnormmsise(j,h)*SINI-VnZ(h)*SINI*COSI)/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
       %ViXnop=(-VnZ(j,h)*SINI*COSI-VnX(j,h)*SINI*kappaspnormmsise(j,h)+VnX(j,h)*(kappaspnormmsise(j,h)*kappaspnormmsise(j,h)+COSI*COSI))/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
       %ViXnop=(VnX(j,h)*SINI+VnX(j,h)*kappaspnormmsise(j,h)+VnZ(j,h)*COSI)*kappaspnormmsise(j,h)/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
           
       ViZnop(j,h)=(-VnX(h)*SINI*COSI-VnY(h)*kappaspnormmsise(j,h)*COSI+VnZ(h)*(kappaspnormmsise(j,h)*kappaspnormmsise(j,h)+SINI*SINI))/(1+kappaspnormmsise(j,h)* kappaspnormmsise(j,h));
  
 

   end
end
for j=Jmin(labindex):Jmax(labindex)
   for h=Kmin(labindex):Kmax(labindex)+1   
      ViXj2(j,h)=0.5*(ViXnop(j,h)+ViXnop(j+1,h));
      
      if ViXj2(j,h)>0
          tetaXj2(j,h)=1;
          rXj2(j,h)=(NdiffSTS2(j,h)-NdiffSTS2(j-1,h))/(NdiffSTS2(j+1,h)-NdiffSTS2(j,h));
                else
          rXj2(j,h)=(NdiffSTS2(j+2,h)-NdiffSTS2(j+1,h))/(NdiffSTS2(j+1,h)-NdiffSTS2(j,h));
          tetaXj2(j,h)=-1;
      end
      psiXj2(j,h)=0;
   end
end
for j=Jmin(labindex):Jmax(labindex)+1
   for h=Kmin(labindex):Kmax(labindex) 
       ViZh2(j,h)=0.5*(ViZnop(j,h)+ViZnop(j,h+1));
        if ViZh2(j,h)>0
          tetaZh2(j,h)=1;
          rZh2(j,h)=(NdiffSTS2(j,h)-NdiffSTS2(j,h-1))/(NdiffSTS2(j,h+1)-NdiffSTS2(j,h));
        else
           rZh2(j,h)=(NdiffSTS2(j,h+2)-NdiffSTS2(j,h+1))/(NdiffSTS2(j,h+1)-NdiffSTS2(j,h));
          tetaZh2(j,h)=-1;
        end
      psiZh2(j,h)=0;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% fluxes


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% psiflux(y/z)
psiflux=0.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% flows(y/z)
for j=Jmin(labindex):Jmax(labindex)-1   
   for k=Kmin(labindex):Kmax(labindex) 
       FX2t2(j,k)=0.5*ViXj2(j,k)*((1+tetaXj2(j,k))*NdiffSTS2(j,k)+(1-tetaXj2(j,k))*NdiffSTS2(j+1,k))+0.5*abs(ViXj2(j,k))*(1-abs(ViXj2(j,k)*0.5*dt/dx))*psiXj2(j,k)*(NdiffSTS2(j+1,k)-NdiffSTS2(j,k));
     
       %{
      if ViXj2(j,k)>0
         FX2t(j,k)=(NdiffSTS2(j,k)+psiflux*(NdiffSTS2(j+1,k)-NdiffSTS2(j,k)))*ViXj2(j,k);
      else 
         FX2t(j,k)=(NdiffSTS2(j+1,k)+psiflux*(NdiffSTS2(j,k)-NdiffSTS2(j+1,k)))*ViXj2(j,k);
      end
       %}
   end
end

for j=Jmin(labindex):Jmax(labindex)   
   for k=Kmin(labindex):Kmax(labindex)-1 
       FZ2t2(j,k)=0.5*ViZh2(j,k)*((1+tetaZh2(j,k))*NdiffSTS2(j,k)+(1-tetaZh2(j,k))*NdiffSTS2(j,k+1))+0.5*abs(ViZh2(j,k))*(1-abs(ViZh2(j,k)*0.5*dt/(z(k)-z(k-1))))*psiZh2(j,k)*(NdiffSTS2(j,k+1)-NdiffSTS2(j,k));

       %{
      if  ViZh2(j,k)>0
         FZ2t(j,k)=(NdiffSTS2(j,k)+psiflux*(NdiffSTS2(j,k+1)-NdiffSTS2(j,k)))*ViZh2(j,k);
      else
         FZ2t(j,k)=(NdiffSTS2(j,k+1)+psiflux*(NdiffSTS2(j,k)-NdiffSTS2(j,k+1)))*ViZh2(j,k); 
      end
       %}
   end
end 

%Nnegativetotal=0;
%Npositivetotal=0;
for j=Jmin(labindex)+1:Jmax(labindex)-1
for k=Kmin(labindex)+1:Kmax(labindex)-1

N(j,k)=NdiffSTS2(j,k)-((FX2t2(j,k)-FX2t2(j-1,k))*(0.5*dt/dx)+(FZ2t2(j,k)-FZ2t2(j,k-1))*((0.5*dt)/(z(k)-z(k-1))));




if N(j,k)<0
         N(j,k)=0;  
         
         

end

%ratio(j,k)=NdiffSTS2(j,k)/((FX2t(j,k)-FX2t(j-1,k))*(1/dx)+(FZ2t(j,k)-FZ2t(j,k-1))*((1)/(z(k)-z(k-1))));

end
end
   

%if labindex==Nfirstlab
%%%Groundfloor

if labindex==Nfirstlab
 %{   
N(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),N(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
N(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),N(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
N(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),N(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
 %}

N(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)+1:Kmin(labindex)+5),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),N(Jmin(labindex)+1:Jmax(labindex)-1,Kmin(labindex)+1:Kmin(labindex)+5),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
N(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),N(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
N(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex));%interp2(Kdom(Kmin(labindex)-3:Kmin(labindex)),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),N(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)-3:Kmin(labindex)),Ktotal(Kmin(labindex)-3:Kmin(labindex)),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if N(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))<0
        N(Jmin(labindex)-3:Jmax(labindex)+3,Kmin(labindex)-3:Kmin(labindex))=0;    
end

end
%%%%rooof

if labindex==Nlastlab
   %{ 
N(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),N(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
N(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),N(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
N(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),N(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');
%}
N(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex)-5:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmax(labindex)-1),N(Jmin(labindex)+1:Jmax(labindex)-1,Kmax(labindex)-5:Kmax(labindex)-1),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)+1:Jmax(labindex)-1,1),'spline');
N(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmin(labindex)-3:Jmin(labindex),Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),N(Jmin(labindex)+1:Jmin(labindex)+5,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
N(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=Ngauss(Jmax(labindex):Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3);%interp2(Kdom(Kmax(labindex):Kmax(labindex)+3),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),N(Jmax(labindex)-5:Jmax(labindex)-1,Kmax(labindex):Kmax(labindex)+3),Ktotal(Kmax(labindex):Kmax(labindex)+3),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if N(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)<0
         N(Jmin(labindex)-3:Jmax(labindex)+3,Kmax(labindex):Kmax(labindex)+3)=0;    
end

end
%%%%walls
N(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmin(labindex)+1:Jmin(labindex)+5),N(Jmin(labindex)+1:Jmin(labindex)+5,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmin(labindex)-3:Jmin(labindex),1),'spline');
N(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=interp2(Kdom(Kmin(labindex)+1:Kmax(labindex)-1),Jdom(Jmax(labindex)-5:Jmax(labindex)-1),N(Jmax(labindex)-5:Jmax(labindex)-1,Kmin(labindex)+1:Kmax(labindex)-1),Ktotal(Kmin(labindex)+1:Kmax(labindex)-1),Jtotal(Jmax(labindex):Jmax(labindex)+3,1),'spline');

if N(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)<0
         N(Jmin(labindex)-3:Jmin(labindex),Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end

if N(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)<0
 N(Jmax(labindex):Jmax(labindex)+3,Kmin(labindex)+1:Kmax(labindex)-1)=0;    
end









if labindex<Nlastlab


 for j=Jmin(labindex)-3:Jmax(labindex)+3
 for k=1:4 %Kmin(labindex)+1:Kmin(labindex)+4
a7sentup(j,k)=N(j,Kmax(labindex)-k);
end
end


labSend(a7sentup,labindex+1,str2num(strcat(int2str(labindex),int2str(7),int2str(t))))
 end




   if labindex>Nfirstlab

for j=Jmin(labindex)-3:Jmax(labindex)+3
for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
a7sentdown(j,k)=N(j,Kmin(labindex)+k);
end
end


labSend(a7sentdown,labindex-1,str2num(strcat(int2str(labindex),int2str(7),int2str(t))))
    end


 
  
  if labindex<Nlastlab

a7receivedfromupstairs=labReceive(labindex+1,str2num(strcat(int2str(labindex+1),int2str(7),int2str(t))));

for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
N(j,Kmax(labindex)-1+k)=a7receivedfromupstairs(j,k);

          end
end

  end

    if labindex>Nfirstlab

a7receivedfromdownstairs=labReceive(labindex-1,str2num(strcat(int2str(labindex-1),int2str(7),int2str(t))));


for j=Jmin(labindex)-3:Jmax(labindex)+3
          for k=1:4%Kmin(labindex)+1:Kmin(labindex)+4
              
N(j,Kmin(labindex)-k+1)=a7receivedfromdownstairs(j,k);

          end
end
    end



 if mod(t,50)==0
   tindex=t/50;  
for j=1:Jmax(labindex)-Jmin(labindex)-1
for k=1:Kmax(labindex)-Kmin(labindex)-1
    
      Nwithtime(j,k,tindex)=N(j+Jmin(labindex),k+Kmin(labindex));
     % filename=strcat(rejimissaxeli,'lab',int2str(labindex),'t',int2str(t));
      %parsave(filename,N);
     
end
end
 end

end

filename=strcat(rejimissaxeli,'lab',int2str(labindex));
 
parsave(filename,Nwithtime);
end





   % delete(gcp)