% Gulf Stream region for Peter Cornillon
% extract SST
% lats 25N to 50N, lons 278E to 320E
% (example extraction on face 5, i.e., rotated UV fields)

% define desired region
nx=1080;
prec='real*4';
region_name='GulfStream2';
minlat=25;
maxlat=50;
minlon=278;
maxlon=320;

% extract indices for desired region
gdir='~dmenemen/llc_1080/grid/';
fnam=[gdir 'Depth.data'];
[fld fc ix jx] = ...
    quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);
quikpcolor(fld'); caxis([0 1])
kx=1;

% Get and save grid information
close all
pin='~dmenemen/llc_1080/grid/';
pout=['~dmenemen/llc_1080/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
suf1=['_' int2str(length(ix)) 'x' int2str(length(jx))];
suf2=[suf1 'x' int2str(length(kx))];

% Grid cell center
for fnm={'Depth','RAC','XC','YC','hFacC'}
    fin=[pin fnm{1} '.data'];
    switch fnm{1}
      case{'hFacC'}
        fld=read_llc_fkij(fin,nx,fc,kx,ix,jx);
        fout=[fnm{1} suf2];
      otherwise
        fld=read_llc_fkij(fin,nx,fc,1,ix,jx);
        fout=[fnm{1} suf1];
    end
    writebin(fout,fld);
end

% Southwest corner (vorticity) points, no direction
for fnm={'XG','YG','RAZ'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} suf1];
    for f=1:length(fc)
        switch fc(f)
          case {1,2}
            fld=read_llc_fkij(fin,nx,fc,1,ix,jx);
          case {4,5}
            fld=read_llc_fkij(fin,nx,fc,1,ix,jx-1); % <<<<<<<<
        end
    end
    writebin(fout,fld);
end

% West edge points, no direction
fnx='DXC';
fny='DYC';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx=read_llc_fkij(finx,nx,fc,1,ix,jx);
        fldy=read_llc_fkij(finy,nx,fc,1,ix,jx);
      case {4,5}
        fldx=read_llc_fkij(finy,nx,fc,1,ix,jx);
        fldy=read_llc_fkij(finx,nx,fc,1,ix,jx-1); % <<<<<<<<
    end
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest edge points, no direction
fnx='DXG';
fny='DYG';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx=read_llc_fkij(finx,nx,fc,1,ix,jx);
        fldy=read_llc_fkij(finy,nx,fc,1,ix,jx);
      case {4,5}
        fldx=read_llc_fkij(finy,nx,fc,1,ix,jx-1); % <<<<<<<<
        fldy=read_llc_fkij(finx,nx,fc,1,ix,jx);
    end
end
writebin(foutx,fldx);
writebin(fouty,fldy);