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

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