% 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);