internal_wave/input/gendata.m

% This is a matlab script that generates the input data
% variable x resolution
prec='real*8';
ieee='b';


% Dimensions of grid
nx=60;
ny=3;
nz=20;
% Nominal depth of model (meters)
H=200.0;
% Size of domain
Lx=13.3e3;

% Horizontal resolution (m)
dx=zeros(nx,1);
for i=1:nx
dx(i) = Lx/(nx+1);
end

dy = Lx/nx
% Stratification
gravity=9.81;
talpha=2.0e-4;
N2=1e-6;

Tz=N2/(gravity*talpha);

dz=H/nz;
sprintf('delZ = %d * %7.6g,',nz,dz)

x=zeros(nx,1);
x(1) = dx(1);
for i=2:nx
x(i)=x(i-1) + dx(i);
end
z=-dz/2:-dz:-H;

%[Y,X]=meshgrid(y,x);

% Temperature profile
Tref=Tz*z-mean(Tz*z);
[sprintf('Tref =') sprintf(' %8.6g,',Tref)]
t=0.0*rand([nx,ny,nz]);
for k=1:nz
t(:,:,k) = t(:,:,k) + Tref(k);
end
fid=fopen('T.init','w',ieee); fwrite(fid,t,prec); fclose(fid);

% Sloping channel 
slope=0.0
offset=2.5e3;
dmax=-40.0;
d=0.0*rand([nx,ny]);
for i=1:nx
for j=1:ny
d(i,j) = -H;
d(i,j) = d(i,j) + slope*(x(i) - offset);
if d(i,j) < -H ;
d(i,j) = -H;
end
if d(i,j) > dmax;
d(i,j) = dmax;
end
end
end
%d(nx,:)=0.0;
fid=fopen('topog.slope','w',ieee); fwrite(fid,d,prec); fclose(fid);
plot(x,d(:,2))

fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);

%convex slope
nxslope=(dmax + H)/(slope)
d1=zeros(nx,ny);
hamp=(H-dmax)/5.0
pi=4.0*atan(1.0)
for i=1:nx
for j=1:ny
if x(i) < (offset + nxslope)
if x(i) < offset
d1(i,j) = d(i,j);
else
d1(i,j) = -H;
d1(i,j) = d(i,j) + hamp*sin(pi*(x(i)-offset)/nxslope);
if d1(i,j) < -H ;
d1(i,j) = -H;
end
if d1(i,j) > dmax;
d1(i,j) = dmax;
end
end
else
d1(i,j) = d(i,j);
end
end
end
%d1(end-1:end,:)=d1(1:2,:); % debug by aja
fid=fopen('topog.convex','w',ieee); fwrite(fid,d1,prec); fclose(fid);
hold on
plot(x,d1(:,2),'g')

%convex slope
d2=zeros(nx,ny);
for i=1:nx
for j=1:ny
if x(i) < (offset + nxslope)
if x(i) < offset
d2(i,j) = d(i,j);
else
d2(i,j) = -H;
d2(i,j) = d(i,j) - hamp*sin(pi*(x(i)-offset)/nxslope);
if d2(i,j) < -H ;
d2(i,j) = -H;
end
if d2(i,j) > dmax;
d2(i,j) = dmax;
end
end
else
d2(i,j) = d(i,j);
end
end
end
%d2(end-1:end,:)=d2(1:2,:); % debug by aja
fid=fopen('topog.concave','w',ieee); fwrite(fid,d2,prec); fclose(fid);
hold on
plot(x,d2(:,2),'r')


fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);
1	% This is a matlab script that generates the input data
2	% variable x resolution
3	prec='real*8';
4	ieee='b';
5
6
7	% Dimensions of grid
8	nx=60;
9	ny=3;
10	nz=20;
11	% Nominal depth of model (meters)
12	H=200.0;
13	% Size of domain
14	Lx=13.3e3;
15
16	% Horizontal resolution (m)
17	dx=zeros(nx,1);
18	for i=1:nx
19	dx(i) = Lx/(nx+1);
20	end
21
22	dy = Lx/nx
23	% Stratification
24	gravity=9.81;
25	talpha=2.0e-4;
26	N2=1e-6;
27
28	Tz=N2/(gravity*talpha);
29
30	dz=H/nz;
31	sprintf('delZ = %d * %7.6g,',nz,dz)
32
33	x=zeros(nx,1);
34	x(1) = dx(1);
35	for i=2:nx
36	x(i)=x(i-1) + dx(i);
37	end
38	z=-dz/2:-dz:-H;
39
40	%[Y,X]=meshgrid(y,x);
41
42	% Temperature profile
43	Tref=Tzz-mean(Tzz);
44	[sprintf('Tref =') sprintf(' %8.6g,',Tref)]
45	t=0.0*rand([nx,ny,nz]);
46	for k=1:nz
47	t(:,:,k) = t(:,:,k) + Tref(k);
48	end
49	fid=fopen('T.init','w',ieee); fwrite(fid,t,prec); fclose(fid);
50
51	% Sloping channel
52	slope=0.0
53	offset=2.5e3;
54	dmax=-40.0;
55	d=0.0*rand([nx,ny]);
56	for i=1:nx
57	for j=1:ny
58	d(i,j) = -H;
59	d(i,j) = d(i,j) + slope*(x(i) - offset);
60	if d(i,j) < -H ;
61	d(i,j) = -H;
62	end
63	if d(i,j) > dmax;
64	d(i,j) = dmax;
65	end
66	end
67	end
68	%d(nx,:)=0.0;
69	fid=fopen('topog.slope','w',ieee); fwrite(fid,d,prec); fclose(fid);
70	plot(x,d(:,2))
71
72	fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);
73
74	%convex slope
75	nxslope=(dmax + H)/(slope)
76	d1=zeros(nx,ny);
77	hamp=(H-dmax)/5.0
78	pi=4.0*atan(1.0)
79	for i=1:nx
80	for j=1:ny
81	if x(i) < (offset + nxslope)
82	if x(i) < offset
83	d1(i,j) = d(i,j);
84	else
85	d1(i,j) = -H;
86	d1(i,j) = d(i,j) + hampsin(pi(x(i)-offset)/nxslope);
87	if d1(i,j) < -H ;
88	d1(i,j) = -H;
89	end
90	if d1(i,j) > dmax;
91	d1(i,j) = dmax;
92	end
93	end
94	else
95	d1(i,j) = d(i,j);
96	end
97	end
98	end
99	%d1(end-1:end,:)=d1(1:2,:); % debug by aja
100	fid=fopen('topog.convex','w',ieee); fwrite(fid,d1,prec); fclose(fid);
101	hold on
102	plot(x,d1(:,2),'g')
103
104	%convex slope
105	d2=zeros(nx,ny);
106	for i=1:nx
107	for j=1:ny
108	if x(i) < (offset + nxslope)
109	if x(i) < offset
110	d2(i,j) = d(i,j);
111	else
112	d2(i,j) = -H;
113	d2(i,j) = d(i,j) - hampsin(pi(x(i)-offset)/nxslope);
114	if d2(i,j) < -H ;
115	d2(i,j) = -H;
116	end
117	if d2(i,j) > dmax;
118	d2(i,j) = dmax;
119	end
120	end
121	else
122	d2(i,j) = d(i,j);
123	end
124	end
125	end
126	%d2(end-1:end,:)=d2(1:2,:); % debug by aja
127	fid=fopen('topog.concave','w',ieee); fwrite(fid,d2,prec); fclose(fid);
128	hold on
129	plot(x,d2(:,2),'r')
130
131
132	fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);