
clear

%--------
nx = 2000;
nz = 350;

dx = 30;
dz = 10;

dz = dz * ones(nz,1);
dz(251:350) = dz(250) * 1.02.^(1:100);
%--------

%--------
Lx = nx*dx;
Hz = sum(dz);

x  = (dx/2:dx:Lx-dx/2);
z  = - Hz + cumsum(dz) - dz/2;

[x,z] = meshgrid(x,z);
%--------

%--------
n = 10; %nb de bosses
xbot = x(1,:);
Hbot = 110*(sin(n*2*pi.*xbot./(60000))+1)- Hz;
%Hbot = 110*(sin(2*pi*0.5/(6*10^3).*xbot(1:end))+1)- Hz;
%--------

%--------
U = 0.025*ones(nz,nx);
V = zeros(nz,nx);
T = (1e-3)^2 *(z+Hz) /9.81/2e-4;
R = T(:,1);
%--------

%--------
U  = flipud(U)';
V  = flipud(V)';
T  = flipud(T)';
R  = flipud(R);
dz = flipud(dz);

fid = fopen('U.init','w','b'); 
fwrite(fid,U,'real*8'); 
fclose(fid);

fid = fopen('V.init','w','b'); 
fwrite(fid,V,'real*8'); 
fclose(fid);

fid = fopen('T.init','w','b'); 
fwrite(fid,T,'real*8'); 
fclose(fid);

save('R.init','R','-ascii');

fid = fopen('topog.init','w','b'); 
fwrite(fid,Hbot(1,:)','real*8'); 
fclose(fid);

fid=fopen('delZ.init','w','b'); 
fwrite(fid,dz,'real*8'); 
fclose(fid);

U  = flipud(U');
V  = flipud(V');
T  = flipud(T');
R  = flipud(R);
dz = flipud(dz);
%--------

