1 |
jmc |
1.2 |
% $Header: /u/gcmpack/MITgcm/verification/advect_xz/input/gendata.m,v 1.1 2006/01/03 19:04:51 jmc Exp $ |
2 |
jmc |
1.1 |
% $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 |
jmc |
1.2 |
dx=1.e4; dy=dx; dz=1.e2; |
16 |
jmc |
1.1 |
% 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 |
jmc |
1.2 |
psi1=sin(pi*xf'/Lx)*ones(1,nz+1); |
50 |
jmc |
1.1 |
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 |
jmc |
1.2 |
psi=H*psi1.*psi2; |
56 |
jmc |
1.1 |
|
57 |
|
|
uTrans=psi(:,[1:nz])-psi(:,[2:nz+1]); |
58 |
|
|
u0=uTrans(1:nx,:).*rhFc/dz; |
59 |
|
|
|
60 |
jmc |
1.2 |
%- add small, vertically uniform, divergent component: |
61 |
|
|
ud=psi1(:,1).*rHu; |
62 |
|
|
ud=ud*H/160.; |
63 |
|
|
u1=ud(1:nx,1)*ones(1,nz); |
64 |
|
|
u1(find(hFac==0.))=0.; |
65 |
|
|
|
66 |
jmc |
1.1 |
if w2file, |
67 |
|
|
uname='Uvel.bin'; |
68 |
|
|
fid=fopen(uname,'w',ieee); fwrite(fid,u0,prec); fclose(fid); |
69 |
|
|
fprintf([wf,uname,'\n']); |
70 |
jmc |
1.2 |
uname='Udiv.bin'; |
71 |
|
|
fid=fopen(uname,'w',ieee); fwrite(fid,u0+u1,prec); fclose(fid); |
72 |
|
|
fprintf([wf,uname,'\n']); |
73 |
jmc |
1.1 |
end |
74 |
|
|
|
75 |
|
|
%- Initial Temperature : Gaussian bump : |
76 |
|
|
t0=zeros(nx,nz); |
77 |
|
|
% position of the center |
78 |
|
|
xM=xf(5); zM=zf(6); |
79 |
|
|
% size of the patch |
80 |
|
|
dX=2*dx; dH=2*dz; |
81 |
|
|
% amplitude: |
82 |
|
|
tA=exp(-0.5*(35/20)^2); |
83 |
|
|
% |
84 |
|
|
r1=(xc'-xM)/dX; r2=(zc-zM)/dH; |
85 |
|
|
r1=r1.*r1; r2=r2.*r2; |
86 |
|
|
rD=r1*ones(1,nz)+ones(nx,1)*r2; |
87 |
|
|
t0=tA*exp(-rD/2); |
88 |
|
|
|
89 |
|
|
if w2file, |
90 |
|
|
tname='Tini_G.bin'; |
91 |
|
|
fid=fopen(tname,'w',ieee); fwrite(fid,t0,prec); fclose(fid); |
92 |
|
|
fprintf([wf,tname,'\n']); |
93 |
|
|
end |
94 |
|
|
|
95 |
|
|
%return |
96 |
|
|
|
97 |
|
|
%- make plots to check: |
98 |
|
|
figure(1);clf; |
99 |
|
|
subplot(211) |
100 |
jmc |
1.2 |
%imagesc(xf,zf,psi1'); set(gca,'YDir','normal'); colorbar; |
101 |
jmc |
1.1 |
imagesc(xf,zf,psi'); set(gca,'YDir','normal'); colorbar; |
102 |
|
|
hold on ; plot(xc,zb,'b-'); hold off ; grid |
103 |
|
|
title('Stream-Function [m^2/s]'); |
104 |
|
|
subplot(212) |
105 |
jmc |
1.2 |
%imagesc(xf(1:nx),zc,u1'); set(gca,'YDir','normal'); colorbar; |
106 |
jmc |
1.1 |
imagesc(xf(1:nx),zc,u0'); set(gca,'YDir','normal'); colorbar; |
107 |
|
|
hold on ; plot(xc,zb,'r-'); hold off ; grid |
108 |
|
|
title('U velocity [m/s]'); |
109 |
|
|
|
110 |
|
|
figure(2);clf; |
111 |
|
|
subplot(211) |
112 |
|
|
imagesc(xc,zc,rhFc'); set(gca,'YDir','normal'); colorbar; |
113 |
|
|
hold on ; plot(xc,zb,'k-'); hold off ; grid |
114 |
|
|
title('Recip hFac'); |
115 |
|
|
subplot(212) |
116 |
|
|
imagesc(xc,zc,t0'); set(gca,'YDir','normal'); colorbar; |
117 |
|
|
hold on ; plot(xc,zb,'r-'); hold off ; grid |
118 |
|
|
title('Initial Theta'); |
119 |
|
|
|
120 |
|
|
return |