Surface gravity water waves 151
c=sqrt(g*h(1));
zb=0;
zbn=0;
zbo=0;
k=0;
% Main program;
for k=1:nt
k=k+1;
zbn=zbmax*k*Dt/nd;
if k*Dt>nd
zbn=zbmax;
end
% Movement pattern of the sea floor;
for i=nx2-2:nx2+2
bed(i)=(zbn-2*zb+zbo);
end
for i=2:nx-1
tt=r*((h(i+1)+h(i))*(z(i+1)-z(i))-(h(i)+h(i-1))*
(z(i)-z(i-1)));
zn(i)=2*z(i)-zo(i)+tt-Dt*fb(i)*(z(i)-zo(i))+bed(i);
end
zn(nx)=z(nx)+sqrt(g*h(nx))*Dt/Dx*(z(nx-1)-z(nx));
zn(1)=z(1)+sqrt(g*h(1))*Dt/Dx*(z(2)-z(1));
for i=1:nx
zo(i)=z(i);
z(i)=zn(i);
% Wave data at half-time of the simulation period;
% Note that z1 is recorded only if nt is an even
number;
if k==nt/2
for i=1:nx
z1(i)=z(i);
end
end
end
zbo=zb;
zb=zbn;
end
plot(1:nx,z,'','LineWidth',1.5)
hold on
plot(1:nx,z1,'r','LineWidth',1.5)
xlabel('Number of spatial steps x 2 [m]')
ylabel('Wave height [m]')
text(70,0.41,'Time=nt*Dt/2')
text(31,0.43,'Time=nt*Dt')
text(58,-0.035,'Location of sea floor movement')