% define desired region
nx=2160;
region_name='Antarctic';
minlat=-80;
maxlat=-50;

% extract indices for desired region
prec='real*4';
gdir='~/llc_2160/grid/';
fnam=[gdir 'Depth.data'];
[tmp fc ix jx]=quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat);

% make sure jx is same for all faces
jx=min([jx{[1:5]}]):max([jx{[1:5]}]);
n=length(jx);

% map field on array and plot
m(1)=0;
for f=1:length(fc)
    m(f+1)=length(ix{fc(f)});
end
fld=zeros(sum(m),n);
for f=1:length(fc)
    fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx);
end
quikpcolor(fld')

% Get and save grid information
pout=['~/llc_2160/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
for fnm={'AngleCS','AngleSN','DXC','DXG','DYC','DYG','Depth', ...
         'RAC','RAS','RAW','RAZ','U2zonDir','V2zonDir', ...
         'XC','XG','YC','YG','hFacC','hFacS','hFacW'}
    fin=[gdir fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    switch fnm{1}
      case{'hFacC','hFacS','hFacW'}
        fld=zeros(sum(m),n,length(kx));
        for f=1:length(fc)
            fld((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
                read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx);
        end
        fout=[fout 'x' int2str(length(kx))];
      otherwise
        fld=zeros(sum(m),n);
        for f=1:length(fc)
            fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx);
        end
    end
    writebin(fout,fld);
end