sigma_testing/input-sigma/topo_bump.m

% 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


1	% Set up bump spread over 200km in 2000km channel going from -2km to 400m
2	% with 1km horizontal resolution.
3	Lc=2000e3;
4	Lx=1000e3;
5	Lb=200e3;
6	Hlow=-2000;
7	Htop=-1800;
8	%Htop=-1000;
9	s_top=-800;
10	%Htop=-400;
11	nx=1000;
12	clear d
13	dx=Lx/nx.*ones(nx,1);
14	Bmid=Lx/2;
15	offset=Bmid;
16	x=zeros(nx,1);
17	x(1)=dx(1)*0.5;
18	for i=2:nx
19	t(i-1)=tanh( (Lx-x(i-1)-offset)/(Lb*0.4) );
20	d(i-1) = -(Htop-Hlow)/2*(t(i-1) + 1) - Hlow;
21	x(i)=x(i-1)+dx(i-1)0.5+dx(i)0.5;
22	end
23	t(nx)=tanh( (Lx-x(nx)-offset)/(Lb*0.5) );
24	d(nx) = -(Htop-Hlow)/2*(t(nx) + 1) - Hlow;
25	dfull=[-d(end:-1:1) -d];
26	d=dfull;
27	x=ones(1,nx2)dx(1);x=cumsum(x)-dx(1)/2;
28	plot(x,d);axis([0 max(x) -2000 0]);
29
30	fid=fopen('topog_bump.data','w','ieee-be');
31	fwrite(fid,d,'float64');
32	fclose(fid);
33
34	% Define flow
35	vmax=0.1;
36	[mv,im]=max(d);
37	tm=abs(d(im))*vmax;
38	vofx=tm./abs(d);
39	fid=fopen('vbaro_bump.data','w','ieee-be');
40	fwrite(fid,vofx,'float64');
41	fclose(fid);
42
43	% Now try a couple of plain sigma surfaces
44	% sigma=z/-H, sigma = 0 @ z = 0; sigma = 1 @ z = -H\|ps;
45	sigvals=[ 0 0.2 0.4 0.6 0.8 1];
46	sigvals=[0:0.1:1];
47	clear rC, rW;
48	nlevs=length(sigvals)-1;
49	rC=zeros(length(x),nlevs);
50	for i=1:nlevs
51	s=(sigvals(i)+sigvals(i+1))*0.5;
52	rC(:,i)=s.*(d);
53	end
54	rW=zeros(length(x),length(sigvals));
55	for i=1:length(sigvals)
56	rW(:,i)=sigvals(i).*(d);
57	end
58
59	% Now try some hybrid sigma - z\|p surfaces (actually its not hybrid here!)
60	H0=-s_top;
61	b=[ 1 0.8 0.6 0.4 0.2 0 0 0 0 0 ];
62	a=[ 0 0 0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 ]*H0;
63	nlevs=length(b)-1;
64	rCh=zeros(length(x),nlevs);
65	for i=1:nlevs
66	s_b=0.5*(b(i)+b(i+1));
67	s_a=0.5*(a(i)+a(i+1));
68	RCh(:,i)=s_b.(d)+(1-s_b)s_top+s_a;
69	end
70	for i=1:length(b)
71	RWh(:,i)=b(i).(d)+(1-b(i))s_top+a(i);
72	end
73
74