| 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 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph