/[MITgcm]/MITgcm/verification/advect_xz/input/gendata.m
ViewVC logotype

Contents of /MITgcm/verification/advect_xz/input/gendata.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (show annotations) (download)
Mon Dec 5 22:03:52 2011 UTC (12 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63g, checkpoint64, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.1: +15 -4 lines
- scale down the depth by 1/100 (from 200 km to just 2.km)
- add a 2nd initial velocity file with divergent transport

1 % $Header: /u/gcmpack/MITgcm/verification/advect_xz/input/gendata.m,v 1.1 2006/01/03 19:04:51 jmc Exp $
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.e2;
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=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=H*psi1.*psi2;
56
57 uTrans=psi(:,[1:nz])-psi(:,[2:nz+1]);
58 u0=uTrans(1:nx,:).*rhFc/dz;
59
60 %- 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 if w2file,
67 uname='Uvel.bin';
68 fid=fopen(uname,'w',ieee); fwrite(fid,u0,prec); fclose(fid);
69 fprintf([wf,uname,'\n']);
70 uname='Udiv.bin';
71 fid=fopen(uname,'w',ieee); fwrite(fid,u0+u1,prec); fclose(fid);
72 fprintf([wf,uname,'\n']);
73 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 %imagesc(xf,zf,psi1'); set(gca,'YDir','normal'); colorbar;
101 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 %imagesc(xf(1:nx),zc,u1'); set(gca,'YDir','normal'); colorbar;
106 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

  ViewVC Help
Powered by ViewVC 1.1.22