% NWAustralia (Northwest Australian Shelf) SWOT Crossover Region
% oceFWflx, oceQnet, oceQsw, oceSflux, oceTAUX, oceTAUY
% Eta, KPPhbl, PhiBot, Salt, Theta, U, V, W
% extract 13-Sep-2011 to 15-Nov-2012
% lats -15 to -11N, lons 120.5 to 124.5E
% (example extraction on face 1 or 2)

% {{{ define desired region
region_name='NWAustralia';
minlat=-15;
maxlat=-11;
minlon=120.5;
maxlon=124.5;
mints=dte2ts('13-Sep-2011',25,2011,9,10);
maxts=dte2ts('15-Nov-2012',25,2011,9,10);
nx=4320;
prec='real*4';
% }}}

% {{{ 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');
load ~dmenemen/llc_4320/grid/thk90.mat
bot90=dpt90+thk90/2; bot90(end)=bot90(end)+1;
kx=1:min(find(bot90>mmax(fld)));
suf1=['_' int2str(length(ix)) 'x' int2str(length(jx))];
suf2=[suf1 'x' int2str(length(kx))];
close all
% }}}

% {{{ get and save grid information
pin='~dmenemen/llc_4320/grid/';
pout=['~dmenemen/llc_4320/regions/Crossover/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
% {{{ All grid except angles
for fnm={'Depth','RAC','XC','YC','hFacC','XG','YG', ...
         'RAZ','DXC','DYC','DXV','DYU','DXG','DYG'}
    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
% }}}
% {{{ AngleCS and AngleSN at grid cell centers
fnx='AngleCS';
fny='AngleSN';
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);
% }}}
% }}}

% {{{ create commands for extracting model output fields
pout=['~dmenemen/llc_4320/regions/Crossover/' region_name '/'];
extract='/home4/bcnelson/MITgcm/extract/latest/extract4320 -g ';
timesteps=[int2str(mints) '-' int2str(maxts) ' '];
startPoint=[int2str((fc-1)*nx+min(ix)) ',' int2str(min(jx)) ',1 '];

% {{{ get and save regional 2D fields
extent=[int2str(length(ix)) ',' int2str(length(jx)) ',1'];
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw', ...
         'oceSflux','oceTAUX','oceTAUY'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pout fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end
% }}}

% {{{ get and save regional 3D fields
extent=[int2str(length(ix)) ',' int2str(length(jx)) ',' int2str(kx(end))];
for fnm={'Salt','Theta','U','V','W'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pout fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end
% }}}
% }}}