% extract 15-Oct-2011 to 15-Nov-2012
% 4S-4N, 216E-234E

% define desired region
nx=4320;
prec='real*4';
region_name='Box11';
minlat=-4;
maxlat= 4;
minlon=216;
maxlon=234;
mints=dte2ts('15-Oct-2011',25,2011,9,10);
maxts=dte2ts('15-Nov-2012',25,2011,9,10);

% 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')
load ~dmenemen/llc_4320/grid/thk90.mat
bot90=dpt90+thk90/2; bot90(end)=bot90(end)+1;
kx=1:min(find(bot90>mmax(fld)));

% Get and save grid information
close all
pin='~dmenemen/llc_4320/grid/';
pout=['~dmenemen/llc_4320/regions/Boxes/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
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'}
        fout=[fout 'x' int2str(length(kx))];
        fldd=zeros(sum(m),n,length(kx));
        for f=1:length(fc)
            fldd((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
                read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx{fc(f)});
        end
        writebin(fout,fldd);
      otherwise
        fld=zeros(sum(m),n);
        for f=1:length(fc)
            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(fout,fld);
    end
end

% get and save regional Eta and PhiBot
pin='~dmenemen/llc_4320/MITgcm/run';
pout=['~dmenemen/llc_4320/regions/Boxes/' region_name '/'];
for fnm={'Eta','PhiBot'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=mints:144:maxts, mydisp(ts)
        if ts<=485568
            fin=[pin '/' fnm{1} '.' myint2str(ts,10) '.data'];
        else
            fin=[pin '_485568/' fnm{1} '.' myint2str(ts,10) '.data'];
        end
        dy=ts2dte(ts,25,2011,9,10,30);
        fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) '.' dy];
        for f=1:length(fc)
            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(fout,fld);
    end
end

% get and save regional S/T

for fnm={'Salt'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=mints:144:maxts, mydisp(ts)
        if ts<=485568
            fin=[pin '/' fnm{1} '.' myint2str(ts,10) '.data'];
        else
            fin=[pin '_485568/' fnm{1} '.' myint2str(ts,10) '.data'];
        end
        dy=ts2dte(ts,25,2011,9,10,30);
        fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
        for f=1:length(fc)
            fldd((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
                read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx{fc(f)});
        end
        writebin(fout,fldd);
    end
end

for fnm={'Theta'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=mints:144:maxts, mydisp(ts)
        if ts<=485568
            fin=[pin '/' fnm{1} '.' myint2str(ts,10) '.data'];
        else
            fin=[pin '_485568/' fnm{1} '.' myint2str(ts,10) '.data'];
        end
        dy=ts2dte(ts,25,2011,9,10,30);
        fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
        for f=1:length(fc)
            fldd((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
                read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx{fc(f)});
        end
        writebin(fout,fldd);
    end
end

% 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(['cd ' pout])
for ts=mints:144:maxts, mydisp(ts)
    if ts<485568
        finv=[pin '/V.' myint2str(ts,10) '.data'];
    else
        finv=[pin '_485568/V.' myint2str(ts,10) '.data'];
    end
    dy=ts2dte(ts,25,2011,9,10,30);
    foutu=['U/U_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    for f=1:length(fc)
        fldd((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
            read_llc_fkij(finv,nx,fc(f),kx,ix{fc(f)},jx{fc(f)});
    end
    writebin(foutu,fldd);
end

eval(['mkdir ' pout 'V'])
eval(['cd ' pout])
for ts=mints:144:maxts, mydisp(ts)
    if ts<485568
        finu=[pin '/U.' myint2str(ts,10) '.data'];
    else
        finu=[pin '_485568/U.' myint2str(ts,10) '.data'];
    end
    dy=ts2dte(ts,25,2011,9,10,30);
    foutv=['V/V_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    for f=1:length(fc)
        fldd((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
            read_llc_fkij(finu,nx,fc(f),kx,ix{fc(f)},jx{fc(f)}-1);
    end
    writebin(foutv,fldd);
end