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

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

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


Revision 1.7 - (show annotations) (download)
Mon May 2 19:35:38 2011 UTC (13 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62x
Changes since 1.6: +19 -6 lines
input OB files with 2 identical time-records, except U-West/U-East
 with +/- 1 per-cent.

1 % This is a matlab script that generates the input data
2
3 % $Header: /u/gcmpack/MITgcm/verification/exp4/input/gendata.m,v 1.6 2010/04/06 21:04:45 jmc Exp $
4 % $Name: $
5
6 % Dimensions of grid
7 nx=80;
8 ny=42;
9 nz=8;
10 % Nominal depth of model (meters)
11 H=4500;
12 % Scale of bump (m)
13 L=25e3;
14 % Height of bump (m)
15 dh=0.90*H;
16 % Horizontal resolution (m)
17 dx=5e3;
18 % Rotation
19 f=1e-4;
20 % Stratification
21 N=1.5 * f*L/H;
22
23 % Gravity
24 g=9.81;
25 % E.O.S.
26 alpha=2.e-4;
27
28 Tz=N^2/(g*alpha);
29 fprintf(' Tz= %e ;',Tz);
30
31 dz=H/nz;
32 fprintf(' delZ = %d * %7.6g\n',nz,dz);
33
34 x=(1:nx)*dx;x=x-mean(x);
35 y=(1:ny)*dx;y=y-mean(y);
36 z=-dz/2:-dz:-H;
37
38 [Y,X]=meshgrid(y,x);
39
40 % Temperature profile
41 fprintf('Tref ='); fprintf(' %8.6g,',Tz*z-mean(Tz*z)); fprintf('\n');
42
43 ieee='b';
44 prec='real*8';
45
46 % Gaussian bump
47 h=-H+dh*exp( -(X.^2+Y.^2)/(2*(L^2)) );
48 fid=fopen('topog.bump','w',ieee); fwrite(fid,h,prec); fclose(fid);
49
50 % $$$ % Side walls + bump
51 % $$$ h(:,1)=0;
52 % $$$ h(:,ny)=0;
53 % $$$ fid=fopen('topog.bumpchannel','w',ieee); fwrite(fid,h,prec); fclose(fid);
54
55 % $$$ % Simple channel
56 % $$$ h(:,1)=0;
57 % $$$ h(:,2:ny-1)=-H;
58 % $$$ h(:,ny)=0;
59 % $$$ fid=fopen('topog.channel','w',ieee); fwrite(fid,h,prec); fclose(fid);
60
61 % initial fields for salinity
62 si = 35;
63 %fid=fopen('S.init','w',ieee); fwrite(fid,si*ones(nx,ny,nz),prec); fclose(fid);
64
65 % open boundary conditions;
66 u0 = .25;
67 s0 = si+1;
68 w0= 1.e-3;
69
70 % create two time slabs for testing
71 uMerid = cat(3,u0*ones(nx,nz),u0*ones(nx,nz));
72 uZonal = cat(3,u0*ones(ny,nz), zeros(ny,nz));
73 sZonal = cat(3,s0*ones(ny,nz),s0*ones(ny,nz));
74
75 %- time varying fraction = 1 % of full velocity
76 du=u0*0.01;
77 uWest = cat(3,(u0+du)*ones(ny,nz),(u0-du)*ones(ny,nz));
78 uEast = cat(3,(u0-du)*ones(ny,nz),(u0+du)*ones(ny,nz));
79
80 % to test Non-Hydrostatic OBCS:
81 w1=[0:nz-1]*pi/nz; w1=-w0*sin(w1);
82 wZonal = cat (3, ones(ny,1)*w1, zeros(ny,nz));
83
84 fid=fopen('OBmeridU.bin','w',ieee); fwrite(fid,uMerid,prec); fclose(fid);
85 %fid=fopen('OBzonalU.bin','w',ieee); fwrite(fid,uZonal,prec); fclose(fid);
86 fid=fopen('OB_WestU.bin','w',ieee); fwrite(fid,uWest ,prec); fclose(fid);
87 fid=fopen('OB_EastU.bin','w',ieee); fwrite(fid,uEast ,prec); fclose(fid);
88 fid=fopen('OBzonalS.bin','w',ieee); fwrite(fid,sZonal,prec); fclose(fid);
89 fid=fopen('OBzonalW.bin','w',ieee); fwrite(fid,wZonal,prec); fclose(fid);
90
91 %- rbcs mask & restauring tracer field:
92 msk=ones(nx,ny,nz);
93 xMx=max(x);
94 shapeX=(x-xMx)/dx;
95 shapeX=exp(shapeX*2/3);
96
97 [I]=find(shapeX < 5.e-3); fprintf('zero out rbc-mask up to i= %i\n',max(I));
98 shapeX(I)=0.;
99 var=shapeX'*ones(1,ny*nz);
100 fid=fopen('rbcs_mask.bin','w',ieee); fwrite(fid,var,prec); fclose(fid);
101
102 tr1=(si+s0)/2;
103 var=tr1*ones(nx,ny,nz);
104 fid=fopen('rbcs_Tr1_fld.bin','w',ieee); fwrite(fid,var,prec); fclose(fid);
105

  ViewVC Help
Powered by ViewVC 1.1.22