| 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 |
Parent Directory
| 
Revision Log
| 
Revision Graph