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

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

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


Revision 1.1 - (show annotations) (download)
Wed Feb 12 23:15:51 2014 UTC (8 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64u, checkpoint64w, checkpoint64v, 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, checkpoint65, HEAD
remove "genit.m" (does not correspond to this set-up) and move
 ../input.0.smag/genit.m to gendata.m

1 %function geninit()
2 % This function generates initial condition files for the mitgcm.
3 % if init_vel=1, then the front problem is used
4 % if init_vel=0, then resting (other than small symmetry-breaking
5 % random noise) initial conditions are used.
6
7 nx=50;
8 ny=nx/2+1;
9 nz=40;
10 init_vel=1;
11 dxspacing=1000;
12 dyspacing=dxspacing;
13
14 Lx=dxspacing*nx;
15 Ly=dyspacing*ny;
16
17 fprintf(' nx= %i , ny= %i , nz= %i ; dx=%6.1f , dy=%6.1f\n', ...
18 nx,ny,nz,dxspacing,dxspacing)
19
20 %-- Params
21 g=9.81;
22 tAlpha=-2e-4;
23 f0=7.29e-5;
24 rho=1035;
25 day=24*60^2;
26 prec='real*8';
27 ieee='b';
28
29 H=200; %Max depth
30
31 fprintf(' Lx=%6.1f , Ly=%6.1f , H=%6.1f\n',Lx,Ly,H);
32
33 %-- Grid: x
34 dx=ones(1,nx); % uniform resolution
35 dx=dx*Lx/sum(dx);
36 xf=cumsum([0 dx]); % Face x points
37 xc=(xf(1:end-1)+xf(2:end))/2; % Centered x points
38
39 %-- Grid: y
40 dy=ones(1,ny); % uniform resolution
41 dy=dy*Ly/sum(dy);
42 yf=cumsum([0 dy]); % Face y-points
43 yc=(yf(1:end-1)+yf(2:end))/2; %Centered y-points
44 L=yc(end)-yc(1); % this takes into account the wall of topography!!!
45
46 %-- Grid: z
47 dh=H/nz*ones(1,nz);
48 zf=-round(cumsum([0 dh])/1)*1; % Face z points
49 dh=-diff(zf);
50 zc=(zf(1:end-1)+zf(2:end))/2; % centered z points
51 nz=length(dh);
52 H=sum(dh);
53
54 [XT,YT,ZT]=ndgrid(xc,yc,zc); % This is the centered, temperature grid.
55 [XU,YU,ZU]=ndgrid(xf,yc,zc); % This is the U grid.
56 [XV,YV,ZV]=ndgrid(xc,yf,zc); % This is the V grid.
57 [XW,YW,ZW]=ndgrid(xc,yc,zf); % This is the W grid.
58 [XB,YB]=ndgrid(xc,yc); % This is the Bathymetry grid.
59
60 % Stratification Params
61 Dml=200; % Mixed-layer Depth
62 D=Dml;
63 Lf=2*1e3; % Half-Width of Front
64 y0=yc(round(length(yc)/2)); % centered location of Front
65
66 Ms=(2e-4)^2 ; %db/dy 2e-4 is about 1 K/50 km
67 Ns=9*Ms^2/f0^2 ; %db/dz
68
69 deltatheta=Lf*Ms/g/tAlpha;
70 T0=17;
71 dthetadz=-Ns/g/tAlpha;
72
73 fprintf(' Ms=%15.6e , Ns=%15.6e , T0= %6.3f\n',Ms,Ns,T0);
74 fprintf(' deltatheta=%15.6e , dthetadz=%15.6e\n',deltatheta,dthetadz);
75
76 %-- Bathymetry
77 % Bathymetry is on Xc, Yc
78 hh=ones(size(XB));
79 hh(:,end)=0*hh(:,end);
80 hh=-H*hh;
81
82 figure(1); clf
83 subplot(221)
84 pcolor(XB/1e3,YB/1e3,hh)
85 axis equal
86
87 subplot(223)
88 plot(xc/1e3,dx/1e3,'.'); xlabel('X (km)');ylabel('\Delta x (km)')
89 title('\Delta x')
90
91 subplot(222)
92 plot(dy/1e3,yc/1e3,'.'); ylabel('Y (km)');xlabel('\Delta y (km)')
93 title('\Delta y')
94
95 subplot(224)
96 plot(dh,zc,'.'); ylabel('Z (m)');xlabel('\Delta z (m)')
97 title('\Delta z')
98
99 fid=fopen('topo_sl.bin','w',ieee); fwrite(fid,hh,prec); fclose(fid);
100
101 % Initial Temp Profile
102 figure(2); clf
103
104 theta=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*tanh((YT-y0)/Lf)/2;
105
106 % Impose a strong initial and restoring mixed layer to Dml
107 i=1;
108 iml=1;
109
110 while ((ZT(1,1,iml)>-Dml)&(iml<length(ZT(1,1,:))))
111 iml=iml+1;
112 end
113
114 while (i<iml)
115 theta(:,:,i)=theta(:,:,iml);
116 i=i+1;
117 end
118
119 fprintf(' Z , Theta: j=1 , j=end\n');
120 fprintf(' %10.4f %10.4f %10.4f\n', ...
121 [ZT(1,1,1) theta(1,1,1) theta(1,end,1) ; ...
122 ZT(1,1,i) theta(1,1,i) theta(1,end,i); ...
123 ZT(1,1,end) theta(1,1,end) theta(1,end,end)]' ...
124 );
125
126 subplot(3,2,1)
127 [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
128 axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
129 colorbar
130 title('Potl Temp')
131 xlabel('y (km)');ylabel('z (m)')
132
133 subplot(3,2,2)
134 [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
135 colorbar
136 title('Potl Temp')
137 xlabel('y (km)');ylabel('z (m)')
138
139
140 subplot(3,2,3)
141 dens=theta*rho*tAlpha+rho;
142 [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
143 title('Density')
144 xlabel('y (km)');ylabel('z (m)')
145 axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
146
147 %Spice
148
149 x1=xc(round(length(xc)/4)); % centered location of Front #1
150 x2=xc(round(3*length(xc)/4)); % centered location of Front #2
151 spice=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;
152
153 i=1;
154 while (i<iml)
155 spice(:,:,i)=spice(:,:,iml);
156 i=i+1;
157 end
158
159 subplot(3,2,4)
160 [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
161 title('Density')
162 xlabel('y (km)');ylabel('z (m)')
163
164 subplot(3,2,5)
165 [h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(spice(:,1,:)));
166 title('Spice')
167 xlabel('x (km)');ylabel('z (m)')
168
169 subplot(3,2,6)
170 plot(YB(1,:)/1e3,hh(1,:))
171 title('topography')
172 xlabel('x (km)');ylabel('value')
173
174 %-- Perturb Initial temperature
175 pert=rand(size(theta(:,:,1)));
176 pert=1e-5*(pert-0.5);
177 for i=1:length(theta(1,1,:))
178 theta(:,:,i)=theta(:,:,i)+pert;
179 end
180
181 fid=fopen('thetaInitial.bin','w',ieee);
182 fwrite(fid,theta,prec); fclose(fid);
183
184 fid=fopen('spiceInitial.bin','w',ieee);
185 fwrite(fid,spice,prec); fclose(fid);
186

  ViewVC Help
Powered by ViewVC 1.1.22