Contaminant and sediment transport by advection and diffusion 209
else
advx=ur*(c(i+1,j)-c(i,j))/Dd;
end
advy=0;
if v+vo>0
advy=v*(c(i,j)-c(i,j-1))/Dd;
else
advy=vo*(c(i,j+1)-c(i,j))/Dd;
end
co=c(i,j+1);
if ho<0.1
co=c(i,j);
end
cu=c(i,j-1);
if hu<0.1
cu=c(i,j);
end
cr=c(i+1,j);
if hr<0.1
cr=c(i,j);
end
cl=c(i-1,j);
if hl<0.1
cl=c(i,j);
end
diff1=ed/h*((co-c(i,j))*(h+ho)-(c(i,j)cu)*(h+hu))/2/Dd^2;
diff2=ed/h*((cr-c(i,j))*(hr+h)-(c(i,j)cl)*(h+hl))/2/Dd^2;
cn(i,j)=c(i,j)+Dt*(-advx-advy+diff1+diff2dc*c(i,j));
end
end
end
% Boundary conditions;
for i=2:nx-1
cn(i,2)=2*cn(i,3)-c(i,4);
end
%
for j=2:ny-1
%
cn(2,j)=2*cn(3,j)-cn(4,j);
%
cn(nx-2,j)=2*cn(nx-3,j)-cn(nx-4,j);
%
end
% Point source pollution;
cn(imm,jmm)=csource;
% Updating the contaminant concentration values;
for j=1:ny
for i=1:nx
c(i,j)=cn(i,j);
if cn(i,j)>cmax(i,j)