% extract Eta/S/T/U/V/W in Kuroshio Region
% 140-160E, 30-40N

% define desired region
nx=4320;
prec='real*4';
region_name='Kuroshio3';
minlat=30;
maxlat=40;
minlon=140;
maxlon=160;
kx=1:90;

% extract indices for desired region
gdir='/nobackupp2/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')
colorbar

% get and save grid information
close all
pin='/nobackupp2/dmenemen/llc_4320/grid/';
pout=['/nobackupp2/dmenemen/llc_4320/regions/' region_name '/grid/'];
eval(['mkdir ' pout])
eval(['cd ' pout])
for fnm={'Depth','RAC','XC','XG','YC','YG'}
    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))];
        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{fc(f)});
        end
      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
    end
    writebin(fout,fld);
end
for fnm={'hFacC'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    for k=1:length(kx); mydisp(k)
    switch fnm{1}
      case{'hFacC'}
        fout=[fout 'x' int2str(length(kx))];
        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{fc(f)});
        end
      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
    end
    writebin(fout,fld);
end
end

% get and save regional 2D fields
pin='~dmenemen/llc_4320/MITgcm/run/';
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
kx=1;
fld=zeros(sum(m),n);
for fnm={'Eta','PhiBot'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=10368:144:485568, disp(ts)
        fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
        dy=ts2dte(ts,25,2011,9,10,30);
        fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n) '.' 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 regional 3D fields
pin='~dmenemen/llc_4320/MITgcm/run/';
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
kx=1:90;
fld=zeros(sum(m),n);
for fnm={'Salt','Theta'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=10368:144:485568, disp(ts)
        fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
        dy=ts2dte(ts,25,2011,9,10,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 bottom temperature
region_name='Kuroshio';
pin=['~dmenemen/llc_4320/regions/' region_name '/'];
hFacC=readbin([pin 'grid/hFacC_3662x392x90'],[3662,392,90]);
bottom=zeros(3662,392);
for i=1:3662, mydisp(i)
    for j=1:392
        tmp=squeeze(hFacC(i,j,:));
        in=find(tmp);
        if ~isempty(in)
            bottom(i,j)=max(in);
        end
    end
end
eval(['mkdir ' pin 'ThetaBot'])
eval(['cd '  pin 'ThetaBot'])
for ts=10368:144:485568, disp(ts)
    dy=ts2dte(ts,25,2011,9,10,30);
    fin=[pin 'Theta/Theta_3662x392x90.' dy];
    fout=['ThetaBot_3662x392.' dy];
    tmp=readbin(fin,[3662,392,90]);
    theta=0*bottom;
    for i=1:3662
        for j=1:392
            in=bottom(i,j);
            if in
                theta(i,j)=tmp(i,j,in);
            end
        end
    end
    writebin(fout,theta);
end

%%%%%%

cd ~dmenemen/llc_4320/regions/Kuroshio/ThetaBot
ts=10368;
dy=ts2dte(ts,25,2011,9,10,30);
fin=['ThetaBot_3662x392.' dy];
tmp1=readbin(fin,[3662,392]);
for ts=10368:144:485568
    dy=ts2dte(ts,25,2011,9,10,30);
    fin=['ThetaBot_3662x392.' dy];
    tmp=readbin(fin,[3662,392]);
    clf
    mypcolor(tmp'-tmp1');
    caxis([-1 1]/50);
    thincolorbar('horiz')
    title(dy)
    pause(.1)
end