mlosch/interp_llc/topo_s2004.m

% 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
  
1	mlosch	1.1	% This is a script that executes a series of commands for
2			% extracting GEBCO data, interpolating, editing and then
3			% generating a mask-file.
4			%
5			% Uses data files: HGRID.mat VGRID.mat
6			% Creates data files: bathymetry.bin
7			% Creates arrays: n/a
8			% Uses m-scripts: setgrid grid_sphere finegrid gen_hxhy xyrecur_* edtopo
9			%
10			% Created 08/15/99 by adcroft@mit.edu
11			% Modified 11/09/99 by adcroft@mit.edu
12			% Maintained by adcroft@mit.edu, abiastoch@ucsd.edu
13			% Modifield to read etopo2 01/10/04 by mlosch@awi-bremerhaven.de
14
15			%clear all
16			format compact
17			more off
18
19			skipfv=0
20
21			load HGRID.mat nxc nyc lon_lo lon_hi lat_lo lat_hi periodic
22
23			if ~skipfv
24
25			load HFGRID xcf
26			% Get the GEBCO data on a conveniently sized grid
27			disp('Extracting S2004 data ... ');
28			[He2,xe2,ye2]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi);
29			i0=find(xe2>=360+s_lon);
30			xe2(i0) = xe2(i0)-360;
31			if xe2(end)==xe2(1); xe2(end)=xe2(1)+360; end
32			% make the domain bigger so that it covers everything
33			[nxe,nye]=size(He2);
34			imax = min(find(xe2+360>max(xcf(:))));
35			imin = max(find(xe2-360<min(xcf(:))));
36			He2 = [He2(imin+1:end-1,:); He2; He2(2:imax,:)];
37			xe2 = [xe2(imin+1:end-1)-360 xe2 xe2(2:imax)+360];
38
39			disp(sprintf('S2004 dimensions are %i x %i',size(He2)));
40			save S2004.mat He2 xe2 ye2
41
42			disp('Extracting BRIOS data ... ');
43			[Hb2,Zb2,xb2,yb2]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi);
44			i0 = find(xb2>=360+s_lon);
45			xb2 = [xb2(i0)-360; xb2];
46			Hb2 = [Hb2(i0,:); Hb2];
47			Zb2 = [Zb2(i0,:); Zb2];
48
49			disp(sprintf('BRIOS dimensions are %i x %i',size(Hb2)));
50			save BRIOS.mat Hb2 Zb2 xb2 yb2
51
52			% Interpolate to fine grid
53			load HFGRID xcf ycf nfinepow
54			disp('Interpolating to fine grid ... ');
55			Hfine = 0*xcf;
56			Bfine = 0*xcf;
57			Zfine = 0*xcf;
58			% this may be replace with something more sophisticated (and expensive)
59			imeth = 'nearest';
60			imeth = 'linear';
61			disp(sprintf('interpolation method: %s',imeth))
62			largeproblem = 1;
63			if largeproblem
64			% divide and conquer
65			[nxt nyt]=size(xcf);
66			nt = 2.^max(nfinepow-2,0);
67			nxt = nxt./nt;
68			nyt = nyt./nt;
69			for j=1:nt
70			jy = nyt*(j-1)+[1:nyt];
71			for i=1:nt
72			ix = nxt*(i-1)+[1:nxt];
73			disp(sprintf('tile number = (%u,%u)',i,j));
74			Hfine(ix,jy)=interp2(xe2(:),ye2(:)',He2',xcf(ix,jy),ycf(ix,jy),imeth);
75			Bfine(ix,jy)=interp2(xb2(:),[yb2;-yb2(end:-1:1)]', ...
76			[Hb2 zeros(size(Hb2))]',xcf(ix,jy),ycf(ix,jy),imeth);
77			Zfine(ix,jy)=interp2(xb2(:),[yb2;-yb2(end:-1:1)]', ...
78			[Zb2 zeros(size(Zb2))]', xcf(ix,jy),ycf(ix,jy),imeth);
79			end
80			end
81			else
82			% do everything in one go
83			Hfine=interp2(xe2(:),ye2(:)',He2',xcf,ycf,imeth);
84			Bfine=interp2(xb2(:),[yb2;-yb2(end:-1:1)]',[Hb2 zeros(size(Hb2))]', ...
85			xcf,ycf,imeth);
86			Zfine=interp2(xb2(:),[yb2;-yb2(end:-1:1)]',[Zb2 zeros(size(Zb2))]', ...
87			xcf,ycf,'nearest');
88			end
89
90			Bfine(isnan(Bfine)) = 0;
91			Zfine(isnan(Zfine)) = 0;
92			% add bathymetry under shelfice
93			msk0 = Hfine;
94			msk0(find(msk0>0))=1;
95			sizmsk = Zfine;
96			sihmsk = Bfine;
97			sihmsk(find(sihmsk~=0))=1;
98			sizmsk(find(sizmsk~=0))=1;
99			%sizmsk = sizmsk.*sihmsk;
100			% introduce wall at the south (affects only Ross Sea)
101			Zfine = Zfine.*sizmsk;
102			% under ice bottom topography
103			i1 = find(sizmsk == 1);
104			Hfine(i1) = Bfine(i1);
105			%% Handarbeit:
106			%x0=
107			%y0=
108			ierr = find(Zfine(i1)<Hfine(i1));
109
110			clear Bfine Hb2 xe2 ye2
111
112			%load HFINE
113			disp(sprintf('Fine grid dimensions are %i x %i x %i',size(Hfine)));
114			save HFINE.mat Hfine xcf ycf
115			save ZFINE.mat Zfine xcf ycf
116			clear He2 xe2 ye2 xcf ycf
117
118			disp('Adding halo region ...');
119			load HFGRID nfinepow
120			Hfine=exch_llc(Hfine,2^(nfinepow));
121			disp('Limiting to ocean (h<=0) ...');
122			Hfine=min(Hfine,0);
123			disp(sprintf('Fine grid+halo dimensions are %i x %i x %i',size(Hfine)));
124
125			% Create U/V point depths
126			disp('Initializing thin walls on fine grid ... ');
127			Hx=[];
128			Hy=[];
129			%[Hx,Hy]=gen_hxhy(Hfine);
130
131			% Recursively coarsen grid
132			disp('Recursively coarsening from fine grid to target ...');
133			%[H2,x2,y2]=xyrecur_max(Hfine,xcfine,ycfine,gridscalepow);
134			%[H2]=xyrecur_ave(Hfine,nfinepow);
135			clear H2 Hx2 Hy2
136			for k=1:size(Hfine,3);
137			[H2(:,:,k),Hx2(:,:,k),Hy2(:,:,k)]=xyrecur_fv(Hfine(:,:,k),Hx,Hy,nfinepow);
138			end
139			clear Hfine Hx Hy
140			save H2.mat H2 Hx2 Hy2
141
142			% repeat for shelf ice topography Z
143			disp('Shelfice: Adding halo region ...');
144			Zfine=exch_llc(Zfine,2^(nfinepow));
145			disp('Limiting to ocean (h<=0) ...');
146			Zfine=min(Zfine,0);
147			disp(sprintf('Fine grid+halo dimensions are %i x %i x %i',size(Zfine)));
148			% Create U/V point depths
149			disp('Initializing thin walls on fine grid ... ');
150			Zx=[];
151			Zy=[];
152			% Recursively coarsen grid
153			disp('Recursively coarsening from fine grid to target ...');
154			clear ZZ2 Z2 Zx2 Zy2
155			% hack to be able to use xyrecur_fv for shelfice topography: turn
156			% everything upside down
157			ZLARGE = 10000;
158			ZZ = (-ZLARGE-Zfine).*mask(Zfine);
159			for k=1:size(Zfine,3);
160			% [Z2(:,:,k),Zx2(:,:,k),Zy2(:,:,k)]=xyrecur_fv(Zfine(:,:,k),Zx,Zy,nfinepow);
161			[ZZ2(:,:,k),Zx2(:,:,k),Zy2(:,:,k)]=xyrecur_fv(ZZ(:,:,k),Zx,Zy,nfinepow);
162			end
163			zmsk = mask(ZZ2);
164			Z2 = -(ZZ2+ZLARGE).*zmsk;
165			clear Zfine Zx Zy
166			save Z2.mat Z2 Zx2 Zy2
167
168			% ~skipfv
169			else
170			load H2
171			load Z2
172			end
173
174			% Drop halo
175			disp(sprintf('Removing halo: remaining grid is (%i,%i)',nxc,nyc));
176			H2=H2(2:end-1,2:end-1,:);
177			Hx2=Hx2(2:end-1,2:end-1,:);
178			Hy2=Hy2(2:end-1,2:end-1,:);
179
180			load VGRID.mat dz
181			hmax=sum(dz)
182			H2=max(H2,-hmax);
183			Hx2=max(Hx2,-hmax);
184			Hy2=max(Hy2,-hmax);
185			save HH2.mat H2 Hx2 Hy2
186
187			% Fill or empty top N layers
188			load VGRID.mat dz
189			N=2;
190			Dz=[sum(dz(1:N)) dz((N+1):end)]
191			zround=70; %dz(1);
192			zround=Dz(1);
193			[msk,ddz] = hfac(Dz,Hx2,0,0);
194			if ndims(ddz)==4
195			ddz(:,:,:,1)=min(1,round( ddz(:,:,:,1)/zround ))*Dz(1);
196			else
197			ddz(:,:,1)=min(1,round( ddz(:,:,1)/zround ))*Dz(1);
198			end
199			Hx2=squeeze(-sum(ddz,ndims(ddz)));
200			[msk,ddz] = hfac(Dz,Hy2,0,0);
201			if ndims(ddz)==4
202			ddz(:,:,:,1)=min(1,round( ddz(:,:,:,1)/zround ))*Dz(1);
203			else
204			ddz(:,:,1)=min(1,round( ddz(:,:,1)/zround ))*Dz(1);
205			end
206			Hy2=squeeze(-sum(ddz,ndims(ddz)));
207
208			% Remove inlets
209			[H2,Hx2,Hy2]=remove_inlets(H2,Hx2,Hy2);
210
211			% Fill pot holes
212			%[H2,Hx2,Hy2]=fill_potholes(H2,Hx2,Hy2);
213
214			save h2 H2 Hx2 Hy2
215
216			% Re-combine thin wall data (not used if keeping thin walls)
217			if exist('Hx2')
218			hn=drop_thinwalls(H2,Hx2,Hy2,periodic);
219			else
220			hn=H2;
221			end
222
223			% for shelf ice do not do as much
224			Z2=Z2(2:end-1,2:end-1,:);
225			Zx2=Zx2(2:end-1,2:end-1,:);
226			Zy2=Zy2(2:end-1,2:end-1,:);
227			if exist('Zx2')
228			zn=drop_thinwalls(Z2,Zx2,Zy2,periodic);
229			else
230			zn=Z2;
231			end
232
233			% Round to nearest "res" metres
234			%% res=5;
235			%% hn=res*round(hn/res);
236
237			save HNo.mat hn
238
239			% Limit deepest point
240			load VGRID.mat dz
241			hn(find(hn<-sum(dz)))=-sum(dz);
242			% Level land to zero
243			hn(find(hn>0))=0;
244
245			% Fill or empty top N layers
246			N=2;
247			Dz=[sum(dz(1:N)) dz((N+1):end)]
248			zround=25; %dz(1);
249			zround=Dz(1);
250			[msk,ddz] = hfac(Dz,hn,0,0);
251			if ndims(ddz)==4
252			ddz(:,:,:,1)=min(1,round( ddz(:,:,:,1)/zround ))*Dz(1);
253			else
254			ddz(:,:,1)=min(1,round( ddz(:,:,1)/zround ))*Dz(1);
255			end
256			hn=squeeze(-sum(ddz,ndims(ddz)));
257
258			% Remove single point holes (new depth matches nearest depth)
259			hn=max(hn, min( min(hn([end 1:end-1],:,:),hn([2:end 1],:,:)) ,...
260			min(hn(:,[1 1:end-1],:),hn(:,[2:end end],:)) ) );
261
262			% Automatically discard lakes and unattached oceans
263			%%clear all
264			%load HNo
265			msk=mask(hn);
266			if size(hn,1) ==360
267			nstartpts = 10;
268			else
269			nstartpts = 5;
270			end
271			for t=1:size(hn,3);
272			msk(:,:,t)=findocean( msk(:,:,t) ,nstartpts, periodic);
273			end
274
275			msk=findocean_fast( msk );
276			hn=hn.*msk;
277
278			% now that we have the basic topography we need to decide which portions
279			% of the cap-faces are used to generate the final cap-face
280			% the idea is to use the triangle that is adjacent to the side face:
281			[nx,ny]=size(hn);
282			nf = nx/4;
283			newface = 0;
284			clear sh
285			for k=1:4
286			face{k} = hn([1:nf] + (k-1)*nf,end-nf+1:end);
287			tmp = face{k}*0;
288			imin = 3;
289			for j=1:(nf/2+1)
290			i = max(1,j-imin):min(nf-j+1+imin,nf);
291			tmp(i,j) = face{k}(i,j);
292			end
293			tmp = rot90(tmp,k-1);
294			% tmp= rot90(face{k},k-1);
295			% sh(k)=subplot(1,4,k);pcol(sq(tmp)'); caxis([-4250 0]); axis image
296			newface = min(newface,tmp);
297			end
298			hn0=hn;
299			h00=hn;
300			for k=1:4
301			hn0([1:nf] + (k-1)*nf,end-nf+1:end)=rot90(newface,1-k);
302			fac = 0; if k==2; fac = 1; end
303			h00([1:nf] + (k-1)nf,end-nf+1:end)=rot90(newface,1-k)fac;
304			end
305			hn=hn0;
306
307			% Edit
308			% Fill oceans to sea-level
309			%hn(find(hn<0))=0;
310			% remove artic from 106 to 270 degrees
311			%hn0=hn;
312
313			save HN.mat hn
314			save ZN.mat zn
315
316
317			% Save bathymetry in bathymetry.bin
318			load FMT
319			%delete bathymetry.bin
320			%wrda('bathymetry.bin',permute(hn,[1 3 2]),1,fmt,Ieee);
321			%wrda('bathymetryHx.bin',Hx2,1,fmt,Ieee);
322			%wrda('bathymetryHy.bin',Hy2,1,fmt,Ieee);
323
324			gen_masks
325
326			% correct shelfice topography
327			load ZN zn
328			zn = zn.*msk(:,:,1);
329			izn=find((zn-hn).*msk(:,:,1) < 0);
330			zn(izn) = hn(izn);
331			save ZN zn
332
333			% divide:
334			hnz = hn;
335			hnz(find(zn==0))=0;
336			hn(find(zn~=0))=0;
337			save HN hn
338			save ZN zn hnz
339
340			gen_masks
341
342			load HGRID xg yg
343			xxg=xg(1:end-1,2:end-1);
344			yyg=yg(1:end-1,2:end-1);
345			pcol(xxg',yyg',sq(hn)')
346			axis([min(xg(:)) max(xg(:)) min(yg(:)) max(yg(:))])
347
348			set(gcf,'renderer','painters');
349			if size(hn,1) == 180
350			title('2x2 lat-lon-cap grid')
351			print -depsc grid2x2.eps
352			end
353			if size(hn,1) == 360
354			title('1x1 lat-lon-cap grid')
355			print -depsc grid1x1.eps
356			end
357