% extract SST in northern extratropical Atlantic
% 288-344E, 9-52N

% define desired region
nx=2160;
prec='real*4';
region_name='NAtlan';
minlat=9;
maxlat=52;
minlon=288;
maxlon=344;

% extract indices for desired region
gdir='~dmenemen/llc_2160/grid/';
fnam=[gdir 'Depth.data'];
[tmp fc ix jx] = ...
    quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);
m=[0 length(ix{fc(2)}) length(ix{fc(1)})];
n=length(jx{fc(1)});
fld=zeros(sum(m),n);
fld(1:m(2),:)=tmp{fc(2)};
fld((m(2)+1):sum(m),:)=tmp{fc(1)};
quikpcolor(fld')
thincolorbar

% get and save grid information
close all
pin='~dmenemen/llc_2160/grid/';
pout=['~dmenemen/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=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    fld(1:m(2),:)         =read_llc_fkij(fin,nx,fc(2),1,ix{fc(2)},jx{fc(2)});
    fld((m(2)+1):sum(m),:)=read_llc_fkij(fin,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
    writebin(fout,fld);
end

% get and save regional 2D fields
pn='~dmenemen/llc_2160/MITgcm/run';
pout=['~dmenemen/llc_2160/regions/' region_name '/'];
fld=zeros(sum(m),n);
for fnm={'Theta'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=92160:80:1586400; disp(ts)
        if (ts<1198080)
            pin=[pn '_day49_624/'];
        else
            pin=[pn '/'];
        end
        fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
        dy=ts2dte(ts,45,2011,1,17,30);
        fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) '.' dy];
        fld(1:m(2),:)         =read_llc_fkij(fin,nx,fc(2),1,ix{fc(2)},jx{fc(2)});
        fld((m(2)+1):sum(m),:)=read_llc_fkij(fin,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
        writebin(fout,fld);
    end
end