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

Contents of /MITgcm_contrib/high_res_cube/matlab/mk_bathy3.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:06 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 % version 3 is a straight 20x20-km bin average
2
3 % load etopo2 bathymetry
4 clear all, close all
5 nx=510;
6 LATC=readbin('LATC.bin',[nx 6 nx],1,'real*8');
7 LONC=readbin('LONC.bin',[nx 6 nx],1,'real*8');
8 ix=find(LONC<0); LONC(ix)=LONC(ix)+360;
9 f=netcdf('/host/nemo/tmp2/odap/dimitri/iter69/etopo2.nc');
10 lat=f{1}(:); lon=f{2}(:);
11 ix=find(lon<0); lon(ix)=lon(ix)+360;
12 [lon ix]=sort(lon); [lat iy]=sort(lat);
13 etopo2=f{3}(:); etopo2=etopo2(iy,ix);
14
15 % pad etopo2 bathymetry to the north
16 iy2=5400:-1:5350;
17 ix=1:10800; ix2=ix+5400;
18 ix=find(ix2>10800); ix2(ix)=ix2(ix)-10800;
19 etopo2=[etopo2; etopo2(iy2,ix2)];
20 for i=1:length(iy2)
21 lat=[lat; max(lat)+lat(2)-lat(1)];
22 end
23 ix=1:10800; iy=5350:5452;
24
25 % filter and subsample to LONC, LATC
26 x=11; dy=111.195;
27 lon2=[lon-360; lon; lon+360];
28 TOPO=LONC;
29 disp(length(LONC(:)))
30 for i=1:length(LONC(:)), mydisp(i)
31 iy=find(abs(lat-LATC(i))*dy<=x);
32 dx=dy*cos(LATC(i)*pi/180);
33 ix=find(abs(lon2-LONC(i))*dx<=x);
34 ix=ix-10800;
35 ix2=find(ix<1); ix(ix2)=ix(ix2)+10800;
36 ix2=find(ix>10800); ix(ix2)=ix(ix2)-10800;
37 topo=etopo2(iy,ix);
38 TOPO(i)=mean(topo(:));
39 end
40 writebin('ETOPO2_510x6x510_ver3.bin',TOPO)
41
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44
45 % fix bathymetry
46 clear all, close all
47 nx=510;
48 fn='ETOPO2_510x6x510_ver3.bin';
49 fld=readbin(fn,[nx 6 nx],1);
50 fld=permute(fld,[1 3 2]);
51
52 % do not allow any numbers between 0 and 10
53 fld(find(fld>-1.5))=0;
54 fld(find(fld<0&fld>-10))=-10;
55
56 % make sure there are no lakes
57 % requires image processing toolbox
58 path('/home/dimitri/matlab/images/images',path);
59 fld2=fld;
60 fld2(find(fld>=0))=0;
61 fld2(find(fld<0))=1;
62
63 tl=1;
64 tmp=fld2(:,:,tl);
65 clf,mypcolor(tmp');
66 cm=colormap; cm(1,:)=[0 0 0];
67 cm(64,:)=[1 1 1]; colormap(cm)
68 tmp(225,469,tl)=1;
69 clf,mypcolor(tmp');
70 tmp2=imfill(logical(~tmp),[1 1;476 410]);
71 clf,mypcolor(tmp2'); colorbar
72 tmp(find(tmp2==0))=0;
73 clf
74 subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
75 subplot(312),mypcolor(tmp'); colorbar
76 subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
77 fld2(:,:,tl)=tmp;
78
79 tl=2;
80 tmp=fld2(:,:,tl);
81 clf,mypcolor(tmp');
82 tmp2=imfill(logical(~tmp),[1 1]);
83 clf,mypcolor(tmp2'); colorbar
84 tmp(find(tmp2==0))=0;
85 clf
86 subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
87 subplot(312),mypcolor(tmp'); colorbar
88 subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
89 fld2(:,:,tl)=tmp;
90
91 tl=3;
92 tmp=fld2(:,:,tl);
93 clf,mypcolor(tmp');
94 tmp2=imfill(logical(~tmp),[1 500;1 190]);
95 clf,mypcolor(tmp2'); colorbar
96 tmp(find(tmp2==0))=0;
97 clf
98 subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
99 subplot(312),mypcolor(tmp'); colorbar
100 subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
101 fld2(:,:,tl)=tmp;
102
103 tl=4;
104 tmp=fld2(:,:,tl);
105 clf,mypcolor(tmp');
106 tmp2=imfill(logical(~tmp),[1 500;510 1;508 31]);
107 clf,mypcolor(tmp2'); colorbar
108 tmp(find(tmp2==0))=0;
109 clf
110 subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
111 subplot(312),mypcolor(tmp'); colorbar
112 subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
113 fld2(:,:,tl)=tmp;
114
115 tl=5;
116 tmp=fld2(:,:,tl);
117 clf,mypcolor(tmp');
118 tmp2=imfill(logical(~tmp),[1 1;1 510;510 510]);
119 clf,mypcolor(tmp2'); colorbar
120 tmp(find(tmp2==0))=0;
121 clf
122 subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
123 subplot(312),mypcolor(tmp'); colorbar
124 subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
125 fld2(:,:,tl)=tmp;
126
127 tl=6;
128 tmp=fld2(:,:,tl);
129 clf,mypcolor(tmp');
130 tmp2=imfill(logical(~tmp),[1 1]);
131 clf,mypcolor(tmp2'); colorbar
132 tmp(find(tmp2==0))=0;
133 clf
134 subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
135 subplot(312),mypcolor(tmp'); colorbar
136 subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
137 fld2(:,:,tl)=tmp;
138
139 iz=find(fld2==0);
140 fld2=fld;
141 fld2(iz)=0;
142 fld2=permute(fld2,[1 3 2]);
143 fn='ETOPO2_510x6x510_ver3_filled.bin';
144 writebin(fn,fld2,1);
145
146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
148
149 % look at bathymetry
150 clear all, close all
151 nx=510; genx
152 fn='ETOPO2_510x6x510_ver3_filled.bin';
153 fld=readbin(fn,[nx 6 nx],1);
154 fld=permute(fld,[1 3 2]);
155 tmp=nan*ones(nx*4,nx*3);
156 tmp(ix1,iy1)=fld(:,:,1);
157 tmp(ix2,iy1)=fld(:,:,2);
158 tmp(ix2,iy2)=fld(:,:,3);
159 tmp(ix3,iy2)=fld(:,:,4);
160 tmp(ix3,iy3)=fld(:,:,5);
161 tmp(ix4,iy3)=fld(:,:,6);
162 clf, mypcolor(tmp');
163 orient landscape, wysiwyg
164
165 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
167
168 % convert to proper format for model input
169 clear all, close all
170 tx=85;ty=85;nt=216;cx=510;cy=510;
171
172 te=readbin('ETOPO2_510x6x510_ver3_filled.bin',[510 6 510],1);
173 te=permute(te,[1 3 2]);
174 nz=1; anti_foo
175 writebin('ETOPO2_18360x85_ver3_filled.bin',phi,1)
176
177
178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180
181 % fill in a hole
182 clear all, close all
183 tx=85;ty=85;nt=216;cx=510;cy=510;
184
185 te=readbin('ETOPO2_510x6x510_ver3_filled.bin',[510 6 510],1);
186 te=permute(te,[1 3 2]);
187
188 te0=te;
189 n=4;
190 ix=170; iy=340;
191 te(ix,iy,n) = ...
192 0.25*te0(ix,iy,n) + ...
193 0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
194 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));
195 disp([te0(ix,iy,n) te(ix,iy,n)])
196 ix=170; iy=341;
197 te(ix,iy,n) = ...
198 0.25*te0(ix,iy,n) + ...
199 0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
200 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));
201 disp([te0(ix,iy,n) te(ix,iy,n)])
202 ix=171; iy=341;
203 te(ix,iy,n) = ...
204 0.25*te0(ix,iy,n) + ...
205 0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
206 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));
207 disp([te0(ix,iy,n) te(ix,iy,n)])
208
209 n=4; ix=169:172; iy=339:342;
210 subplot(211),mypcolor(te0(ix,iy,n)'),colorbar
211 title('fails at day 417')
212 subplot(212),mypcolor(te(ix,iy,n)'),colorbar
213 title('fails at day 522')
214 orient tall, wysiwyg
215
216 nz=1; anti_foo
217 writebin('ETOPO2_18360x85_ver3b_filled.bin',phi,1)
218
219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221
222 % fill in all single-point wells
223 clear all, close all
224 tx=85;ty=85;nt=216;cx=510;cy=510;
225 nx=510; genx
226
227 phi=readbin('ETOPO2_18360x85_ver3b_filled.bin',[18360 85],1);
228 nz=1; foo; te2=te;
229
230 tmp=nan*ones(nx*4,nx*3);
231 tmp(ix1,iy1)=te(:,:,1);
232 tmp(ix2,iy1)=te(:,:,2);
233 tmp(ix2,iy2)=te(:,:,3);
234 tmp(ix3,iy2)=te(:,:,4);
235 tmp(ix3,iy3)=te(:,:,5);
236 tmp(ix4,iy3)=te(:,:,6);
237 clf, mypcolor(tmp');
238 orient landscape, wysiwyg
239
240 % cube face 1
241 mkb1; mkb7;
242 te2(:,:,1)=tmp2(ix2,iy2);
243
244 % cube face 2
245 mkb2; mkb7;
246 te2(:,:,2)=tmp2(ix2,iy2);
247
248 % cube face 3
249 mkb3; mkb7;
250 te2(:,:,3)=tmp2(ix2,iy2);
251
252 % cube face 4
253 mkb4; mkb7;
254 te2(:,:,4)=tmp2(ix2,iy2);
255
256 % cube face 5
257 mkb5; mkb7;
258 te2(:,:,5)=tmp2(ix2,iy2);
259
260 % cube face 6
261 mkb6; mkb7;
262 te2(:,:,6)=tmp2(ix2,iy2);
263
264 % plot new bathymetry
265 tmp=nan*ones(nx*4,nx*3);
266 tmp(ix1,iy1)=te2(:,:,1);
267 tmp(ix2,iy1)=te2(:,:,2);
268 tmp(ix2,iy2)=te2(:,:,3);
269 tmp(ix3,iy2)=te2(:,:,4);
270 tmp(ix3,iy3)=te2(:,:,5);
271 tmp(ix4,iy3)=te2(:,:,6);
272 clf, mypcolor(tmp');
273 orient landscape, wysiwyg
274
275 te=te2; nz=1; anti_foo
276 writebin('ETOPO2_18360x85_ver3c_filled.bin',phi,1)
277
278
279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281
282 % extract Baltic
283 clear all, close all
284 tx=85; ty=85; nt=216; cx=510; cy=510;
285 phi=readbin('ETOPO2_18360x85_ver3c_filled.bin',[18360 85],1);
286 nz=1; foo; bathy=te;
287 LATC=readbin('LATC.bin',[cx 6 cy],1,'real*8');
288 LONC=readbin('LONC.bin',[cx 6 cy],1,'real*8');
289 lat=permute(LATC,[1 3 2]);
290 lon=permute(LONC,[1 3 2]);
291 ix=20:120; iy=150:270;
292 lon=lon(ix,iy,3);
293 lat=lat(ix,iy,3);
294 bathy=bathy(ix,iy,3);
295 clear L* c* i* j* n* p* t*
296
297 pcolor(lon',lat',bathy'), shading flat
298 caxis([-200 0]), colorbar
299 save bathy
300
301
302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
304
305 % fill in all North-South two-point wells
306 clear all, close all
307 tx=85;ty=85;nt=216;cx=510;cy=510;
308 nx=510; genx
309
310 phi=readbin('ETOPO2_18360x85_ver3c_filled.bin',[18360 85],1);
311 nz=1; foo; te2=te;
312
313 tmp=nan*ones(nx*4,nx*3);
314 tmp(ix1,iy1)=te(:,:,1);
315 tmp(ix2,iy1)=te(:,:,2);
316 tmp(ix2,iy2)=te(:,:,3);
317 tmp(ix3,iy2)=te(:,:,4);
318 tmp(ix3,iy3)=te(:,:,5);
319 tmp(ix4,iy3)=te(:,:,6);
320 clf, mypcolor(tmp');
321 orient landscape, wysiwyg
322
323 % cube face 1
324 mkb1, mkb8
325 te2(:,:,1)=tmp2(ix2,iy2);
326
327 % cube face 2
328 mkb2, mkb8
329 te2(:,:,2)=tmp2(ix2,iy2);
330
331 % cube face 3
332 mkb3, mkb8
333 te2(:,:,3)=tmp2(ix2,iy2);
334
335 % cube face 4
336 mkb4, mkb8
337 te2(:,:,4)=tmp2(ix2,iy2);
338
339 % cube face 5
340 mkb5, mkb8
341 te2(:,:,5)=tmp2(ix2,iy2);
342
343 % cube face 6
344 mkb6, mkb8
345 te2(:,:,6)=tmp2(ix2,iy2);
346
347 % plot new bathymetry
348 tmp=nan*ones(nx*4,nx*3);
349 tmp(ix1,iy1)=te2(:,:,1);
350 tmp(ix2,iy1)=te2(:,:,2);
351 tmp(ix2,iy2)=te2(:,:,3);
352 tmp(ix3,iy2)=te2(:,:,4);
353 tmp(ix3,iy3)=te2(:,:,5);
354 tmp(ix4,iy3)=te2(:,:,6);
355 clf, mypcolor(tmp');
356 orient landscape, wysiwyg
357
358 te=te2; nz=1; anti_foo
359 writebin('ETOPO2_18360x85_ver3d_filled.bin',phi,1)
360
361
362 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
364
365 % fill in all East-West two-point wells
366 clear all, close all
367 tx=85;ty=85;nt=216;cx=510;cy=510;
368 nx=510; genx
369
370 phi=readbin('ETOPO2_18360x85_ver3d_filled.bin',[18360 85],1);
371 nz=1; foo; te2=te;
372
373 tmp=nan*ones(nx*4,nx*3);
374 tmp(ix1,iy1)=te(:,:,1);
375 tmp(ix2,iy1)=te(:,:,2);
376 tmp(ix2,iy2)=te(:,:,3);
377 tmp(ix3,iy2)=te(:,:,4);
378 tmp(ix3,iy3)=te(:,:,5);
379 tmp(ix4,iy3)=te(:,:,6);
380 clf, mypcolor(tmp');
381 orient landscape, wysiwyg
382
383 % cube face 1
384 mkb1, mkb9
385 te2(:,:,1)=tmp2(ix2,iy2);
386
387 % cube face 2
388 mkb2, mkb9
389 te2(:,:,2)=tmp2(ix2,iy2);
390
391 % cube face 3
392 mkb3, mkb9
393 te2(:,:,3)=tmp2(ix2,iy2);
394
395 % cube face 4
396 mkb4, mkb9
397 te2(:,:,4)=tmp2(ix2,iy2);
398
399 % cube face 5
400 mkb5, mkb9
401 te2(:,:,5)=tmp2(ix2,iy2);
402
403 % cube face 6
404 mkb6, mkb9
405 te2(:,:,6)=tmp2(ix2,iy2);
406
407 % plot new bathymetry
408 tmp=nan*ones(nx*4,nx*3);
409 tmp(ix1,iy1)=te2(:,:,1);
410 tmp(ix2,iy1)=te2(:,:,2);
411 tmp(ix2,iy2)=te2(:,:,3);
412 tmp(ix3,iy2)=te2(:,:,4);
413 tmp(ix3,iy3)=te2(:,:,5);
414 tmp(ix4,iy3)=te2(:,:,6);
415 clf, mypcolor(tmp');
416 orient landscape, wysiwyg
417
418 te=te2; nz=1; anti_foo
419 writebin('ETOPO2_18360x85_ver3e_filled.bin',phi,1)
420
421
422 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
424 % compare ver3 to ver3e
425
426 % fill in all East-West two-point wells
427 clear all, close all
428 tx=85;ty=85;nt=216;cx=510;cy=510;
429 nx=510; genx
430 phi=readbin('ETOPO2_18360x85_ver3_filled.bin',[18360 85],1);
431 nz=1; foo; te1=te;
432 phi=readbin('ETOPO2_18360x85_ver3e_filled.bin',[18360 85],1);
433 nz=1; foo; te2=te;
434 clear c* i* j* n* p* te tep* tx ty
435 for f=1:6, clf
436 subplot(311), mypcolor(te1(:,:,f)'),caxis([-5 0]),colorbar
437 subplot(312), mypcolor(te2(:,:,f)'),caxis([-5 0]),colorbar
438 subplot(313), mypcolor(te2(:,:,f)'-te1(:,:,f)'),colorbar
439 pause
440 end

  ViewVC Help
Powered by ViewVC 1.1.22