% LLC cutout aligned with latitude-longitude lines
% surface U/V for surface trajectory computations (Tom/Dustin/Kyle)
% PhiBot, bottom U/V, oceTAUX/Y for bottom drag computations (Ali/Bowen)
% extract 13-Sep-2011 to 15-Nov-2012
% lats 70S to 57N
% (example extraction on face 1+2+4+5)
% extract faces 1-2 and 4-5 separately
% then use matlab to concatenate

% {{{ define desired region
nx=4320;
prec='real*4';
region_name='latlon';
minlat=-70;
maxlat=57;
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 '/'];
pn1=[pout 'pn1/'];
pn2=[pout 'pn2/'];
eval(['mkdir ' pn1])
eval(['mkdir ' pn2])
extract='/home4/bcnelson/MITgcm/extract/latest/extract4320 -g ';
% }}}

% {{{ 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);
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')
kx=1:90;
close all
% }}}

% {{{ extract faces 1 and 2
% {{{ get and save grid information
eval(['mkdir ' pn1 'grid/'])
eval(['cd ' pn1 'grid/'])
fld=zeros(nx*4,nx*2);
for fnm={'Depth','RAC','XC','YC','hFacC','XG','YG', ...
         'RAZ','DXC','DYC','DXV','DYU','DXG','DYG', ...
         'AngleCS','AngleSN','DXF','DYF','IDX','LandMask', ...
         'RAS','RAW','U2zonDir','V2zonDir','hFacS','hFacW'}    
    fin=[gdir fnm{1} '.data'];
    switch fnm{1}
      case{'hFacC','hFacS','hFacW'}
        for k=1:length(kx); mydisp(k)
            for f=1:2
                fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(fin,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
            end
            writebin(fnm{1},fld,1,prec,k-1);
        end
      otherwise
        for f=1:2
            fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        end
        writebin(fnm{1},fld);
    end
end
% }}}
% {{{ create commands for extracting model output fields
startPoint=[int2str((fc(1)-1)*nx+min(ix{fc(1)})) ',' int2str(min(jx{fc(1)})) ',1 '];
% }}}
% {{{ get and save regional 2D fields
extent=[int2str(sum(m(2:3))) ',' int2str(n) ',1'];
for fnm={'PhiBot','oceTAUX','oceTAUY'}
    eval(['mkdir ' pn1 fnm{1}])
    eval(['cd ' pn1 fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pn1 fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end
% }}}
% {{{ get and save regional 3D fields
extent=[int2str(sum(m(2:3))) ',' int2str(n) ',1'];
for fnm={'U','V'}
    eval(['mkdir ' pn1 fnm{1}])
    eval(['cd ' pn1 fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pn1 fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end
% }}}
% }}}

% {{{ extract faces 4 and 5
% {{{ get and save grid information

eval(['mkdir ' pn2 'grid/'])
eval(['cd ' pn2 'grid/'])
fld=zeros(nx*2);
fldx=zeros(nx*2);
fldy=zeros(nx*2);
fc=[4 5];
% {{{ Grid cell center
for fnm={'Depth','RAC','XC','YC','hFacC'}
    fin=[gdir fnm{1} '.data'];
    switch fnm{1}
      case{'hFacC'}
        for k=1:length(kx); mydisp(k)
            for f=1:2
                fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(fin,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
            end
            writebin(fnm{1},fld,1,prec,k-1);
        end
      otherwise
        for f=1:2
            fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        end
        writebin(fnm{1},fld);
    end
end
% }}}
% {{{ Southwest corner (vorticity) points, no direction
for fnm={'XG','YG','RAZ'}
    fin=[gdir fnm{1} '.data'];
    for f=1:2
        fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
    end
    writebin(fnm{1},fld);
end
% }}}
% {{{ West edge points, no direction
fnx='DXC';
fny='DYC';
finx=[gdir fnx '.data'];
finy=[gdir fny '.data'];
for f=1:2
    fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
    fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
end
writebin(fnx,fldx);
writebin(fny,fldy);
% }}}
% {{{ Southwest corner (vorticity) points, no direction
fnx='DXV';
fny='DYU';
finx=[gdir fnx '.data'];
finy=[gdir fny '.data'];
for f=1:2
    fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
    fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
end
writebin(fnx,fldx);
writebin(fny,fldy);
% }}}
% {{{ Southwest edge points, no direction
fnx='DXG';
fny='DYG';
finx=[gdir fnx '.data'];
finy=[gdir fny '.data'];
for f=1:2
    fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1); % <<<<<<<<
    fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
        read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
end
writebin(fnx,fldx);
writebin(fny,fldy);
% }}}

% }}}
% {{{ create commands for extracting model output fields
startPoint=[int2str((fc(3)-2)*nx+min(ix{fc(3)})) ',' int2str(min(jx{fc(3)})) ',1 '];
% }}}
% {{{ get and save regional 2D tracer fields
extent=[int2str(sum(m(4:5))) ',' int2str(n) ',1'];
for fnm={'PhiBot'}
    eval(['mkdir ' pn2 fnm{1}])
    eval(['cd ' pn2 fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pn2 fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
end
% }}}
% {{{ get and save regional 3D tracer fields
extent=[int2str(sum(m(4:5))) ',' int2str(n) ',1'];
for fnm={}
    eval(['mkdir ' pn2 fnm{1}])
    eval(['cd ' pn2 fnm{1}])
    fieldNames=[fnm{1} ' '];
    system([extract timesteps  fieldNames  startPoint  extent '  > joblist']);
    disp(['cd ' pn2 fnm{1}])
    disp('parallel --slf $PBS_NODEFILE -j2 -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
startPoint=[int2str((fc(3)-2)*nx+min(ix{fc(3)})) ',' int2str(min(jx{fc(3)})) ',1 '];
% {{{ oceTAUX
extent=[int2str(sum(m(4:5))) ',' int2str(n) ',1'];
eval(['mkdir ' pn2 'oceTAUX'])
eval(['cd ' pn2 'oceTAUX'])
fieldNames=['oceTAUY' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_oceTAUY_/_oceTAUX_/'' joblist');
disp(['cd ' pn2 'oceTAUX'])
disp('parallel --slf $PBS_NODEFILE -j2 -a joblist')
% }}}
% {{{ U
extent=[int2str(sum(m(4:5))) ',' int2str(n) ',1'];
eval(['mkdir ' pn2 'U'])
eval(['cd ' pn2 'U'])
fieldNames=['V' ' '];
system([extract timesteps  fieldNames  startPoint  extent ' > joblist']);
system('sed -i ''s/_V_/_U_/'' joblist');
disp(['cd ' pn2 'U'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
startPoint=[int2str((fc(3)-2)*nx+min(ix{fc(3)})) ',' int2str(min(jx{fc(3)})-1) ',1 '];
% {{{ oceTAUY
extent=[int2str(sum(m(4:5))) ',' int2str(n) ',1'];
eval(['mkdir ' pn2 'oceTAUY'])
eval(['cd ' pn2 '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 ' pn2 'oceTAUY'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% {{{ V
extent=[int2str(sum(m(4:5))) ',' int2str(n) ',1'];
eval(['mkdir ' pn2 'V'])
eval(['cd ' pn2 '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 ' pn2 'V'])
disp('parallel --slf $PBS_NODEFILE -a joblist')
% }}}
% }}}
% }}}

% {{{ combine face 1-2 and face 4-5
% {{{ grid
eval(['mkdir ' pout 'grid/']) 
for fnm={'Depth','RAC','XC','YC','XG','YG', ...
         'RAZ','DXC','DYC','DXV','DYU','DXG','DYG'}
    fn1=[pn1 'grid/' fnm{1}];
    fn2=[pn2 'grid/' fnm{1}];
    fno=[pout 'grid/' fnm{1} '_17280x8640'];
    fl1=readbin(fn1,[nx*2 nx*2]);
    fl2=readbin(fn2,[nx*2 nx*2]);
    fld=[fl1; fl2];
    writebin(fno,fld);
    clf, mypcolor(fld'); title(fnm{1}), pause
end
fnm={'hFacC'};
fn1=[pn1 'grid/' fnm{1}];
fn2=[pn2 'grid/' fnm{1}];
fno=[pout 'grid/' fnm{1} '_17280x8640x90'];
for k=1:90
    fl1=readbin(fn1,[nx*2 nx*2],1,prec,k-1);
    fl2=readbin(fn2,[nx*2 nx*2],1,prec,k-1);
    fld=[fl1; fl2];
    writebin(fno,fld,1,prec,k-1);
    clf, mypcolor(fld'); title(k), pause(.01)
end
% }}}
% {{{ 2D fields
for fnm={'PhiBot','oceTAUX','oceTAUY'}
    eval(['mkdir ' pout fnm{1} '/'])
    for ts=mints:144:maxts
        sts=myint2str(ts,10);
        fn1=[pn1 fnm{1} '/' sts '_' fnm{1} '_1.2881.1_8640.8640.1'];
        fn2=[pn2 fnm{1} '/' sts '_' fnm{1} '_8641.2881.1_8640.8640.1'];
        if strcmp(fnm{1},'oceTAUY')
            fn2=[pn2 fnm{1} '/' sts '_' fnm{1} '_8641.2880.1_8640.8640.1'];
        end
        fno=[pout fnm{1} '/' fnm{1} '_17280x8640x1.' ts2dte(ts,25,2011,9,10,30)];
        fl1=readbin(fn1,[nx*2 nx*2]);
        fl2=readbin(fn2,[nx*2 nx*2]);
        fld=[fl1; fl2];
        writebin(fno,fld);
    end
end
% }}}
% {{{ 3D fields
for fnm={'U','V'}
    eval(['mkdir ' pout fnm{1} '/'])
    for ts=mints:144:maxts
        sts=myint2str(ts,10);
        fn1=[pn1 fnm{1} '/' sts '_' fnm{1} '_1.2881.1_8640.8640.1'];
        fn2=[pn2 fnm{1} '/' sts '_' fnm{1} '_8641.2881.1_8640.8640.1'];
        if strcmp(fnm{1},'V')
            fn2=[pn2 fnm{1} '/' sts '_' fnm{1} '_8641.2880.1_8640.8640.1'];
        end
        fno=[pout fnm{1} '/' fnm{1} '_17280x8640x1.' ts2dte(ts,25,2011,9,10,30)];
        fl1=readbin(fn1,[nx*2 nx*2]);
        fl2=readbin(fn2,[nx*2 nx*2]);
        fld=[fl1; fl2];
        writebin(fno,fld);
    end
end
% }}}
% }}}