/[MITgcm]/MITgcm_contrib/high_res_cube/matlab/mk_bathy5.m
ViewVC logotype

Contents of /MITgcm_contrib/high_res_cube/matlab/mk_bathy5.m

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


Revision 1.1 - (show annotations) (download)
Thu Aug 10 22:31:39 2006 UTC (18 years, 11 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
Added files used to generate bathymetry and Levitus input files for cs510
integrations.

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3
4 % version 5 is a straight 20x20-km bin average
5 % same as version 3, but includes Black Sea
6
7 % read-in bin-averaged bathymetry
8 clear all, close all
9 nx=510;
10 fn='ETOPO2_510x6x510_ver3.bin';
11 fld=readbin(fn,[nx 6 nx],1);
12 fld=permute(fld,[1 3 2]);
13
14 % do not allow any numbers between 0 and 10
15 fld(find(fld>-1.5))=0;
16 fld(find(fld<0&fld>-10))=-10;
17
18 % dig connection to Black Sea
19 fld(389:397,505,1)=-10;
20 fld(402,510,1)=-10;
21 fld(33:36,84,3)=-10;
22
23 % make sure there are no lakes
24 % requires image processing toolbox
25 path('/home/dimitri/matlab/images/images',path);
26 tx=85;ty=85;nt=216;cx=510;cy=510;nx=510;genx;
27 te=fld;
28 te(find(fld>=0))=0;
29 te(find(fld<0))=1;
30 for tl=1:6
31 figure(1), clf
32 eval(['mkb' int2str(tl) '; caxis([0 1])'])
33 if tl==6
34 tmp2=imfill(logical(~tmp),[511 511]);
35 else
36 tmp2=imfill(logical(~tmp),[765 765]);
37 end
38 tmp(find(tmp2==0))=0;
39 figure(2), clf, orient tall, wysiwyg
40 subplot(311),mypcolor(te(:,:,tl)');colorbar
41 subplot(312),mypcolor(tmp(ix2,iy2)');colorbar
42 subplot(313),mypcolor(te(:,:,tl)'-tmp(ix2,iy2)');colorbar
43 te(:,:,tl)=tmp(ix2,iy2); pause(.01)
44 end
45
46 % save filled bathymetry
47 fld(find(te==0))=0;
48 fld=permute(fld,[1 3 2]);
49 fn='ETOPO2_510x6x510_ver5_filled.bin';
50 writebin(fn,fld,1);
51
52 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54
55 % look at bathymetry
56 clear all, close all
57 nx=510; genx
58 fn='ETOPO2_510x6x510_ver5_filled.bin';
59 fld=readbin(fn,[nx 6 nx],1);
60 fld=permute(fld,[1 3 2]);
61 tmp=nan*ones(nx*4,nx*3);
62 tmp(ix1,iy1)=fld(:,:,1);
63 tmp(ix2,iy1)=fld(:,:,2);
64 tmp(ix2,iy2)=fld(:,:,3);
65 tmp(ix3,iy2)=fld(:,:,4);
66 tmp(ix3,iy3)=fld(:,:,5);
67 tmp(ix4,iy3)=fld(:,:,6);
68 clf, mypcolor(tmp');
69 orient landscape, wysiwyg
70
71 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73
74 % fill in a hole that caused trouble
75 clear all, close all
76 tx=85;ty=85;nt=216;cx=510;cy=510;
77
78 te=readbin('ETOPO2_510x6x510_ver5_filled.bin',[510 6 510],1);
79 te=permute(te,[1 3 2]);
80
81 te0=te;
82 n=4;
83 ix=170; iy=340;
84 te(ix,iy,n) = ...
85 0.25*te0(ix,iy,n) + ...
86 0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
87 0.0625*(te0(ix-1,iy-1,n)+te0(ix-1,iy+1,n)+te0(ix+1,iy-1,n)+te0(ix+1,iy+1,n));
88 disp([te0(ix,iy,n) te(ix,iy,n)])
89 ix=170; iy=341;
90 te(ix,iy,n) = ...
91 0.25*te0(ix,iy,n) + ...
92 0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
93 0.0625*(te0(ix-1,iy-1,n)+te0(ix-1,iy+1,n)+te0(ix+1,iy-1,n)+te0(ix+1,iy+1,n));
94 disp([te0(ix,iy,n) te(ix,iy,n)])
95 ix=171; iy=341;
96 te(ix,iy,n) = ...
97 0.25*te0(ix,iy,n) + ...
98 0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
99 0.0625*(te0(ix-1,iy-1,n)+te0(ix-1,iy+1,n)+te0(ix+1,iy-1,n)+te0(ix+1,iy+1,n));
100 disp([te0(ix,iy,n) te(ix,iy,n)])
101
102 n=4; ix=169:172; iy=339:342;
103 subplot(211),mypcolor(te0(ix,iy,n)'),colorbar
104 title('fails at day 417')
105 subplot(212),mypcolor(te(ix,iy,n)'),colorbar
106 title('fails at day 522')
107 orient tall, wysiwyg
108
109 nz=1; anti_foo
110 writebin('ETOPO2_18360x85_ver5b_filled.bin',phi,1)
111
112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
114
115 % fill in all single-point wells
116 clear all, close all
117 tx=85;ty=85;nt=216;cx=510;cy=510;
118 nx=510; genx
119
120 phi=readbin('ETOPO2_18360x85_ver5b_filled.bin',[18360 85],1);
121 nz=1; foo; te2=te;
122
123 tmp=nan*ones(nx*4,nx*3);
124 tmp(ix1,iy1)=te(:,:,1);
125 tmp(ix2,iy1)=te(:,:,2);
126 tmp(ix2,iy2)=te(:,:,3);
127 tmp(ix3,iy2)=te(:,:,4);
128 tmp(ix3,iy3)=te(:,:,5);
129 tmp(ix4,iy3)=te(:,:,6);
130 clf, mypcolor(tmp');
131 orient landscape, wysiwyg
132
133 % cube face 1
134 mkb1; mkb7;
135 te2(:,:,1)=tmp2(ix2,iy2);
136
137 % cube face 2
138 mkb2; mkb7;
139 te2(:,:,2)=tmp2(ix2,iy2);
140
141 % cube face 3
142 mkb3; mkb7;
143 te2(:,:,3)=tmp2(ix2,iy2);
144
145 % cube face 4
146 mkb4; mkb7;
147 te2(:,:,4)=tmp2(ix2,iy2);
148
149 % cube face 5
150 mkb5; mkb7;
151 te2(:,:,5)=tmp2(ix2,iy2);
152
153 % cube face 6
154 mkb6; mkb7;
155 te2(:,:,6)=tmp2(ix2,iy2);
156
157 te=te2; nz=1; anti_foo
158 writebin('ETOPO2_18360x85_ver5c_filled.bin',phi,1)
159
160 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
162
163 % fill in all North-South two-point wells
164 clear all, close all
165 tx=85;ty=85;nt=216;cx=510;cy=510;
166 nx=510; genx
167
168 phi=readbin('ETOPO2_18360x85_ver5c_filled.bin',[18360 85],1);
169 nz=1; foo; te2=te;
170
171 tmp=nan*ones(nx*4,nx*3);
172 tmp(ix1,iy1)=te(:,:,1);
173 tmp(ix2,iy1)=te(:,:,2);
174 tmp(ix2,iy2)=te(:,:,3);
175 tmp(ix3,iy2)=te(:,:,4);
176 tmp(ix3,iy3)=te(:,:,5);
177 tmp(ix4,iy3)=te(:,:,6);
178 clf, mypcolor(tmp');
179 orient landscape, wysiwyg
180
181 % cube face 1
182 mkb1, mkb8
183 te2(:,:,1)=tmp2(ix2,iy2);
184
185 % cube face 2
186 mkb2, mkb8
187 te2(:,:,2)=tmp2(ix2,iy2);
188
189 % cube face 3
190 mkb3, mkb8
191 te2(:,:,3)=tmp2(ix2,iy2);
192
193 % cube face 4
194 mkb4, mkb8
195 te2(:,:,4)=tmp2(ix2,iy2);
196
197 % cube face 5
198 mkb5, mkb8
199 te2(:,:,5)=tmp2(ix2,iy2);
200
201 % cube face 6
202 mkb6, mkb8
203 te2(:,:,6)=tmp2(ix2,iy2);
204
205 % plot new bathymetry
206 tmp=nan*ones(nx*4,nx*3);
207 tmp(ix1,iy1)=te2(:,:,1);
208 tmp(ix2,iy1)=te2(:,:,2);
209 tmp(ix2,iy2)=te2(:,:,3);
210 tmp(ix3,iy2)=te2(:,:,4);
211 tmp(ix3,iy3)=te2(:,:,5);
212 tmp(ix4,iy3)=te2(:,:,6);
213 clf, mypcolor(tmp');
214 orient landscape, wysiwyg
215
216 te=te2; nz=1; anti_foo
217 writebin('ETOPO2_18360x85_ver5d_filled.bin',phi,1)
218
219
220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
222
223 % fill in all East-West two-point wells
224 clear all, close all
225 tx=85;ty=85;nt=216;cx=510;cy=510;
226 nx=510; genx
227
228 phi=readbin('ETOPO2_18360x85_ver5d_filled.bin',[18360 85],1);
229 nz=1; foo; te2=te;
230
231 tmp=nan*ones(nx*4,nx*3);
232 tmp(ix1,iy1)=te(:,:,1);
233 tmp(ix2,iy1)=te(:,:,2);
234 tmp(ix2,iy2)=te(:,:,3);
235 tmp(ix3,iy2)=te(:,:,4);
236 tmp(ix3,iy3)=te(:,:,5);
237 tmp(ix4,iy3)=te(:,:,6);
238 clf, mypcolor(tmp');
239 orient landscape, wysiwyg
240
241 % cube face 1
242 mkb1, mkb9
243 te2(:,:,1)=tmp2(ix2,iy2);
244
245 % cube face 2
246 mkb2, mkb9
247 te2(:,:,2)=tmp2(ix2,iy2);
248
249 % cube face 3
250 mkb3, mkb9
251 te2(:,:,3)=tmp2(ix2,iy2);
252
253 % cube face 4
254 mkb4, mkb9
255 te2(:,:,4)=tmp2(ix2,iy2);
256
257 % cube face 5
258 mkb5, mkb9
259 te2(:,:,5)=tmp2(ix2,iy2);
260
261 % cube face 6
262 mkb6, mkb9
263 te2(:,:,6)=tmp2(ix2,iy2);
264
265 % plot new bathymetry
266 tmp=nan*ones(nx*4,nx*3);
267 tmp(ix1,iy1)=te2(:,:,1);
268 tmp(ix2,iy1)=te2(:,:,2);
269 tmp(ix2,iy2)=te2(:,:,3);
270 tmp(ix3,iy2)=te2(:,:,4);
271 tmp(ix3,iy3)=te2(:,:,5);
272 tmp(ix4,iy3)=te2(:,:,6);
273 clf, mypcolor(tmp');
274 orient landscape, wysiwyg
275
276 te=te2; nz=1; anti_foo
277 writebin('ETOPO2_18360x85_ver5e_filled.bin',phi,1)
278
279
280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
282 % compare ver5 to ver5e
283
284 % fill in all East-West two-point wells
285 clear all, close all
286 tx=85;ty=85;nt=216;cx=510;cy=510;
287 fn='ETOPO2_510x6x510_ver5_filled.bin';
288 te1=readbin('ETOPO2_510x6x510_ver5_filled.bin',[cx 6 cy],1);
289 te1=permute(te1,[1 3 2]);
290 phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85],1);
291 nz=1; foo; te2=te;
292 clear c* i* j* n* p* te tep* tx ty
293 for f=1:6, clf
294 subplot(311), mypcolor(te1(:,:,f)');caxis([-10 0]),colorbar
295 subplot(312), mypcolor(te2(:,:,f)');caxis([-10 0]),colorbar
296 subplot(313), mypcolor(te2(:,:,f)'-te1(:,:,f)');colorbar
297 pause
298 end
299
300
301 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303
304 % write out as 205 tiles with blanks
305 clear all, close all
306 load blanklist.s205t_85x85
307
308 tx=85;ty=85;nt=216;cx=510;cy=510;nz=1;
309 phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85 nz],1);
310 foo;
311 anti_foo2
312 writebin('ETOPO2_17425x85_ver5e_filled.bin',phi)
313
314 nz=50;
315 phi=readbin('LEVS01_JAN_18360x85x50.bin',[18360 85 nz],1);
316 foo;
317 anti_foo2
318 writebin('LEVS01_JAN_17425x85x50.bin',phi)
319
320 phi=readbin('LEVT01_JAN_18360x85x50.bin',[18360 85 nz],1);
321 foo;
322 anti_foo2
323 writebin('LEVT01_JAN_17425x85x50.bin',phi)
324
325
326 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
327 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
328
329 % write out in 510x6x510 format
330 clear all, close all
331 tx=85;ty=85;nt=216;cx=510;cy=510;nz=1;
332
333 phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85 nz],1);
334 foo; te=permute(te,[1 3 2]);
335 writebin('ETOPO2_510x6x510_ver5e_filled.bin',te)
336
337 nz=50;
338 phi=readbin('LEVS01_JAN_18360x85x50.bin',[18360 85 nz],1);
339 foo; te=permute(te,[1 3 2 4]);
340 writebin('LEVS01_JAN_510x6x510x50.bin',te)
341
342 phi=readbin('LEVT01_JAN_18360x85x50.bin',[18360 85 nz],1);
343 foo; te=permute(te,[1 3 2 4]);
344 writebin('LEVT01_JAN_510x6x510x50.bin',te)
345
346
347 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
348 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
349
350 % compare Baltic bathymetry
351 clear all, close all
352 tx=85;ty=85;nt=216;cx=510;cy=510;nz=1;
353 phi=readbin('../output/data/grid/XC.data',[18360 85]); foo; xc=te(:,:,3);
354 phi=readbin('../output/data/grid/YC.data',[18360 85]); foo; yc=te(:,:,3);
355 phi=readbin('ETOPO2_510x6x510_ver5e_filled.bin',[510 6 510]);
356 phi=permute(phi,[1 3 2]); phi=phi(:,:,3);
357
358 load bathy_baltic_lon_lat_depth.txt
359 lon=bathy_baltic_lon_lat_depth(:,1);
360 lat=bathy_baltic_lon_lat_depth(:,2);
361 dpt=bathy_baltic_lon_lat_depth(:,3);
362 clear b* c* i* j* n* t*
363 pcolor(lon,lat,dpt)
364
365 DPT=0*xc; disp(length(lat))
366 for i=1:length(lat), mydisp(i)
367 ix=find(abs(xc-lon(i))<.01&abs(yc-lat(i))<.01);
368 DPT(ix)=dpt(i);
369 end
370 DPT(find(DPT<-1e30))=0;
371
372 ix=40:120; iy=155:225; cx=[-50 0];
373 clf,subplot(311),mypcolor(phi(ix,iy)'); caxis(cx), colorbar
374 subplot(312),mypcolor(DPT(ix,iy)'); caxis(cx), colorbar
375 subplot(313),mypcolor(phi(ix,iy)'-DPT(ix,iy)'); caxis([-1 1]*20), colorbar
376
377
378
379
380
381
382 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
383 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
384
385
386 % smooth a seamount that caused trouble
387 clear all, close all
388 tx=85;ty=85;nt=216;cx=510;cy=510;
389
390 phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85],1);
391 nz=1; foo;
392
393 mypcolor(te(:,:,4)')
394
395
396
397 nz=1; anti_foo
398 writebin('ETOPO2_18360x85_ver5f_filled.bin',phi,1)

  ViewVC Help
Powered by ViewVC 1.1.22