% initialize
pin='~dmenemen/llc_2160/regions/latlon/'; % model output files location
eval(['cd ' pin 'drift'])
nx=8640; ny=4320;                         % model grid dimensions
suf=['_' int2str(nx) 'x' int2str(ny)];
for fld={'XC','YC','hFacC'}               % load model grid info arrays
  fnm=[pin 'grid/' fld{1} suf];
  eval([fld{1} '=readbin(fnm,[nx ny]);'])
end
XC(find(XC<XC(1)))=XC(find(XC<XC(1)))+360;% unwrap longitude array
[J I]=meshgrid(1:ny,1:nx);

% create matrix of barrier zones as per Duke et al. 2002, Fig. 2
Regions=int8(ones(nx,ny));
Regions(find(hFacC==0))=-1;
% zone 2
Regions(7150:end,2980:end)=2;
Regions(7370:end,2910:2980)=2;
Regions(7530:end,2780:2910)=2;
Regions(7555:end,2760:2780)=2;
Regions(7575:7626,2749:2760)=2;
Regions(7680:7860,2725:2760)=2;
Regions(7880:end,1:2760)=2;
Regions(1:400,1:2500)=2;
% zone 3
ix=find(I>1&I<1400&(Regions==1|Regions==-1));
Regions(ix)=3;
% zone 4
ix=find(I>1400&I<2380&(Regions==1|Regions==-1));
Regions(ix)=4;
% zone 5
Regions(2380:3100,1:2000)=5;
Regions(2380:3300,2001:2060)=5;
Regions(2380:3750,2061:2150)=5;
Regions(2380:3800,2151:2190)=5;
Regions(2380:3970,2191:2350)=5;
Regions(2380:3910,2351:2450)=5;
Regions(2380:4010,2451:2929)=5;
Regions(2380:4260,2930:3100)=5;
Regions(2380:5100,3101:end)=5;
% zone 6
ix=find(I<5250&J>2500&(Regions==1|Regions==-1));
Regions(ix)=6;
ix=find(I<6400&J<2501&(Regions==1|Regions==-1));
Regions(ix)=6;
% zone 1
ix=find((Regions<2));
Regions(ix)=1;

clf, mypcolor(Regions'), thincolorbar
save global/DukeRegions2002 Regions