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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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 labindexNfirstlab 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)