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 |