% Labrador Sea domain for Zhiyu Liu and Weiwei
% oceFWflx, oceQnet, oceQsw, oceTAUX, oceTAUY,
% Eta, KPPhbl, PhiBot, Salt, Theta, U, V, W
% extract 01-Nov-2011 to 31-May-2012
% lats 51N to 66N, lons 64W to 36W
% (example extraction on face 5+1)
% extract face 1 and face 5 separately
% then use cat to concatenate

% {{{ define desired region
nx=4320;
prec='real*4';
region_name='LabSea';
minlon=360-64;
maxlon=360-36;
minlat=51;
maxlat=66;
mints=dte2ts('01-Nov-2011',25,2011,9,10);
maxts=dte2ts('31-May-2012',25,2011,9,10);
timesteps=[int2str(mints) '-' int2str(maxts) ' '];
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
pin1=[pout 'fc5/'];
pin2=[pout 'fc1/'];
eval(['mkdir ' pin1])
eval(['mkdir ' pin2])
% }}}

% {{{ 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);
FC=[5 1];
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;
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*
% }}}

% {{{ extract face 5
% {{{ Get and save grid information
eval(['mkdir ' pin1 'grid/'])
eval(['cd ' pin1 'grid/'])
suf1=['_' int2str(m1) '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,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
% }}}
% {{{ 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,fc1,1,ix1,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,fc1,1,ix1,jx);
fldy=read_llc_fkij(finx,nx,fc1,1,ix1,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,fc1,1,ix1,jx-1); % <<<<<<<<
fldy=read_llc_fkij(finx,nx,fc1,1,ix1,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,fc1,1,ix1,jx-1); % <<<<<<<<
fldy=read_llc_fkij(finx,nx,fc1,1,ix1,jx);
writebin(foutx,fldx);
writebin(fouty,fldy);
% }}}
% }}}
% {{{ create commands for extracting model output fields
startPoint=[int2str((fc1-2)*nx+min(ix1)) ',' int2str(min(jx)) ',1 '];
extract='/home4/bcnelson/MITgcm/extract/latest/extract4320 -g ';
% }}}
% {{{ get and save regional 2D tracer fields
extent=[int2str(m1) ',' int2str(n) ',1'];
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw', ...
         'SIarea','SIheff','SIhsalt','SIhsnow','oceSflux'}
    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 tracer fields
extent=[int2str(m1) ',' int2str(n) ',' int2str(z)];
for fnm={'Salt','Theta','W'}
    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 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 5
% {{{ oceTAUX
extent=[int2str(m1) ',' int2str(n) ',1'];
eval(['mkdir ' pin1 'oceTAUX'])
eval(['cd ' pin1 'oceTAUX'])
fieldNames=['oceTAUY' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_oceTAUY_/_oceTAUX_/'' joblist');
disp(['cd ' pin1 'oceTAUX'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ SIuice
extent=[int2str(m1) ',' int2str(n) ',1'];
eval(['mkdir ' pin1 'SIuice'])
eval(['cd ' pin1 'SIuice'])
fieldNames=['SIvice' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_SIvice_/_SIuice_/'' joblist');
disp(['cd ' pin1 'SIuice'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ U
extent=[int2str(m1) ',' int2str(n) ',' int2str(z)];
eval(['mkdir ' pin1 'U'])
eval(['cd ' pin1 'U'])
fieldNames=['V' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_V_/_U_/'' joblist');
disp(['cd ' pin1 'U'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ oceTAUY
startPoint=[int2str((fc1-2)*nx+min(ix1)) ',' int2str(min(jx)-1) ',1 '];
extent=[int2str(m1) ',' int2str(n) ',1'];
eval(['mkdir ' pin1 'oceTAUY'])
eval(['cd ' pin1 '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 ' pin1 'oceTAUY'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ SIvice
startPoint=[int2str((fc1-2)*nx+min(ix1)) ',' int2str(min(jx)-1) ',1 '];
extent=[int2str(m1) ',' int2str(n) ',1'];
eval(['mkdir ' pin1 'SIvice'])
eval(['cd ' pin1 'SIvice'])
fieldNames=['SIuice' ' '];
system([extract '-n ' timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_SIuice_/_SIvice_/'' joblist');
system('sed -i ''s/_Neg//'' joblist');
disp(['cd ' pin1 'SIvice'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ V
extent=[int2str(m1) ',' int2str(n) ',' int2str(z)];
eval(['mkdir ' pin1 'V'])
eval(['cd ' pin1 '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 ' pin1 'V'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% }}}
% }}}

% {{{ extract face 1
% {{{ get and save grid information
eval(['mkdir ' pin2 'grid/'])
eval(['cd ' pin2 'grid/'])
suf1=['_' int2str(m2) '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,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
% }}}
% {{{ create commands for extracting model output fields
startPoint=[int2str(min(ix2)) ',' int2str(min(jx)) ',1 '];
extract='/home4/bcnelson/MITgcm/extract/latest/extract4320 -g ';
% }}}
% {{{ get and save regional 2D fields
extent=[int2str(m2) ',' int2str(n) ',1'];
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw','oceTAUX','oceTAUY', ...
         'SIarea','SIheff','SIhsalt','SIhsnow','SIuice','SIvice','oceSflux'}
    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 fields
extent=[int2str(m2) ',' int2str(n) ',' int2str(z)];
for fnm={'Salt','Theta','U','V','W'}
    eval(['mkdir ' pin2 fnm{1}])
    eval(['cd ' pin2 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
% }}}
% }}}

% {{{ combine face 5 and face 1
% {{{ 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} '_1304x1600'];
    fn2=[pin2 'grid/' fnm{1} '_124x1600'];
    fno=[pout 'grid/' fnm{1} '_1428x1600'];
    fl1=readbin(fn1,[1304 1600]);
    fl2=readbin(fn2,[124 1600]);
    fld=[fl1; fl2];
    writebin(fno,fld);
    clf, mypcolor(fld'); title(fnm{1}), pause
end
fnm={'hFacC'};
fn1=[pin1 'grid/' fnm{1} '_1304x1600x84'];
fn2=[pin2 'grid/' fnm{1} '_124x1600x84'];
fno=[pout 'grid/' fnm{1} '_1428x1600x84'];
fl1=readbin(fn1,[1304 1600 84]);
fl2=readbin(fn2,[124 1600 84]);
fld=[fl1; fl2];
writebin(fno,fld);
% }}}
% {{{ 2D fields
for fnm={'Eta','PhiBot','KPPhbl','oceFWflx','oceQnet','oceQsw','oceTAUX','oceTAUY', ...
         'SIarea','SIheff','SIhsalt','SIhsnow','SIuice','SIvice','oceSflux'}
    eval(['mkdir ' pout fnm{1} '/'])
    for ts=mints:144:maxts
        sts=myint2str(ts,10);
        fn1=[pin1 fnm{1} '/' sts '_' fnm{1} '_15977.11006.1_1304.1600.1'];
        fn2=[pin2 fnm{1} '/' sts '_' fnm{1} '_1.11006.1_124.1600.1'];
        if strcmp(fnm{1},'oceTAUY') | strcmp(fnm{1},'SIvice')
           fn1=[pin1 fnm{1} '/' sts '_' fnm{1} '_15977.11005.1_1304.1600.1'];
        end
        fno=[pout fnm{1} '/' fnm{1} '_1428x1600.' 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} '_15977.11006.1_1304.1600.84'];
        fn2=[pin2 fnm{1} '/' sts '_' fnm{1} '_1.11006.1_124.1600.84'];
        if strcmp(fnm{1},'V')
            fn1=[pin1 fnm{1} '/' sts '_' fnm{1} '_15977.11005.1_1304.1600.84'];
        end
        fno=[pout fnm{1} '/' fnm{1} '_1428x1600x84.' ts2dte(ts,25,2011,9,10,30)];
        if ~exist(fn1) | ~exist(fn2)
            break
        end
        d1=dir(fn1); d2=dir(fn2);
        if d1.bytes==701030400 & d2.bytes==66662400
            fl1=readbin(fn1,[m1 n z]);
            fl2=readbin(fn2,[m2 n z]);
            fld=[fl1; fl2];
            writebin(fno,fld);
        else
            break
        end
    end
end
% }}}
% }}}