% Set up bump spread over 200km in 2000km channel going from -2km to 400m 
% with 1km horizontal resolution.
Lc=2000e3;
Lx=1000e3;
Lb=200e3;
Hlow=-2000;
Htop=-1800;
%Htop=-1000;
s_top=-800;
%Htop=-400;
nx=1000;
clear d
dx=Lx/nx.*ones(nx,1);
Bmid=Lx/2;
offset=Bmid;
x=zeros(nx,1);
x(1)=dx(1)*0.5;
for i=2:nx
    t(i-1)=tanh( (Lx-x(i-1)-offset)/(Lb*0.4) );
    d(i-1) = -(Htop-Hlow)/2*(t(i-1) + 1) - Hlow;
    x(i)=x(i-1)+dx(i-1)*0.5+dx(i)*0.5;
end
t(nx)=tanh( (Lx-x(nx)-offset)/(Lb*0.5) );
d(nx) = -(Htop-Hlow)/2*(t(nx) + 1) - Hlow;
dfull=[-d(end:-1:1) -d];
d=dfull;
x=ones(1,nx*2)*dx(1);x=cumsum(x)-dx(1)/2;
plot(x,d);axis([0 max(x) -2000 0]);

fid=fopen('topog_bump.data','w','ieee-be');
fwrite(fid,d,'float64');
fclose(fid);

% Define flow
vmax=0.1;
[mv,im]=max(d);
tm=abs(d(im))*vmax;
vofx=tm./abs(d);
fid=fopen('vbaro_bump.data','w','ieee-be');
fwrite(fid,vofx,'float64');
fclose(fid);

% Now try a couple of plain sigma surfaces
% sigma=z/-H, sigma = 0 @ z = 0; sigma = 1 @ z = -H|ps;
sigvals=[ 0 0.2 0.4 0.6 0.8 1];
sigvals=[0:0.1:1];
clear rC, rW;
nlevs=length(sigvals)-1;
rC=zeros(length(x),nlevs);
for i=1:nlevs
    s=(sigvals(i)+sigvals(i+1))*0.5;
    rC(:,i)=s.*(d);
end
rW=zeros(length(x),length(sigvals));
for i=1:length(sigvals)
    rW(:,i)=sigvals(i).*(d);
end

% Now try some hybrid sigma - z|p surfaces (actually its not hybrid here!)
H0=-s_top;
b=[ 1  0.8  0.6  0.4  0.2 0   0   0   0   0   ];
a=[ 0    0    0  0.0  0.0 0.2 0.4 0.6 0.8 1.0 ]*H0;
nlevs=length(b)-1;
rCh=zeros(length(x),nlevs);
for i=1:nlevs
    s_b=0.5*(b(i)+b(i+1));
    s_a=0.5*(a(i)+a(i+1));
    RCh(:,i)=s_b.*(d)+(1-s_b)*s_top+s_a;
end
for i=1:length(b)
    RWh(:,i)=b(i).*(d)+(1-b(i))*s_top+a(i);
end