% extract Eta and S/T/U/V/W for Kuroshio Extension
% 25N to 40N, 155E to 175E
% define desired region
region_name='KuroshioExt';
minlat=25;
maxlat=40;
minlon=155;
maxlon=175;
nx=4320;
prec='real*4';
kx=1:40;

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

% Get and save grid information
close all
pin='~dmenemen/llc_4320/grid/';
pout=['~dmenemen/llc_4320/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, no direction
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

% Grid cell center, with direction
fnx='DXF';
fny='DYF';
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);
    end
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest edge points, with 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, with direction
fnx='RAW';
fny='RAS';
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, with direction
fnx='hFacW';
fny='hFacS';
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,kx,ix,jx);
        fldy=read_llc_fkij(finy,nx,fc,kx,ix,jx);
      case {4,5}
        fldx=read_llc_fkij(finy,nx,fc,kx,ix,jx);
        fldy=read_llc_fkij(finx,nx,fc,kx,ix,jx-1); % <<<<<<<<
    end
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest corner (vorticity) points, with direction
fnx='DXV';
fny='DYU';
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-1); % <<<<<<<<
    end
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest edge points, with 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);

% get and save regional Eta
pin='~dmenemen/llc_4320/MITgcm/run';
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
for fnm={'Eta'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=10368:144:1366128, mydisp(ts)
        if ts<485568
            fin=[pin '/' fnm{1} '.' myint2str(ts,10) '.data'];
        else
            fin=[pin '_485568/' fnm{1} '.' myint2str(ts,10) '.data'];
        end
        dy=ts2dte(ts,25,2011,9,10,30);
        fout=[fnm{1} '_' int2str(length(ix)) 'x' int2str(length(jx)) '.' dy];
        fld=read_llc_fkij(fin,nx,fc,1,ix,jx);
        writebin(fout,fld);
    end
end

% get and save regional S/T/W
for fnm={'Salt','Theta','W'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=10368:144:1366128, mydisp(ts)
        if ts<485568
            fin=[pin '/' fnm{1} '.' myint2str(ts,10) '.data'];
        else
            fin=[pin '_485568/' fnm{1} '.' myint2str(ts,10) '.data'];
        end
        dy=ts2dte(ts,25,2011,9,10,30);
        fout=[fnm{1} '_' int2str(length(ix)) 'x' int2str(length(jx)) 'x' int2str(length(kx)) '.' dy];
        for k=1:length(kx); mydisp(k)
            fld=read_llc_fkij(fin,nx,fc,kx(k),ix,jx);
            writebin(fout,fld,1,'real*4',k-1);
        end
    end
end

% get and save regional fields of U and V
% note that zonal velocity is U in faces 1/2 and V in faces 4/5
% and meridional velocity is V in faces 1/2 and -U in faces 4/5
eval(['mkdir ' pout 'U'])
eval(['mkdir ' pout 'V'])
eval(['cd ' pout])
for ts=10368:144:1366128, mydisp(ts)
    if ts<485568
        finu=[pin '/U.' myint2str(ts,10) '.data'];
        finv=[pin '/V.' myint2str(ts,10) '.data'];
    else
        finu=[pin '_485568/U.' myint2str(ts,10) '.data'];
        finv=[pin '_485568/V.' myint2str(ts,10) '.data'];
    end
    dy=ts2dte(ts,25,2011,9,10,30);
    foutu=['U/U_' int2str(length(ix)) 'x' int2str(length(jx)) 'x' int2str(length(kx)) '.' dy];
    foutv=['V/V_' int2str(length(ix)) 'x' int2str(length(jx)) 'x' int2str(length(kx)) '.' dy];
    for k=1:length(kx); mydisp(k)
        fld=read_llc_fkij(finv,nx,fc,kx(k),ix,jx);
        writebin(foutu,fld,1,'real*4',k-1);
        fld=-read_llc_fkij(finu,nx,fc,kx(k),ix,jx-1);
        writebin(foutv,fld,1,'real*4',k-1);
    end
end