% extract variables in California Coast
% 90-130E, -5 to 25N

% define desired region
nx=2160;
prec='real*4';
region_name='CalCoast';
minlat=10;
maxlat=40;
minlon=-130;
maxlon=-110;
kx=1:88;

% extract indices for desired region
gdir='~dmenemen/llc_2160/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')

% get and save 2D grid information
close all
pin='~dmenemen/llc_2160/grid/';
pout=['~dmenemen/llc_2160/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'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(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

% get and save 3D grid information
for fnm={'hFacC','hFacS','hFacW'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx))];
    for k=1:length(kx); mydisp(k)
        for f=1:length(fc)
            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(fout,fld,1,'real*4',k-1);
    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
pn='~dmenemen/llc_2160/MITgcm/run';
pout=['~dmenemen/llc_2160/regions/' region_name '/'];
eval(['mkdir ' pout 'U'])
eval(['mkdir ' pout 'V'])
eval(['cd ' pout])
for ts=92160:80:1586400; disp(ts)
    if (ts<1198080)
        pin=[pn '_day49_624/'];
    else
        pin=[pn '/'];
    end
    finu=[pin 'U.' myint2str(ts,10) '.data'];
    finv=[pin 'V.' myint2str(ts,10) '.data'];
    dy=ts2dte(ts,45,2011,1,17,30);
    foutu=['U/U_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    foutv=['V/V_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
    for k=1:length(kx); mydisp(k)
        for f=1:length(fc)
            switch fc(f)
              case {1,2}
                fldu((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
              case {4,5}
                fldu((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
                    read_llc_fkij(finv,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)});
                fldv((sum(m(1:f))+1):sum(m(1:(f+1))),:) = - ...
                    read_llc_fkij(finu,nx,fc(f),kx(k),ix{fc(f)},jx{fc(f)}-1);
            end
        end
        writebin(foutu,fldu,1,'real*4',k-1);
        writebin(foutv,fldv,1,'real*4',k-1);
    end
end

% get and save S,T,W
for fnm={'Salt','Theta','W'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=92160:80:1586400; disp(ts)
        if (ts<1198080)
            pin=[pn '_day49_624/'];
        else
            pin=[pn '/'];
        end
        fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
        dy=ts2dte(ts,45,2011,1,17,30);
        fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) 'x' int2str(length(kx)) '.' dy];
        for k=1:length(kx); mydisp(k)
            for f=1:length(fc)
                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(fout,fld,1,'real*4',k-1);
        end
    end
end

% get and save KPPhbl,oceTAUX,oceTAUY,Eta,PhiBot
pn='~dmenemen/llc_2160/MITgcm/run';
pout=['~dmenemen/llc_2160/regions/' region_name '/'];
for fnm={'KPPhbl','oceTAUX','oceTAUY','Eta','PhiBot'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=92160:80:1586400; disp(ts)
        if (ts<1198080)
            pin=[pn '_day49_624/'];
        else
            pin=[pn '/'];
        end
        fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
        dy=ts2dte(ts,45,2011,1,17,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