% daily-averaged for one year
% define desired region
more off
nx=2160;
region_name='Weddell/Region1';
minlat=-78;
maxlat=-60;
minlon=298;
maxlon=345;
kx=1:90;

% extract indices for desired region
prec='real*4';
gdir=['~/llc_' int2str(nx) '/grid/'];
fin=[gdir 'Depth.data'];
[tmp fc ix jx]=quikread_llc(fin,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);

% put face 5 output first
fc=[5 1];
m(1)=0;
for f=1:length(fc)
    m(f+1)=length(ix{fc(f)});
end

% make sure jx is same for all faces
jx=min([jx{1} jx{5}]):max([jx{1} jx{5}]);
n=length(jx);

% map field on array and plot
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);
end
quikpcolor(fld')
close all

% Get and save grid information
pout=['~dmenemen/llc_' int2str(nx) '/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
for fnm={'AngleCS','AngleSN','DXC','DXG','DYC','DYG','Depth', ...
         'RAC','RAS','RAW','RAZ','U2zonDir','V2zonDir', ...
         'XC','XG','YC','YG','hFacC','hFacS','hFacW'}
    fin=[gdir fnm{1} '.data'];
    disp(fin)
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    switch fnm{1}
      case{'hFacC','hFacS','hFacW'}
        fld=zeros(sum(m),n,length(kx));
        for f=1:length(fc)
            fld((sum(m(1:f))+1):sum(m(1:(f+1))),:,:) = ...
                read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx);
        end
        fout=[fout 'x' int2str(length(kx))];
      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);
        end
    end
    writebin(fout,fld);
end

% get and save regional W/S/T
pin=['~dmenemen/llc_' int2str(nx) '/MITgcm/run_day49_624/'];
pout=['~dmenemen/llc_' int2str(nx) '/regions/' region_name '/'];
for fnm={'W','Salt','Theta'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for dy=datenum(2011,9,13):datenum(2012,9,12)
        fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) ...
              '.' datestr(dy,29)];
        for k=1:length(kx)
            fld=zeros(sum(m),n);
            disp([fout ' ' int2str(k) ' ' datestr(now)])
            for hr=1:24
                ts=dte2ts(datestr(dy+hr/24),45,2011,1,17);
                fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
                for f=1:length(fc)
                    fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                        fld((sum(m(1:f))+1):sum(m(1:(f+1))),:) + ...
                        read_llc_fkij(fin,nx,fc(f),kx(k),ix{fc(f)},jx) / 24;
                end
            end
            writebin(fout,fld,1,'real*4',k-1);
        end
    end
end

% get and save regional U/V
clear fld
eval(['mkdir ' pout 'U'])
eval(['mkdir ' pout 'V'])
eval(['cd ' pout])
for dy=datenum(2011,9,13):datenum(2012,9,12)
    foutu=['U/U_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) ...
           '.' datestr(dy,29)];
    foutv=['V/V_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) ...
           '.' datestr(dy,29)];
    for k=1:length(kx)
        u=zeros(sum(m),n);
        v=zeros(sum(m),n);
        disp([foutu ' ' int2str(k) ' ' datestr(now)])
        for hr=1:24
            ts=dte2ts(datestr(dy+hr/24),45,2011,1,17);
            fu=[pin 'U.' myint2str(ts,10) '.data'];
            fv=[pin 'V.' myint2str(ts,10) '.data'];
            f=2;
            u((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                u((sum(m(1:f))+1):sum(m(1:(f+1))),:) + ...
                read_llc_fkij(fu,nx,fc(f),kx(k),ix{fc(f)},jx) / 24;
            v((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                v((sum(m(1:f))+1):sum(m(1:(f+1))),:) + ...
                read_llc_fkij(fv,nx,fc(f),kx(k),ix{fc(f)},jx) / 24;
            f=1;
            u((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                u((sum(m(1:f))+1):sum(m(1:(f+1))),:) + ...
                read_llc_fkij(fv,nx,fc(f),kx(k),ix{fc(f)},jx) / 24;
            v((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                v((sum(m(1:f))+1):sum(m(1:(f+1))),:) - ...
                read_llc_fkij(fu,nx,fc(f),kx(k),ix{fc(f)},jx-1) / 24;
        end
        writebin(foutu,u,1,'real*4',k-1);
        writebin(foutv,v,1,'real*4',k-1);
    end
end