front_relax/input/gendata.m

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

% Dimensions/resolution of grid (resolution/depth in m)
nx=1; dx=10e3;
ny=32; dy=10e3;
dz=[50 50 55 60 65 70 80 95 120 155 200 260 320 400 480];
nz=prod(size(dz));

% Nominal depth of model (meters)
Ho=sum(dz);

% Topography (channel)
H=-Ho*ones(nx,ny);
H(:,end)=0;  % Solid wall in North
fid=fopen('topog.bin','w',ieee); fwrite(fid,H,prec); fclose(fid);

% Size of domain (m)
Lx=dx*nx;
Ly=dy*(ny-1);  % Solid wall in North

% Variable resolution
y=(1:ny)/(ny-0)-0.5;
dy=1-0.3*exp( -(5*y).^2 );
%dy=ones(1,ny);              % Constant resolution
dy=dy/sum(dy(1:ny-1))*Ly;
fid=fopen('dy.bin','w',ieee); fwrite(fid,dy,prec); fclose(fid);

% Coordinates
xc=((1:nx)-0.5)*dx;
yf=-Ly/2+[0 cumsum(dy)];
yc=(yf(1:end-1)+yf(2:end))/2;
zf=[0 -cumsum(dz)];
zc=(zf(1:end-1)+zf(2:end))/2;
[X,Y]=ndgrid(xc,yc);

% Stratification
fo=1e-4;
alpha=2e-4;
gravity=9.81;
N_over_f=20;
N2=(N_over_f*fo)^2;
dTdz=N2/(alpha*gravity);

slope=1e-3;
dTdy=-slope*dTdz;

for k=1:nz,
 T(:,:,k)=20 +dTdz*zc(k) ...
             +dTdy*Ly*sin(pi*Y/Ly)...
             *exp(-(3*zc(k)/Ho)^2);
end
T(:,end,:)=0;
fid=fopen('theta.bin','w',ieee); fwrite(fid,T,prec); fclose(fid);

% Salt (passive tracer)
for k=1:nz,
%S(:,:,k)=ones(nx,ny)+zc(k)/Ho;  % Linear with z
%S(:,:,k)=exp(2*zc(k)/Ho);       % Exponential with z
%S(:,:,k)=exp(2*Y/Ly);           % Exponential with y
 S(:,:,k)=exp(-(2*Y/Ly).^2);           % Exponential with y
end
S(:,end,:)=0;
fid=fopen('salt.bin','w',ieee); fwrite(fid,S,prec); fclose(fid);
1	% This is a matlab script that generates the input data
2	prec='real*8';
3	ieee='b';
4
5	% Dimensions/resolution of grid (resolution/depth in m)
6	nx=1; dx=10e3;
7	ny=32; dy=10e3;
8	dz=[50 50 55 60 65 70 80 95 120 155 200 260 320 400 480];
9	nz=prod(size(dz));
10
11	% Nominal depth of model (meters)
12	Ho=sum(dz);
13
14	% Topography (channel)
15	H=-Ho*ones(nx,ny);
16	H(:,end)=0; % Solid wall in North
17	fid=fopen('topog.bin','w',ieee); fwrite(fid,H,prec); fclose(fid);
18
19	% Size of domain (m)
20	Lx=dx*nx;
21	Ly=dy*(ny-1); % Solid wall in North
22
23	% Variable resolution
24	y=(1:ny)/(ny-0)-0.5;
25	dy=1-0.3exp( -(5y).^2 );
26	%dy=ones(1,ny); % Constant resolution
27	dy=dy/sum(dy(1:ny-1))*Ly;
28	fid=fopen('dy.bin','w',ieee); fwrite(fid,dy,prec); fclose(fid);
29
30	% Coordinates
31	xc=((1:nx)-0.5)*dx;
32	yf=-Ly/2+[0 cumsum(dy)];
33	yc=(yf(1:end-1)+yf(2:end))/2;
34	zf=[0 -cumsum(dz)];
35	zc=(zf(1:end-1)+zf(2:end))/2;
36	[X,Y]=ndgrid(xc,yc);
37
38	% Stratification
39	fo=1e-4;
40	alpha=2e-4;
41	gravity=9.81;
42	N_over_f=20;
43	N2=(N_over_f*fo)^2;
44	dTdz=N2/(alpha*gravity);
45
46	slope=1e-3;
47	dTdy=-slope*dTdz;
48
49	for k=1:nz,
50	T(:,:,k)=20 +dTdz*zc(k) ...
51	+dTdyLysin(pi*Y/Ly)...
52	exp(-(3zc(k)/Ho)^2);
53	end
54	T(:,end,:)=0;
55	fid=fopen('theta.bin','w',ieee); fwrite(fid,T,prec); fclose(fid);
56
57	% Salt (passive tracer)
58	for k=1:nz,
59	%S(:,:,k)=ones(nx,ny)+zc(k)/Ho; % Linear with z
60	%S(:,:,k)=exp(2*zc(k)/Ho); % Exponential with z
61	%S(:,:,k)=exp(2*Y/Ly); % Exponential with y
62	S(:,:,k)=exp(-(2*Y/Ly).^2); % Exponential with y
63	end
64	S(:,end,:)=0;
65	fid=fopen('salt.bin','w',ieee); fwrite(fid,S,prec); fclose(fid);