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	adcroft	1.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);