% This is a script that executes a series of commands for
% extracting GEBCO data, interpolating, editing and then
% generating a mask-file.
%
% Uses data files: HGRID.mat VGRID.mat
% Creates data files:  bathymetry.bin
% Creates arrays: n/a
% Uses m-scripts: setgrid grid_sphere finegrid gen_hxhy xyrecur_* edtopo 
%
% Created  08/15/99 by adcroft@mit.edu
% Modified 11/09/99 by adcroft@mit.edu
% Maintained by adcroft@mit.edu, abiastoch@ucsd.edu
% Modifield to read etopo2 01/10/04 by mlosch@awi-bremerhaven.de

%clear all
format compact
more off

skipfv=0

load HGRID.mat nxc nyc lon_lo lon_hi lat_lo lat_hi periodic

if ~skipfv

load HFGRID xcf
% Get the GEBCO data on a conveniently sized grid
disp('Extracting S2004 data ... ');
[He2,xe2,ye2]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi);
i0=find(xe2>=360+s_lon);
xe2(i0) = xe2(i0)-360;
if xe2(end)==xe2(1); xe2(end)=xe2(1)+360; end
% make the domain bigger so that it covers everything
[nxe,nye]=size(He2);
imax = min(find(xe2+360>max(xcf(:))));
imin = max(find(xe2-360<min(xcf(:))));
He2 = [He2(imin+1:end-1,:); He2; He2(2:imax,:)];
xe2 = [xe2(imin+1:end-1)-360 xe2 xe2(2:imax)+360];

disp(sprintf('S2004 dimensions are %i x %i',size(He2)));
save S2004.mat He2 xe2 ye2

disp('Extracting BRIOS data ... ');
[Hb2,Zb2,xb2,yb2]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi);
i0 = find(xb2>=360+s_lon);
xb2 = [xb2(i0)-360; xb2];
Hb2 = [Hb2(i0,:); Hb2];
Zb2 = [Zb2(i0,:); Zb2];

disp(sprintf('BRIOS dimensions are %i x %i',size(Hb2)));
save BRIOS.mat Hb2 Zb2 xb2 yb2

% Interpolate to fine grid
load HFGRID xcf ycf nfinepow
disp('Interpolating to fine grid ... ');
Hfine = 0*xcf;
Bfine = 0*xcf;
Zfine = 0*xcf;
% this may be replace with something more sophisticated (and expensive)
imeth = 'nearest'; 
imeth = 'linear'; 
disp(sprintf('interpolation method: %s',imeth))
largeproblem = 1;
if largeproblem
  % divide and conquer
  [nxt nyt]=size(xcf);
  nt = 2.^max(nfinepow-2,0);
  nxt = nxt./nt;
  nyt = nyt./nt;
  for j=1:nt
    jy = nyt*(j-1)+[1:nyt];
    for i=1:nt
      ix = nxt*(i-1)+[1:nxt];
      disp(sprintf('tile number = (%u,%u)',i,j));
      Hfine(ix,jy)=interp2(xe2(:),ye2(:)',He2',xcf(ix,jy),ycf(ix,jy),imeth);
      Bfine(ix,jy)=interp2(xb2(:),[yb2;-yb2(end:-1:1)]', ...
		       [Hb2 zeros(size(Hb2))]',xcf(ix,jy),ycf(ix,jy),imeth);
      Zfine(ix,jy)=interp2(xb2(:),[yb2;-yb2(end:-1:1)]', ...
		       [Zb2 zeros(size(Zb2))]', xcf(ix,jy),ycf(ix,jy),imeth);
    end
  end
else
  % do everything in one go
  Hfine=interp2(xe2(:),ye2(:)',He2',xcf,ycf,imeth);
  Bfine=interp2(xb2(:),[yb2;-yb2(end:-1:1)]',[Hb2 zeros(size(Hb2))]', ...
		xcf,ycf,imeth);
  Zfine=interp2(xb2(:),[yb2;-yb2(end:-1:1)]',[Zb2 zeros(size(Zb2))]', ...
		xcf,ycf,'nearest');
end

Bfine(isnan(Bfine)) = 0;
Zfine(isnan(Zfine)) = 0;
% add bathymetry under shelfice
msk0 = Hfine;
msk0(find(msk0>0))=1;
sizmsk = Zfine;
sihmsk = Bfine;
sihmsk(find(sihmsk~=0))=1;
sizmsk(find(sizmsk~=0))=1;
%sizmsk = sizmsk.*sihmsk;
% introduce wall at the south (affects only Ross Sea)
Zfine = Zfine.*sizmsk;
% under ice bottom topography
i1 = find(sizmsk == 1);
Hfine(i1) = Bfine(i1);
%% Handarbeit:
%x0=
%y0=
ierr = find(Zfine(i1)<Hfine(i1));

clear Bfine Hb2 xe2 ye2

%load HFINE
disp(sprintf('Fine grid dimensions are %i x %i x %i',size(Hfine)));
save HFINE.mat Hfine xcf ycf
save ZFINE.mat Zfine xcf ycf
clear He2 xe2 ye2 xcf ycf

disp('Adding halo region ...');
load HFGRID nfinepow
Hfine=exch_llc(Hfine,2^(nfinepow));
disp('Limiting to ocean (h<=0) ...');
Hfine=min(Hfine,0);
disp(sprintf('Fine grid+halo dimensions are %i x %i x %i',size(Hfine)));

% Create U/V point depths
disp('Initializing thin walls on fine grid ... ');
Hx=[];
Hy=[];
%[Hx,Hy]=gen_hxhy(Hfine);

% Recursively coarsen grid
disp('Recursively coarsening from fine grid to target ...');
%[H2,x2,y2]=xyrecur_max(Hfine,xcfine,ycfine,gridscalepow);
%[H2]=xyrecur_ave(Hfine,nfinepow);
clear H2 Hx2 Hy2
for k=1:size(Hfine,3);
 [H2(:,:,k),Hx2(:,:,k),Hy2(:,:,k)]=xyrecur_fv(Hfine(:,:,k),Hx,Hy,nfinepow);
end
clear Hfine Hx Hy 
save H2.mat H2 Hx2 Hy2

% repeat for shelf ice topography Z
disp('Shelfice: Adding halo region ...');
Zfine=exch_llc(Zfine,2^(nfinepow));
disp('Limiting to ocean (h<=0) ...');
Zfine=min(Zfine,0);
disp(sprintf('Fine grid+halo dimensions are %i x %i x %i',size(Zfine)));
% Create U/V point depths
disp('Initializing thin walls on fine grid ... ');
Zx=[];
Zy=[];
% Recursively coarsen grid
disp('Recursively coarsening from fine grid to target ...');
clear ZZ2 Z2 Zx2 Zy2
% hack to be able to use xyrecur_fv for shelfice topography: turn
% everything upside down
ZLARGE = 10000;
ZZ = (-ZLARGE-Zfine).*mask(Zfine);
for k=1:size(Zfine,3);
%  [Z2(:,:,k),Zx2(:,:,k),Zy2(:,:,k)]=xyrecur_fv(Zfine(:,:,k),Zx,Zy,nfinepow);
  [ZZ2(:,:,k),Zx2(:,:,k),Zy2(:,:,k)]=xyrecur_fv(ZZ(:,:,k),Zx,Zy,nfinepow);
end
zmsk = mask(ZZ2);
Z2 = -(ZZ2+ZLARGE).*zmsk;
clear Zfine Zx Zy 
save Z2.mat Z2 Zx2 Zy2

% ~skipfv
else
load H2
load Z2
end

% Drop halo
disp(sprintf('Removing halo: remaining grid is (%i,%i)',nxc,nyc));
H2=H2(2:end-1,2:end-1,:);
Hx2=Hx2(2:end-1,2:end-1,:);
Hy2=Hy2(2:end-1,2:end-1,:);

load VGRID.mat dz
hmax=sum(dz)
H2=max(H2,-hmax);
Hx2=max(Hx2,-hmax);
Hy2=max(Hy2,-hmax);
save HH2.mat H2 Hx2 Hy2

% Fill or empty top N layers
load VGRID.mat dz
N=2;
Dz=[sum(dz(1:N)) dz((N+1):end)]
zround=70; %dz(1);
zround=Dz(1);
[msk,ddz] = hfac(Dz,Hx2,0,0);
if ndims(ddz)==4
 ddz(:,:,:,1)=min(1,round( ddz(:,:,:,1)/zround ))*Dz(1);
else
 ddz(:,:,1)=min(1,round( ddz(:,:,1)/zround ))*Dz(1);
end
Hx2=squeeze(-sum(ddz,ndims(ddz)));
[msk,ddz] = hfac(Dz,Hy2,0,0);
if ndims(ddz)==4
 ddz(:,:,:,1)=min(1,round( ddz(:,:,:,1)/zround ))*Dz(1);
else
 ddz(:,:,1)=min(1,round( ddz(:,:,1)/zround ))*Dz(1);
end
Hy2=squeeze(-sum(ddz,ndims(ddz)));

% Remove inlets
[H2,Hx2,Hy2]=remove_inlets(H2,Hx2,Hy2);

% Fill pot holes
%[H2,Hx2,Hy2]=fill_potholes(H2,Hx2,Hy2);

save h2 H2 Hx2 Hy2

% Re-combine thin wall data (not used if keeping thin walls)
if exist('Hx2')
  hn=drop_thinwalls(H2,Hx2,Hy2,periodic);
else
  hn=H2;
end

% for shelf ice do not do as much
Z2=Z2(2:end-1,2:end-1,:);
Zx2=Zx2(2:end-1,2:end-1,:);
Zy2=Zy2(2:end-1,2:end-1,:);
if exist('Zx2')
  zn=drop_thinwalls(Z2,Zx2,Zy2,periodic);
else
  zn=Z2;
end

% Round to nearest "res" metres
%% res=5;
%% hn=res*round(hn/res);

save HNo.mat hn

% Limit deepest point
load VGRID.mat dz
hn(find(hn<-sum(dz)))=-sum(dz);
% Level land to zero
hn(find(hn>0))=0;

% Fill or empty top N layers
N=2;
Dz=[sum(dz(1:N)) dz((N+1):end)]
zround=25; %dz(1);
zround=Dz(1);
[msk,ddz] = hfac(Dz,hn,0,0);
if ndims(ddz)==4
 ddz(:,:,:,1)=min(1,round( ddz(:,:,:,1)/zround ))*Dz(1);
else
 ddz(:,:,1)=min(1,round( ddz(:,:,1)/zround ))*Dz(1);
end
hn=squeeze(-sum(ddz,ndims(ddz)));

% Remove single point holes (new depth matches nearest depth)
hn=max(hn, min(  min(hn([end 1:end-1],:,:),hn([2:end 1],:,:)) ,...
                 min(hn(:,[1 1:end-1],:),hn(:,[2:end end],:)) ) );

% Automatically discard lakes and unattached oceans
%%clear all
%load HNo
msk=mask(hn);
if size(hn,1) ==360
  nstartpts = 10;
else
  nstartpts = 5;
end
for t=1:size(hn,3);
  msk(:,:,t)=findocean( msk(:,:,t) ,nstartpts, periodic);
end

msk=findocean_fast( msk );
hn=hn.*msk;

% now that we have the basic topography we need to decide which portions
% of the cap-faces are used to generate the final cap-face 
% the idea is to use the triangle that is adjacent to the side face:
[nx,ny]=size(hn);
nf = nx/4;
newface = 0;
clear sh
for k=1:4
  face{k} = hn([1:nf] + (k-1)*nf,end-nf+1:end);
  tmp = face{k}*0;
  imin = 3;
  for j=1:(nf/2+1)
    i = max(1,j-imin):min(nf-j+1+imin,nf);
    tmp(i,j) = face{k}(i,j);
  end
  tmp = rot90(tmp,k-1);
%  tmp= rot90(face{k},k-1);
%  sh(k)=subplot(1,4,k);pcol(sq(tmp)'); caxis([-4250 0]); axis image
  newface = min(newface,tmp);  
end
hn0=hn;
h00=hn;
for k=1:4
  hn0([1:nf] + (k-1)*nf,end-nf+1:end)=rot90(newface,1-k);
  fac = 0; if k==2; fac = 1; end
  h00([1:nf] + (k-1)*nf,end-nf+1:end)=rot90(newface,1-k)*fac;
end
hn=hn0;

% Edit
% Fill oceans to sea-level
%hn(find(hn<0))=0;
% remove artic from 106 to 270 degrees
%hn0=hn;

save HN.mat hn
save ZN.mat zn


% Save bathymetry in bathymetry.bin
load FMT
%delete bathymetry.bin
%wrda('bathymetry.bin',permute(hn,[1 3 2]),1,fmt,Ieee);
%wrda('bathymetryHx.bin',Hx2,1,fmt,Ieee);
%wrda('bathymetryHy.bin',Hy2,1,fmt,Ieee);

gen_masks

% correct shelfice topography
load ZN zn
zn = zn.*msk(:,:,1);
izn=find((zn-hn).*msk(:,:,1) < 0);
zn(izn) = hn(izn);
save ZN zn

% divide:
hnz = hn;
hnz(find(zn==0))=0;
hn(find(zn~=0))=0;
save HN hn
save ZN zn hnz

gen_masks

load HGRID xg yg
xxg=xg(1:end-1,2:end-1);
yyg=yg(1:end-1,2:end-1);
pcol(xxg',yyg',sq(hn)') 
axis([min(xg(:)) max(xg(:)) min(yg(:)) max(yg(:))])

set(gcf,'renderer','painters');
if size(hn,1) == 180
  title('2x2 lat-lon-cap grid')
  print -depsc grid2x2.eps
end
if size(hn,1) == 360
  title('1x1 lat-lon-cap grid')
  print -depsc grid1x1.eps
end