clear all; close all;
Drho=47.;
Rho0=1000.;
ai=0.027;
Dh=0.0035;
H=0.15;
zi=0.024;
dplus=zi;
dminus=H-zi;
W=sqrt(4.*dminus*dminus*dplus*dplus/3./(dminus-dplus)/ai);
slope=0.217;
Ls=H/slope;
Lr=0.05;
disp(Lr+Ls)
N=100;
xs=linspace(0,(Ls+Lr),N);
zs=linspace(0,H,N);
[x,z]=meshgrid(xs,zs);
%x=x';z=z';
%z=z';
%x=x';
    for i=1:N
        for j=1:N
        zeta(i,j)=2.*ai*sech(x(i,j)/2./W)*sech(x(i,j)/2./W);
        rhopri(i,j)=Drho*0.5*(1.+tanh((z(i,j)-zi-zeta(i,j))/Dh));
        Rho(i,j)=Rho0+rhopri(i,j);
        end
    end
Rho=flipud(Rho);
contourf(x,z,Rho);colorbar
axis([0 Ls+Lr 0.05 H])
shading flat
xlabel('x','fontsize',14);ylabel('y','fontsize',14)
title('Initial Condition','fontsize',14)
saveas(gca,'f1','jpg')