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=1; |
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=Tz*z-mean(Tz*z); |
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.03 |
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(:,1)) |
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) + hamp*sin(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 |
plot(x,d1(:,1),'g') |
102 |
|
103 |
%convex slope |
104 |
d2=zeros(nx,ny); |
105 |
for i=1:nx |
106 |
for j=1:ny |
107 |
if x(i) < (offset + nxslope) |
108 |
if x(i) < offset |
109 |
d2(i,j) = d(i,j); |
110 |
else |
111 |
d2(i,j) = -H; |
112 |
d2(i,j) = d(i,j) - hamp*sin(pi*(x(i)-offset)/nxslope); |
113 |
if d2(i,j) < -H ; |
114 |
d2(i,j) = -H; |
115 |
end |
116 |
if d2(i,j) > dmax; |
117 |
d2(i,j) = dmax; |
118 |
end |
119 |
end |
120 |
else |
121 |
d2(i,j) = d(i,j); |
122 |
end |
123 |
end |
124 |
end |
125 |
%d2(end-1:end,:)=d2(1:2,:); % debug by aja |
126 |
fid=fopen('topog.concave','w',ieee); fwrite(fid,d2,prec); fclose(fid); |
127 |
hold on |
128 |
plot(x,d2(:,1),'r') |
129 |
hold off |
130 |
|
131 |
|
132 |
fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid); |