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	cnh	1.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