% version 3 is a straight 20x20-km bin average

% load etopo2 bathymetry
clear all, close all
nx=510;
LATC=readbin('LATC.bin',[nx 6 nx],1,'real*8');
LONC=readbin('LONC.bin',[nx 6 nx],1,'real*8');
ix=find(LONC<0); LONC(ix)=LONC(ix)+360;
f=netcdf('/host/nemo/tmp2/odap/dimitri/iter69/etopo2.nc');
lat=f{1}(:); lon=f{2}(:);
ix=find(lon<0); lon(ix)=lon(ix)+360;
[lon ix]=sort(lon); [lat iy]=sort(lat);
etopo2=f{3}(:); etopo2=etopo2(iy,ix);

% pad etopo2 bathymetry to the north
iy2=5400:-1:5350;
ix=1:10800; ix2=ix+5400;
ix=find(ix2>10800); ix2(ix)=ix2(ix)-10800;
etopo2=[etopo2; etopo2(iy2,ix2)];
for i=1:length(iy2)
  lat=[lat; max(lat)+lat(2)-lat(1)];
end
ix=1:10800; iy=5350:5452;

% filter and subsample to LONC, LATC
x=11; dy=111.195;
lon2=[lon-360; lon; lon+360];
TOPO=LONC;
disp(length(LONC(:)))
for i=1:length(LONC(:)), mydisp(i)
  iy=find(abs(lat-LATC(i))*dy<=x);
  dx=dy*cos(LATC(i)*pi/180);
  ix=find(abs(lon2-LONC(i))*dx<=x);
  ix=ix-10800;
  ix2=find(ix<1); ix(ix2)=ix(ix2)+10800;
  ix2=find(ix>10800); ix(ix2)=ix(ix2)-10800;
  topo=etopo2(iy,ix);
  TOPO(i)=mean(topo(:));
end
writebin('ETOPO2_510x6x510_ver3.bin',TOPO)

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

% fix 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;

% make sure there are no lakes
% requires image processing toolbox
path('/home/dimitri/matlab/images/images',path);
fld2=fld;
fld2(find(fld>=0))=0;
fld2(find(fld<0))=1;

tl=1;
tmp=fld2(:,:,tl);
clf,mypcolor(tmp');
cm=colormap; cm(1,:)=[0 0 0];
cm(64,:)=[1 1 1]; colormap(cm)
tmp(225,469,tl)=1;
clf,mypcolor(tmp');
tmp2=imfill(logical(~tmp),[1 1;476 410]);
clf,mypcolor(tmp2'); colorbar
tmp(find(tmp2==0))=0;
clf
subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
subplot(312),mypcolor(tmp'); colorbar
subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
fld2(:,:,tl)=tmp;

tl=2;
tmp=fld2(:,:,tl);
clf,mypcolor(tmp');
tmp2=imfill(logical(~tmp),[1 1]);
clf,mypcolor(tmp2'); colorbar
tmp(find(tmp2==0))=0;
clf
subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
subplot(312),mypcolor(tmp'); colorbar
subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
fld2(:,:,tl)=tmp;

tl=3;
tmp=fld2(:,:,tl);
clf,mypcolor(tmp');
tmp2=imfill(logical(~tmp),[1 500;1 190]);
clf,mypcolor(tmp2'); colorbar
tmp(find(tmp2==0))=0;
clf
subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
subplot(312),mypcolor(tmp'); colorbar
subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
fld2(:,:,tl)=tmp;

tl=4;
tmp=fld2(:,:,tl);
clf,mypcolor(tmp');
tmp2=imfill(logical(~tmp),[1 500;510 1;508 31]);
clf,mypcolor(tmp2'); colorbar
tmp(find(tmp2==0))=0;
clf
subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
subplot(312),mypcolor(tmp'); colorbar
subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
fld2(:,:,tl)=tmp;

tl=5;
tmp=fld2(:,:,tl);
clf,mypcolor(tmp');
tmp2=imfill(logical(~tmp),[1 1;1 510;510 510]);
clf,mypcolor(tmp2'); colorbar
tmp(find(tmp2==0))=0;
clf
subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
subplot(312),mypcolor(tmp'); colorbar
subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
fld2(:,:,tl)=tmp;

tl=6;
tmp=fld2(:,:,tl);
clf,mypcolor(tmp');
tmp2=imfill(logical(~tmp),[1 1]);
clf,mypcolor(tmp2'); colorbar
tmp(find(tmp2==0))=0;
clf
subplot(311),mypcolor(fld2(:,:,tl)'); colorbar
subplot(312),mypcolor(tmp'); colorbar
subplot(313),mypcolor(fld2(:,:,tl)'-tmp'); colorbar
fld2(:,:,tl)=tmp;

iz=find(fld2==0);
fld2=fld;
fld2(iz)=0;
fld2=permute(fld2,[1 3 2]);
fn='ETOPO2_510x6x510_ver3_filled.bin';
writebin(fn,fld2,1);

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

% look at bathymetry
clear all, close all
nx=510; genx
fn='ETOPO2_510x6x510_ver3_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

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

% convert to proper format for model input
clear all, close all
tx=85;ty=85;nt=216;cx=510;cy=510;

te=readbin('ETOPO2_510x6x510_ver3_filled.bin',[510 6 510],1);
te=permute(te,[1 3 2]);
nz=1; anti_foo
writebin('ETOPO2_18360x85_ver3_filled.bin',phi,1)


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

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

te=readbin('ETOPO2_510x6x510_ver3_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_ver3b_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_ver3b_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);

% 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_ver3c_filled.bin',phi,1)


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

% extract Baltic
clear all, close all
tx=85; ty=85; nt=216; cx=510; cy=510;
phi=readbin('ETOPO2_18360x85_ver3c_filled.bin',[18360 85],1);
nz=1; foo; bathy=te;
LATC=readbin('LATC.bin',[cx 6 cy],1,'real*8');
LONC=readbin('LONC.bin',[cx 6 cy],1,'real*8');
lat=permute(LATC,[1 3 2]);
lon=permute(LONC,[1 3 2]);
ix=20:120; iy=150:270;
lon=lon(ix,iy,3);
lat=lat(ix,iy,3);
bathy=bathy(ix,iy,3);
clear L* c* i* j* n* p* t*

pcolor(lon',lat',bathy'), shading flat
caxis([-200 0]), colorbar
save bathy


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

% 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_ver3c_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_ver3d_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_ver3d_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_ver3e_filled.bin',phi,1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare ver3 to ver3e

% 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_ver3_filled.bin',[18360 85],1);
nz=1; foo; te1=te;
phi=readbin('ETOPO2_18360x85_ver3e_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([-5 0]),colorbar
  subplot(312), mypcolor(te2(:,:,f)'),caxis([-5 0]),colorbar
  subplot(313), mypcolor(te2(:,:,f)'-te1(:,:,f)'),colorbar
  pause
end