%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% version 5 is a straight 20x20-km bin average
% same as version 3, but includes Black Sea

% read-in bin-averaged bathymetry
clear all, close all
nx=510;
fn='ETOPO2_510x6x510_ver3.bin';
fld=readbin(fn,[nx 6 nx],1);
fld=permute(fld,[1 3 2]);

% do not allow any numbers between 0 and 10
fld(find(fld>-1.5))=0;
fld(find(fld<0&fld>-10))=-10;

% dig connection to Black Sea
fld(389:397,505,1)=-10;
fld(402,510,1)=-10;
fld(33:36,84,3)=-10;

% make sure there are no lakes
% requires image processing toolbox
path('/home/dimitri/matlab/images/images',path);
tx=85;ty=85;nt=216;cx=510;cy=510;nx=510;genx;
te=fld;
te(find(fld>=0))=0;
te(find(fld<0))=1;
for tl=1:6
  figure(1), clf
  eval(['mkb' int2str(tl) '; caxis([0 1])'])
  if tl==6
    tmp2=imfill(logical(~tmp),[511 511]);
  else
    tmp2=imfill(logical(~tmp),[765 765]);
  end
  tmp(find(tmp2==0))=0;
  figure(2), clf, orient tall, wysiwyg
  subplot(311),mypcolor(te(:,:,tl)');colorbar
  subplot(312),mypcolor(tmp(ix2,iy2)');colorbar
  subplot(313),mypcolor(te(:,:,tl)'-tmp(ix2,iy2)');colorbar
  te(:,:,tl)=tmp(ix2,iy2); pause(.01)
end

% save filled bathymetry
fld(find(te==0))=0;
fld=permute(fld,[1 3 2]);
fn='ETOPO2_510x6x510_ver5_filled.bin';
writebin(fn,fld,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% look at bathymetry
clear all, close all
nx=510; genx
fn='ETOPO2_510x6x510_ver5_filled.bin';
fld=readbin(fn,[nx 6 nx],1);
fld=permute(fld,[1 3 2]);
tmp=nan*ones(nx*4,nx*3);
tmp(ix1,iy1)=fld(:,:,1);
tmp(ix2,iy1)=fld(:,:,2);
tmp(ix2,iy2)=fld(:,:,3);
tmp(ix3,iy2)=fld(:,:,4);
tmp(ix3,iy3)=fld(:,:,5);
tmp(ix4,iy3)=fld(:,:,6);
clf, mypcolor(tmp');
orient landscape, wysiwyg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fill in a hole that caused trouble
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;

te=readbin('ETOPO2_510x6x510_ver5_filled.bin',[510 6 510],1);
te=permute(te,[1 3 2]);

te0=te;
n=4;
ix=170; iy=340;
te(ix,iy,n) = ...
  0.25*te0(ix,iy,n) + ...
  0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
  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));
disp([te0(ix,iy,n) te(ix,iy,n)])
ix=170; iy=341;
te(ix,iy,n) = ...
  0.25*te0(ix,iy,n) + ...
  0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
  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));
disp([te0(ix,iy,n) te(ix,iy,n)])
ix=171; iy=341;
te(ix,iy,n) = ...
  0.25*te0(ix,iy,n) + ...
  0.125*(te0(ix-1,iy,n)+te0(ix+1,iy,n)+te0(ix,iy-1,n)+te0(ix,iy+1,n)) + ...
  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));
disp([te0(ix,iy,n) te(ix,iy,n)])

n=4; ix=169:172; iy=339:342;
subplot(211),mypcolor(te0(ix,iy,n)'),colorbar
title('fails at day 417')
subplot(212),mypcolor(te(ix,iy,n)'),colorbar
title('fails at day 522')
orient tall, wysiwyg

nz=1; anti_foo
writebin('ETOPO2_18360x85_ver5b_filled.bin',phi,1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fill in all single-point wells
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;
nx=510; genx

phi=readbin('ETOPO2_18360x85_ver5b_filled.bin',[18360 85],1);
nz=1; foo; te2=te;

tmp=nan*ones(nx*4,nx*3);
tmp(ix1,iy1)=te(:,:,1);
tmp(ix2,iy1)=te(:,:,2);
tmp(ix2,iy2)=te(:,:,3);
tmp(ix3,iy2)=te(:,:,4);
tmp(ix3,iy3)=te(:,:,5);
tmp(ix4,iy3)=te(:,:,6);
clf, mypcolor(tmp');
orient landscape, wysiwyg

% cube face 1
mkb1; mkb7;
te2(:,:,1)=tmp2(ix2,iy2);

% cube face 2
mkb2; mkb7;
te2(:,:,2)=tmp2(ix2,iy2);

% cube face 3
mkb3; mkb7;
te2(:,:,3)=tmp2(ix2,iy2);

% cube face 4
mkb4; mkb7;
te2(:,:,4)=tmp2(ix2,iy2);

% cube face 5
mkb5; mkb7;
te2(:,:,5)=tmp2(ix2,iy2);

% cube face 6
mkb6; mkb7;
te2(:,:,6)=tmp2(ix2,iy2);

te=te2; nz=1; anti_foo
writebin('ETOPO2_18360x85_ver5c_filled.bin',phi,1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fill in all North-South two-point wells
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;
nx=510; genx

phi=readbin('ETOPO2_18360x85_ver5c_filled.bin',[18360 85],1);
nz=1; foo; te2=te;

tmp=nan*ones(nx*4,nx*3);
tmp(ix1,iy1)=te(:,:,1);
tmp(ix2,iy1)=te(:,:,2);
tmp(ix2,iy2)=te(:,:,3);
tmp(ix3,iy2)=te(:,:,4);
tmp(ix3,iy3)=te(:,:,5);
tmp(ix4,iy3)=te(:,:,6);
clf, mypcolor(tmp');
orient landscape, wysiwyg

% cube face 1
mkb1, mkb8
te2(:,:,1)=tmp2(ix2,iy2);

% cube face 2
mkb2, mkb8
te2(:,:,2)=tmp2(ix2,iy2);

% cube face 3
mkb3, mkb8
te2(:,:,3)=tmp2(ix2,iy2);

% cube face 4
mkb4, mkb8
te2(:,:,4)=tmp2(ix2,iy2);

% cube face 5
mkb5, mkb8
te2(:,:,5)=tmp2(ix2,iy2);

% cube face 6
mkb6, mkb8
te2(:,:,6)=tmp2(ix2,iy2);

% plot new bathymetry
tmp=nan*ones(nx*4,nx*3);
tmp(ix1,iy1)=te2(:,:,1);
tmp(ix2,iy1)=te2(:,:,2);
tmp(ix2,iy2)=te2(:,:,3);
tmp(ix3,iy2)=te2(:,:,4);
tmp(ix3,iy3)=te2(:,:,5);
tmp(ix4,iy3)=te2(:,:,6);
clf, mypcolor(tmp');
orient landscape, wysiwyg

te=te2; nz=1; anti_foo
writebin('ETOPO2_18360x85_ver5d_filled.bin',phi,1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fill in all East-West two-point wells
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;
nx=510; genx

phi=readbin('ETOPO2_18360x85_ver5d_filled.bin',[18360 85],1);
nz=1; foo; te2=te;

tmp=nan*ones(nx*4,nx*3);
tmp(ix1,iy1)=te(:,:,1);
tmp(ix2,iy1)=te(:,:,2);
tmp(ix2,iy2)=te(:,:,3);
tmp(ix3,iy2)=te(:,:,4);
tmp(ix3,iy3)=te(:,:,5);
tmp(ix4,iy3)=te(:,:,6);
clf, mypcolor(tmp');
orient landscape, wysiwyg

% cube face 1
mkb1, mkb9
te2(:,:,1)=tmp2(ix2,iy2);

% cube face 2
mkb2, mkb9
te2(:,:,2)=tmp2(ix2,iy2);

% cube face 3
mkb3, mkb9
te2(:,:,3)=tmp2(ix2,iy2);

% cube face 4
mkb4, mkb9
te2(:,:,4)=tmp2(ix2,iy2);

% cube face 5
mkb5, mkb9
te2(:,:,5)=tmp2(ix2,iy2);

% cube face 6
mkb6, mkb9
te2(:,:,6)=tmp2(ix2,iy2);

% plot new bathymetry
tmp=nan*ones(nx*4,nx*3);
tmp(ix1,iy1)=te2(:,:,1);
tmp(ix2,iy1)=te2(:,:,2);
tmp(ix2,iy2)=te2(:,:,3);
tmp(ix3,iy2)=te2(:,:,4);
tmp(ix3,iy3)=te2(:,:,5);
tmp(ix4,iy3)=te2(:,:,6);
clf, mypcolor(tmp');
orient landscape, wysiwyg

te=te2; nz=1; anti_foo
writebin('ETOPO2_18360x85_ver5e_filled.bin',phi,1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare ver5 to ver5e

% fill in all East-West two-point wells
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;
fn='ETOPO2_510x6x510_ver5_filled.bin';
te1=readbin('ETOPO2_510x6x510_ver5_filled.bin',[cx 6 cy],1);
te1=permute(te1,[1 3 2]);
phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85],1);
nz=1; foo; te2=te;
clear c* i* j* n* p* te tep* tx ty
for f=1:6, clf
  subplot(311), mypcolor(te1(:,:,f)');caxis([-10 0]),colorbar
  subplot(312), mypcolor(te2(:,:,f)');caxis([-10 0]),colorbar
  subplot(313), mypcolor(te2(:,:,f)'-te1(:,:,f)');colorbar
  pause
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% write out as 205 tiles with blanks
clear all, close all
load blanklist.s205t_85x85

tx=85;ty=85;nt=216;cx=510;cy=510;nz=1;
phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85 nz],1);
foo;
anti_foo2
writebin('ETOPO2_17425x85_ver5e_filled.bin',phi)

nz=50;
phi=readbin('LEVS01_JAN_18360x85x50.bin',[18360 85 nz],1);
foo;
anti_foo2
writebin('LEVS01_JAN_17425x85x50.bin',phi)

phi=readbin('LEVT01_JAN_18360x85x50.bin',[18360 85 nz],1);
foo;
anti_foo2
writebin('LEVT01_JAN_17425x85x50.bin',phi)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% write out in 510x6x510 format
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;nz=1;

phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85 nz],1);
foo; te=permute(te,[1 3 2]);
writebin('ETOPO2_510x6x510_ver5e_filled.bin',te)

nz=50;
phi=readbin('LEVS01_JAN_18360x85x50.bin',[18360 85 nz],1);
foo; te=permute(te,[1 3 2 4]);
writebin('LEVS01_JAN_510x6x510x50.bin',te)

phi=readbin('LEVT01_JAN_18360x85x50.bin',[18360 85 nz],1);
foo; te=permute(te,[1 3 2 4]);
writebin('LEVT01_JAN_510x6x510x50.bin',te)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compare Baltic bathymetry
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;nz=1;
phi=readbin('../output/data/grid/XC.data',[18360 85]); foo; xc=te(:,:,3);
phi=readbin('../output/data/grid/YC.data',[18360 85]); foo; yc=te(:,:,3);
phi=readbin('ETOPO2_510x6x510_ver5e_filled.bin',[510 6 510]);
phi=permute(phi,[1 3 2]); phi=phi(:,:,3);

load bathy_baltic_lon_lat_depth.txt
lon=bathy_baltic_lon_lat_depth(:,1);
lat=bathy_baltic_lon_lat_depth(:,2);
dpt=bathy_baltic_lon_lat_depth(:,3);
clear b* c* i* j* n* t*
pcolor(lon,lat,dpt)

DPT=0*xc; disp(length(lat))
for i=1:length(lat), mydisp(i)
  ix=find(abs(xc-lon(i))<.01&abs(yc-lat(i))<.01);
  DPT(ix)=dpt(i);
end
DPT(find(DPT<-1e30))=0;

ix=40:120; iy=155:225; cx=[-50 0];
clf,subplot(311),mypcolor(phi(ix,iy)'); caxis(cx), colorbar
subplot(312),mypcolor(DPT(ix,iy)'); caxis(cx), colorbar
subplot(313),mypcolor(phi(ix,iy)'-DPT(ix,iy)'); caxis([-1 1]*20), colorbar






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% smooth a seamount that caused trouble
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;

phi=readbin('ETOPO2_18360x85_ver5e_filled.bin',[18360 85],1);
nz=1; foo;

mypcolor(te(:,:,4)')



nz=1; anti_foo
writebin('ETOPO2_18360x85_ver5f_filled.bin',phi,1)