/[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.8 - (show annotations) (download)
Thu May 19 00:03:06 2011 UTC (12 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint62z, checkpoint62y, checkpoint64, checkpoint65, checkpoint63, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint65o, 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, HEAD
Changes since 1.7: +8 -1 lines
add OBC files for Eta at Western & Eastern OB (to test prescribing Eta
 in NonLin-FreeSurf formulation)

1 % This is a matlab script that generates the input data
2
3 % $Header: /u/gcmpack/MITgcm/verification/exp4/input/gendata.m,v 1.7 2011/05/02 19:35:38 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 % to test prescribing Eta in NonLin-FreeSurf formulation:
85 et1=0.1;
86 etWest = cat(2,+et1*ones(ny,1),-et1*ones(ny,1));
87 etEast = cat(2,-et1*ones(ny,1),+et1*ones(ny,1));
88
89 fid=fopen('OBmeridU.bin','w',ieee); fwrite(fid,uMerid,prec); fclose(fid);
90 %fid=fopen('OBzonalU.bin','w',ieee); fwrite(fid,uZonal,prec); fclose(fid);
91 fid=fopen('OB_WestU.bin','w',ieee); fwrite(fid,uWest ,prec); fclose(fid);
92 fid=fopen('OB_EastU.bin','w',ieee); fwrite(fid,uEast ,prec); fclose(fid);
93 fid=fopen('OBzonalS.bin','w',ieee); fwrite(fid,sZonal,prec); fclose(fid);
94 fid=fopen('OBzonalW.bin','w',ieee); fwrite(fid,wZonal,prec); fclose(fid);
95 fid=fopen('OB_WestH.bin','w',ieee); fwrite(fid,etWest,prec); fclose(fid);
96 fid=fopen('OB_EastH.bin','w',ieee); fwrite(fid,etEast,prec); fclose(fid);
97
98 %- rbcs mask & restauring tracer field:
99 msk=ones(nx,ny,nz);
100 xMx=max(x);
101 shapeX=(x-xMx)/dx;
102 shapeX=exp(shapeX*2/3);
103
104 [I]=find(shapeX < 5.e-3); fprintf('zero out rbc-mask up to i= %i\n',max(I));
105 shapeX(I)=0.;
106 var=shapeX'*ones(1,ny*nz);
107 fid=fopen('rbcs_mask.bin','w',ieee); fwrite(fid,var,prec); fclose(fid);
108
109 tr1=(si+s0)/2;
110 var=tr1*ones(nx,ny,nz);
111 fid=fopen('rbcs_Tr1_fld.bin','w',ieee); fwrite(fid,var,prec); fclose(fid);
112

  ViewVC Help
Powered by ViewVC 1.1.22