% extract Eta and S/T/U/V/W for SPURS1
% 23.5N to 25.5N, 39.5W to 36.5W
% define desired region
region_name='SPURS1';
minlat=23.5;
maxlat=25.5;
minlon=-39.5;
maxlon=-36.5;
mints=dte2ts('13-Sep-2011',25,2011,9,10);
maxts=dte2ts('15-Nov-2012',25,2011,9,10);
nx=4320;
prec='real*4';
kx=1:89;

% extract indices for desired region
gdir='~dmenemen/llc_4320/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')
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])

% Grid cell center
for fnm={'Depth','RAC','XC','YC','hFacC'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    switch fnm{1}
      case{'hFacC'}
        fld=zeros(sum(m),n,length(kx));
        fld(1:m(2),:,:)         =read_llc_fkij(fin,nx,fc(2),kx,ix{fc(2)},jx{fc(2)});
        fld((m(2)+1):sum(m),:,:)=read_llc_fkij(fin,nx,fc(1),kx,ix{fc(1)},jx{fc(1)});
        fout=[fout 'x' int2str(length(kx))];
      otherwise
        fld=zeros(sum(m),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)});
    end
    writebin(fout,fld);
end

% Southwest corner (vorticity) points, no direction
fld=zeros(sum(m),n);
for fnm={'XG','YG','RAZ'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    fld((m(2)+1):sum(m),:)=read_llc_fkij(fin,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
    fld(1:m(2),:)         =read_llc_fkij(fin,nx,fc(2),1,ix{fc(2)},jx{fc(2)}-1); % <<<<<<<<
    writebin(fout,fld);
end

% West edge points, no direction
fnx='DXC';
fny='DYC';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
fldx=zeros(sum(m),n);
fldy=zeros(sum(m),n);
fldx((m(2)+1):sum(m),:)=read_llc_fkij(finx,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
fldx(1:m(2),:)         =read_llc_fkij(finy,nx,fc(2),1,ix{fc(2)},jx{fc(2)});
fldy((m(2)+1):sum(m),:)=read_llc_fkij(finy,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
fldy(1:m(2),:)         =read_llc_fkij(finx,nx,fc(2),1,ix{fc(2)},jx{fc(2)}-1); % <<<<<<<<
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest corner (vorticity) points, no direction
fnx='DXV';
fny='DYU';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
fldx((m(2)+1):sum(m),:)=read_llc_fkij(finx,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
fldx(1:m(2),:)         =read_llc_fkij(finy,nx,fc(2),1,ix{fc(2)},jx{fc(2)});
fldy((m(2)+1):sum(m),:)=read_llc_fkij(finy,nx,fc(1),1,ix{fc(1)},jx{fc(1)}-1); % <<<<<<<<
fldy(1:m(2),:)         =read_llc_fkij(finx,nx,fc(2),1,ix{fc(2)},jx{fc(2)}-1); % <<<<<<<<
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 '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
fldx((m(2)+1):sum(m),:)=read_llc_fkij(finx,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
fldx(1:m(2),:)         =read_llc_fkij(finy,nx,fc(2),1,ix{fc(2)},jx{fc(2)});
fldy((m(2)+1):sum(m),:)=read_llc_fkij(finy,nx,fc(1),1,ix{fc(1)},jx{fc(1)}-1); % <<<<<<<<
fldy(1:m(2),:)         =read_llc_fkij(finx,nx,fc(2),1,ix{fc(2)},jx{fc(2)});
writebin(foutx,fldx);
writebin(fouty,fldy);

% extract model output fields
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
extract='/home4/bcnelson/MITgcm/extract/v1.08/extract4320 ';
timesteps=[int2str(mints) '-' int2str(maxts) ' '];

% get and save regional Eta and PhiBot
startPoint=[int2str(3*nx+min(ix{5})) ',' int2str(min(jx{5})) ',1 '];
extent=[int2str(sum(m)) ',' int2str(n) ',1'];
for fnm={'Eta','PhiBot'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    fieldNames=[fnm{1} ' '];
    [s w]=system([extract timesteps  fieldNames  startPoint  extent]);
end

% get and save regional S/T/W
startPoint=[int2str(3*nx+min(ix{5})) ',' int2str(min(jx{5})) ',1 '];
extent=[int2str(sum(m)) ',' int2str(n) ',' int2str(kx(end))];
for fnm={'Salt','Theta','W'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    fieldNames=[fnm{1} ' '];
    [s w]=system([extract timesteps  fieldNames  startPoint  extent]);
end

% example parallel extraction in an interactive queue
% ~bcnelson/MITgcm/extract/v1.08/extract4320 -g 10368-1492992 Salt 17209,9179,1 144,115,89 > joblist
% parallel --slf $PBS_NODEFILE -j2 -a joblist

~bcnelson/MITgcm/extract/v1.08/extract4320 -g 908640-1492992 W 17209,9179,1 144,115,89 > joblist

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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])
fldu=zeros(sum(m),n);
fldv=zeros(sum(m),n);
for ts=10368:144:485568, disp(ts)
    finu=[pin 'U.' myint2str(ts,10) '.data'];
    finv=[pin 'V.' myint2str(ts,10) '.data'];
    dy=ts2dte(ts,25,2011,9,10,30);
    foutu=['U/U_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    foutv=['V/V_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    for k=1:length(kx); mydisp(k)
        fldu(1:m(2),:)         =read_llc_fkij(finv,nx,fc(2),1,ix{fc(2)},jx{fc(2)});
        fldu((m(2)+1):sum(m),:)=read_llc_fkij(finu,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
        fldv(1:m(2),:)        =-read_llc_fkij(finu,nx,fc(2),1,ix{fc(2)},jx{fc(2)}-1);
        fldv((m(2)+1):sum(m),:)=read_llc_fkij(finv,nx,fc(1),1,ix{fc(1)},jx{fc(1)});
        writebin(foutu,fldu,1,'real*4',k-1);
        writebin(foutv,fldv,1,'real*4',k-1);
    end
end