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

Annotation 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 - (hide 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 dimitri 1.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