% extract variables in California Coast
% 130W to 110W, 10N to 40N

% define desired region
nx=4320;
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_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')
close all

% Get and save grid information
pin='~dmenemen/llc_4320/grid/';
pout=['~dmenemen/llc_4320/regions/' 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))];
        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

% Southwest corner (vorticity) points, no direction
fld=zeros(sum(m),n);
for fnm={'XG','YG','RAZ'}
    fin=[pin fnm{1} '.data'];
    fout=[fnm{1} '_' int2str(sum(m)) 'x' int2str(n)];
    for f=1:length(fc)
        switch fc(f)
          case {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)});
          case {4,5}
            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
    end
    writebin(fout,fld);
end

% West edge points, no direction
fnx='DXC';
fny='DYC';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
fldx=zeros(sum(m),n);
fldy=zeros(sum(m),n);
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        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
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest corner (vorticity) points, no direction
fnx='DXV';
fny='DYU';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        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
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% Southwest edge points, no direction
fnx='DXG';
fny='DYG';
finx=[pin fnx '.data'];
finy=[pin fny '.data'];
foutx=[fnx '_' int2str(sum(m)) 'x' int2str(n)];
fouty=[fny '_' int2str(sum(m)) 'x' int2str(n)];
for f=1:length(fc)
    switch fc(f)
      case {1,2}
        fldx((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finx,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
        fldy((sum(m(1:f))+1):sum(m(1:(f+1))),:) = ...
            read_llc_fkij(finy,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
      case {4,5}
        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
end
writebin(foutx,fldx);
writebin(fouty,fldy);

% get and save regional W/S/T
pn='~dmenemen/llc_4320/MITgcm/run';
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
for fnm={'W','Salt','Theta'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=10368:144:1365984, mydisp(ts)
        if ts<485568
            pin=[pn '/'];
        else
            pin=[pn '_485568/'];
        end
        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 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_4320/MITgcm/run';
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
eval(['mkdir ' pout 'U'])
eval(['mkdir ' pout 'V'])
eval(['cd ' pout])
for ts=10368:144:1365984, mydisp(ts)
    if ts<485568
        pin=[pn '/'];
    else
        pin=[pn '_485568/'];
    end
    finu=[pin 'U.' myint2str(ts,10) '.data'];
    finv=[pin 'V.' myint2str(ts,10) '.data'];
    dy=ts2dte(ts,25,2011,9,10,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 KPPhbl, ETAN, PhiBot
pn='~dmenemen/llc_4320/MITgcm/run';
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
for fnm={'KPPhbl','Eta','PhiBot'}
    eval(['mkdir ' pout fnm{1}])
    eval(['cd ' pout fnm{1}])
    for ts=10368:144:1365984, mydisp(ts)
        if ts<485568
            pin=[pn '/'];
        else
            pin=[pn '_485568/'];
        end
        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 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 fields of oceTAUX and oceTAUY
% note that zonal velocity is oceTAUX in faces 1/2 and oceTAUY in faces 4/5
% and meridional velocity is oceTAUY in faces 1/2 and -oceTAUX in faces 4/5
pn='~dmenemen/llc_4320/MITgcm/run';
pout=['~dmenemen/llc_4320/regions/' region_name '/'];
eval(['mkdir ' pout 'oceTAUX'])
eval(['mkdir ' pout 'oceTAUY'])
eval(['cd ' pout])
for ts=10368:144:1365984, mydisp(ts)
    if ts<485568
        pin=[pn '/'];
    else
        pin=[pn '_485568/'];
    end
    finu=[pin 'oceTAUX.' myint2str(ts,10) '.data'];
    finv=[pin 'oceTAUY.' myint2str(ts,10) '.data'];
    dy=ts2dte(ts,25,2011,9,10,30);
    foutu=['oceTAUX/oceTAUX_' int2str(sum(m)) 'x' int2str(n) '.' dy];
    foutv=['oceTAUY/oceTAUY_' int2str(sum(m)) 'x' int2str(n) '.' dy];
    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),1,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),1,ix{fc(f)},jx{fc(f)});
            fldv((sum(m(1:f))+1):sum(m(1:(f+1))),:) = - ...
                read_llc_fkij(finu,nx,fc(f),1,ix{fc(f)},jx{fc(f)}-1);
        end
        writebin(foutu,fldu);
        writebin(foutv,fldv);
    end
end