% Subtropical Countercurrent (STCC) region for Zhiyu Liu
% oceFWflx, oceQnet, oceQsw, oceTAUX, oceTAUY
% Eta, KPPhbl, PhiBot, Salt, Theta, U, V, W
% extract 13-Sep-2011 to 15-Nov-2012
% lats 30N to 46N, lons 120E to 138E
% example extraction on face 2+4
% extract face 2 and face 4 separately then use cat to concatenate

% {{{ define desired region
region_name='STCC';
minlat=18.32;
maxlat=22.32;
minlon=140.73;
maxlon=144.73;
mints=dte2ts('13-Sep-2011',25,2011,9,10);
maxts=dte2ts('15-Nov-2012',25,2011,9,10);
timesteps=[int2str(mints) '-' int2str(maxts) ' '];
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
pin1=[pout 'fc2/'];
pin2=[pout 'fc4/'];
nx=4320;
prec='real*4';
% }}}

% {{{ 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(1)=0;
for f=1:length(FC)
    m(f+1)=length(IX{FC(f)});
end
n=length(JX{FC(1)});
fld=zeros(sum(m),n);
for f=1:length(FC)
    fld((sum(m(1:f))+1):sum(m(1:(f+1))),:)=tmp{FC(f)};
end
quikpcolor(fld')
fc1=FC(1);
fc2=FC(2);
ix1=IX{fc1};
ix2=IX{fc2};
jx=JX{fc1};
m1=length(ix1);
m2=length(ix2);
load ~dmenemen/llc_4320/grid/thk90.mat
bot90=dpt90+thk90/2; bot90(end)=bot90(end)+1;
kx=1:min(find(bot90>mmax(fld)));
z=length(kx);
clear FC IX JX b* d* f fld fn* m maxl* minl* prec r* tm*
close all
% }}}

% {{{ extract face 2
% {{{ get and save grid information
eval(['mkdir ' pin1 'grid/'])
eval(['cd ' pin1 'grid/'])
suf1=['_' int2str(m1) 'x' int2str(n)];
suf2=[suf1 'x' int2str(z)];
for fnm={'Depth','RAC','XC','YC','hFacC','XG','YG', ...
         'RAZ','DXC','DYC','DXV','DYU','DXG','DYG'}
    fin=[gdir fnm{1} '.data'];
    switch fnm{1}
      case{'hFacC'}
        fld=read_llc_fkij(fin,nx,fc1,kx,ix1,jx);
        fout=[fnm{1} suf2];
      otherwise
        fld=read_llc_fkij(fin,nx,fc1,1,ix1,jx);
        fout=[fnm{1} suf1];
    end
    writebin(fout,fld);
end
% }}}
% {{{ create commands for extracting model output fields
extract='/home4/bcnelson/MITgcm/extract/latest/extract4320 -g ';
startPoint=[int2str((fc1-1)*nx+min(ix1)) ',' int2str(min(jx)) ',1 '];
% }}}
% {{{ get and save regional 2D fields
extent=[int2str(m1) ',' int2str(n) ',1'];
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw','oceTAUX','oceTAUY'}
    eval(['mkdir ' pin1 fnm{1}])
    eval(['cd ' pin1 fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pin1 fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -a joblist')
end
% }}}
% {{{ Get and save regional 3D fields
extent=[int2str(m1) ',' int2str(n) ',' int2str(z)];
for fnm={'Salt','Theta','U','V','W'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pin1 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
% }}}
% }}}

% {{{ extract face 4
% {{{ Get and save grid information
eval(['mkdir ' pin2 'grid/'])
eval(['cd ' pin2 'grid/'])
suf1=['_' int2str(m2) 'x' int2str(n)];
suf2=[suf1 'x' int2str(z)];
% {{{ grid cell center
for fnm={'Depth','RAC','XC','YC','hFacC'}
    fin=[gdir fnm{1} '.data'];
    switch fnm{1}
      case{'hFacC'}
        fld=read_llc_fkij(fin,nx,fc2,kx,ix2,jx);
        fout=[fnm{1} suf2];
      otherwise
        fld=read_llc_fkij(fin,nx,fc2,1,ix2,jx);
        fout=[fnm{1} suf1];
    end
    writebin(fout,fld);
end
% }}}
% {{{ southwest corner (vorticity) points, no direction
for fnm={'XG','YG','RAZ'}
    fin=[gdir fnm{1} '.data'];
    fout=[fnm{1} suf1];
    fld=read_llc_fkij(fin,nx,fc2,1,ix2,jx-1); % <<<<<<<<
    writebin(fout,fld);
end
% }}}
% {{{ west edge points, no direction
fnx='DXC';
fny='DYC';
finx=[gdir fnx '.data'];
finy=[gdir fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
fldx=read_llc_fkij(finy,nx,fc2,1,ix2,jx);
fldy=read_llc_fkij(finx,nx,fc2,1,ix2,jx-1); % <<<<<<<<
writebin(foutx,fldx);
writebin(fouty,fldy);
% }}}
% {{{ southwest corner (vorticity) points, no direction
fnx='DXV';
fny='DYU';
finx=[gdir fnx '.data'];
finy=[gdir fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
fldx=read_llc_fkij(finy,nx,fc2,1,ix2,jx-1); % <<<<<<<<
fldy=read_llc_fkij(finx,nx,fc2,1,ix2,jx-1); % <<<<<<<<
writebin(foutx,fldx);
writebin(fouty,fldy);
% }}}
% {{{ southwest edge points, no direction
fnx='DXG';
fny='DYG';
finx=[gdir fnx '.data'];
finy=[gdir fny '.data'];
foutx=[fnx suf1];
fouty=[fny suf1];
fldx=read_llc_fkij(finy,nx,fc2,1,ix2,jx-1); % <<<<<<<<
fldy=read_llc_fkij(finx,nx,fc2,1,ix2,jx);
writebin(foutx,fldx);
writebin(fouty,fldy);
% }}}
% }}}
% {{{ create commands for extracting model output fields
startPoint=[int2str((fc2-2)*nx+min(ix2)) ',' int2str(min(jx)) ',1 '];
% }}}
% {{{ get and save regional 2D tracer fields
extent=[int2str(m2) ',' int2str(n) ',1'];
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw'}
    eval(['mkdir ' pin2 fnm{1}])
    eval(['cd ' pin2 fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pin2 fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -a joblist')
end
% }}}
% {{{ get and save regional 3D tracer fields
extent=[int2str(m2) ',' int2str(n) ',' int2str(z)];
for fnm={'Salt','Theta','W'}
    eval(['mkdir ' pin2 fnm{1}])
    eval(['cd ' pin2 fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pin2 fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -a joblist')
end
% }}}
% {{{ get and save regional vector fields
% 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
% there is a -1 index offset -U meridional velocity in faces 4/5
% example below is for face 4
% {{{ oceTAUX
extent=[int2str(m2) ',' int2str(n) ',1'];
eval(['mkdir ' pin2 'oceTAUX'])
eval(['cd ' pin2 'oceTAUX'])
fieldNames=['oceTAUY' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_oceTAUY_/_oceTAUX_/'' joblist');
disp(['cd ' pin2 'oceTAUX'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ U
extent=[int2str(m2) ',' int2str(n) ',' int2str(z)];
eval(['mkdir ' pin2 'U'])
eval(['cd ' pin2 'U'])
fieldNames=['V' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_V_/_U_/'' joblist');
disp(['cd ' pin2 'U'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ oceTAUY
startPoint=[int2str((fc2-2)*nx+min(ix2)) ',' int2str(min(jx)-1) ',1 '];
extent=[int2str(m2) ',' int2str(n) ',1'];
eval(['mkdir ' pin2 'oceTAUY'])
eval(['cd ' pin2 'oceTAUY'])
fieldNames=['oceTAUX' ' '];
system([extract '-n ' timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_oceTAUX_/_oceTAUY_/'' joblist');
system('sed -i ''s/_Neg//'' joblist');
disp(['cd ' pin2 'oceTAUY'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ V
extent=[int2str(m2) ',' int2str(n) ',' int2str(z)];
eval(['mkdir ' pin2 'V'])
eval(['cd ' pin2 'V'])
fieldNames=['U' ' '];
system([extract '-n ' timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_U_/_V_/'' joblist');
system('sed -i ''s/_Neg//'' joblist');
disp(['cd ' pin2 'V'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% }}}
% }}}

% {{{ combine face 2 and face 4
% {{{ grid
eval(['mkdir ' pout 'grid/'])
for fnm={'Depth','RAC','XC','YC','XG','YG', ...
         'RAZ','DXC','DYC','DXV','DYU','DXG','DYG'}
    fn1=[pin1 'grid/' fnm{1} '_61x222'];
    fn2=[pin2 'grid/' fnm{1} '_131x222'];
    fno=[pout 'grid/' fnm{1} '_192x222'];
    fl1=readbin(fn1,[61 222]);
    fl2=readbin(fn2,[131 222]);
    fld=[fl1; fl2];
    writebin(fno,fld);
    clf, mypcolor(fld'); title(fnm{1}), pause
end
fnm={'hFacC'};
fn1=[pin1 'grid/' fnm{1} '_61x222x87'];
fn2=[pin2 'grid/' fnm{1} '_131x222x87'];
fno=[pout 'grid/' fnm{1} '_192x222x87'];
fl1=readbin(fn1,[61 222 87]);
fl2=readbin(fn2,[131 222 87]);
fld=[fl1; fl2];
writebin(fno,fld);
% }}}
% {{{ 2D fields
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw','oceTAUX','oceTAUY'}
    eval(['mkdir ' pout fnm{1} '/'])
    for ts=mints:144:maxts
        sts=myint2str(ts,10);
        fn1=[pin1 fnm{1} '/' sts '_' fnm{1} '_8580.8890.1_61.222.1'];
        fn2=[pin2 fnm{1} '/' sts '_' fnm{1} '_8641.8890.1_131.222.1'];
        if strcmp(fnm{1},'oceTAUY')
            fn2=[pin2 fnm{1} '/' sts '_' fnm{1} '_8641.8889.1_131.222.1'];
        end
        fno=[pout fnm{1} '/' fnm{1} '_192x222.' ts2dte(ts,25,2011,9,10,30)];
        fl1=readbin(fn1,[m1 n]);
        fl2=readbin(fn2,[m2 n]);
        fld=[fl1; fl2];
        writebin(fno,fld);
    end
end
% }}}
% {{{ 3D fields
for fnm={'Salt','Theta','U','V','W'}
    eval(['mkdir ' pout fnm{1} '/'])
    for ts=mints:144:maxts
        sts=myint2str(ts,10);
        fn1=[pin1 fnm{1} '/' sts '_' fnm{1} '_8580.8890.1_61.222.87'];
        fn2=[pin2 fnm{1} '/' sts '_' fnm{1} '_8641.8890.1_131.222.87'];
        if strcmp(fnm{1},'V')
            fn2=[pin2 fnm{1} '/' sts '_' fnm{1} '_8641.8889.1_131.222.87'];
        end
        fno=[pout fnm{1} '/' fnm{1} '_192x222x87.' ts2dte(ts,25,2011,9,10,30)];
        if ~exist(fn1) | ~exist(fn2)
            break
        end
        d1=dir(fn1); d2=dir(fn2);
        if d1.bytes==4712616 & d2.bytes==10120536
            fl1=readbin(fn1,[m1 n z]);
            fl2=readbin(fn2,[m2 n z]);
            fld=[fl1; fl2];
            writebin(fno,fld);
        else
            break
        end
    end
end
% }}}
% }}}