1 |
jmc |
1.1 |
% $Header: $ |
2 |
|
|
% $Name: $ |
3 |
|
|
|
4 |
|
|
% This is a matlab script that generates the input data |
5 |
|
|
prec='real*8'; |
6 |
|
|
ieee='b'; |
7 |
|
|
w2file=1; %- do write to file (1=y/0=no) |
8 |
|
|
|
9 |
|
|
%- Dimensions of grid |
10 |
|
|
nx=20; |
11 |
|
|
ny=1; |
12 |
|
|
nz=20; |
13 |
|
|
|
14 |
|
|
%- Horizontal & vertical resolution (m): |
15 |
|
|
dx=1.e4; dy=dx; dz=1.e4; |
16 |
|
|
% full size of the domain: |
17 |
|
|
Lx=nx*dx; |
18 |
|
|
H=nz*dz ; |
19 |
|
|
|
20 |
|
|
%- grid point coordinate : |
21 |
|
|
xf=[0:nx]*dx; |
22 |
|
|
xc=(xf(1:nx)+xf(2:nx+1))/2; |
23 |
|
|
zc=-dz/2:-dz:-H; |
24 |
|
|
zf=0:-dz:-H; |
25 |
|
|
|
26 |
|
|
wf=' write file: '; |
27 |
|
|
%- bathymetry : |
28 |
|
|
zb=-H*ones(nx,1); |
29 |
|
|
zb=-H*(1.-xc'/Lx/2); |
30 |
|
|
|
31 |
|
|
if w2file, |
32 |
|
|
bname='bathy_slope.bin'; |
33 |
|
|
fid=fopen(bname,'w',ieee); fwrite(fid,zb,prec); fclose(fid); |
34 |
|
|
fprintf([wf,bname,'\n']); |
35 |
|
|
end |
36 |
|
|
|
37 |
|
|
%- partial cell filling (hFac): |
38 |
|
|
hFac=ones(nx,1)*zf(1,[2:nz+1]); |
39 |
|
|
hFac=max(hFac,zb*ones(1,nz)); |
40 |
|
|
hFac=ones(nx,1)*zf(1,[1:nz])-hFac; |
41 |
|
|
hFac=hFac/dz; hFac=max(0,hFac); |
42 |
|
|
|
43 |
|
|
rhFc=reshape(hFac,[nx*nz 1]); |
44 |
|
|
[IK]=find(rhFc > 0); rhFc(IK)=1./rhFc(IK); |
45 |
|
|
rhFc=reshape(rhFc,[nx nz]); |
46 |
|
|
|
47 |
|
|
%- horizontal flow field is defined from stream-function psi : |
48 |
|
|
psi=zeros(nx+1,nz+1); |
49 |
|
|
psi1=Lx*sin(pi*xf'/Lx)*ones(1,nz+1); |
50 |
|
|
zu=zeros(nx+1,1); zu(2:nx)=max(zb(2:nx),zb(1:nx-1)); |
51 |
|
|
rHu=zeros(nx+1,1); rHu(2:nx)=-1./zu(2:nx); |
52 |
|
|
psi2=max( ones(nx+1,1)*zf, zu*ones(1,nz+1) ); |
53 |
|
|
psi2=psi2.*(rHu*ones(1,nz+1)); |
54 |
|
|
psi2=sin(pi*psi2); |
55 |
|
|
psi=psi1.*psi2; |
56 |
|
|
|
57 |
|
|
uTrans=psi(:,[1:nz])-psi(:,[2:nz+1]); |
58 |
|
|
u0=uTrans(1:nx,:).*rhFc/dz; |
59 |
|
|
|
60 |
|
|
if w2file, |
61 |
|
|
uname='Uvel.bin'; |
62 |
|
|
fid=fopen(uname,'w',ieee); fwrite(fid,u0,prec); fclose(fid); |
63 |
|
|
fprintf([wf,uname,'\n']); |
64 |
|
|
end |
65 |
|
|
|
66 |
|
|
%- Initial Temperature : Gaussian bump : |
67 |
|
|
t0=zeros(nx,nz); |
68 |
|
|
% position of the center |
69 |
|
|
xM=xf(5); zM=zf(6); |
70 |
|
|
% size of the patch |
71 |
|
|
dX=2*dx; dH=2*dz; |
72 |
|
|
% amplitude: |
73 |
|
|
tA=exp(-0.5*(35/20)^2); |
74 |
|
|
% |
75 |
|
|
r1=(xc'-xM)/dX; r2=(zc-zM)/dH; |
76 |
|
|
r1=r1.*r1; r2=r2.*r2; |
77 |
|
|
rD=r1*ones(1,nz)+ones(nx,1)*r2; |
78 |
|
|
t0=tA*exp(-rD/2); |
79 |
|
|
|
80 |
|
|
if w2file, |
81 |
|
|
tname='Tini_G.bin'; |
82 |
|
|
fid=fopen(tname,'w',ieee); fwrite(fid,t0,prec); fclose(fid); |
83 |
|
|
fprintf([wf,tname,'\n']); |
84 |
|
|
end |
85 |
|
|
|
86 |
|
|
%return |
87 |
|
|
|
88 |
|
|
%- make plots to check: |
89 |
|
|
figure(1);clf; |
90 |
|
|
subplot(211) |
91 |
|
|
imagesc(xf,zf,psi'); set(gca,'YDir','normal'); colorbar; |
92 |
|
|
hold on ; plot(xc,zb,'b-'); hold off ; grid |
93 |
|
|
title('Stream-Function [m^2/s]'); |
94 |
|
|
subplot(212) |
95 |
|
|
imagesc(xf(1:nx),zc,u0'); set(gca,'YDir','normal'); colorbar; |
96 |
|
|
hold on ; plot(xc,zb,'r-'); hold off ; grid |
97 |
|
|
title('U velocity [m/s]'); |
98 |
|
|
|
99 |
|
|
figure(2);clf; |
100 |
|
|
subplot(211) |
101 |
|
|
imagesc(xc,zc,rhFc'); set(gca,'YDir','normal'); colorbar; |
102 |
|
|
hold on ; plot(xc,zb,'k-'); hold off ; grid |
103 |
|
|
title('Recip hFac'); |
104 |
|
|
subplot(212) |
105 |
|
|
imagesc(xc,zc,t0'); set(gca,'YDir','normal'); colorbar; |
106 |
|
|
hold on ; plot(xc,zb,'r-'); hold off ; grid |
107 |
|
|
title('Initial Theta'); |
108 |
|
|
|
109 |
|
|
return |