% define desired region
region_name='Chile';
minlat=-35;
maxlat=0;
minlon=-90;
maxlon=-70;

% extract indices for desired region
nx=2160;
prec='real*4';
gdir='~dmenemen/llc_2160/grid/';
fnam=[gdir 'Depth.data'];
[fld fc ix jx] = ...
    quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);
quikpcolor(fld')
close all

% get and save regional fields
pn='~dmenemen/llc_2160/MITgcm/run';
pout=['~dmenemen/llc_2160/regions/' region_name '/'];
kx=1:90;
for fnm={'Salt','Theta','U','V','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
        switch fnm{1}
          case 'U'
            fin=[pin 'V.' myint2str(ts,10) '.data'];
          case 'V'
            fin=[pin 'U.' myint2str(ts,10) '.data'];
          otherwise
            fin=[pin fnm{1} '.' myint2str(ts,10) '.data'];
        end
        dy=ts2dte(ts,45,2011,1,17,30);
        fout=[fnm{1} '_' int2str(length(ix)) 'x' int2str(length(jx)) ...
              'x' int2str(length(kx)) '.' dy];
        for k=1:length(kx); mydisp(k)
            if strcmp(fnm{1},'V')
                fld=-read_llc_fkij(fin,nx,fc,kx(k),ix,jx-1);
            else
                fld=read_llc_fkij(fin,nx,fc,kx(k),ix,jx);
            end
            writebin(fout,fld,1,'real*4',k-1);
        end
    end
end

% daily averages
for fnm={'Salt','Theta','U','V','W'}
    eval(['mkdir ' pout fnm{1} '_daily'])
    eval(['cd ' pout fnm{1} '_daily'])
    for ts=92160:(80*24):1586400
        dy=ts2dte(ts,45,2011,1,17,29); disp(dy);
        fout=[fnm{1} '_' int2str(length(ix)) 'x' int2str(length(jx)) ...
              'x' int2str(length(kx)) '.' dy];
        fld=0; n=0;
        for h=0:23, mydisp(h)
            dy=ts2dte(ts+h*80,45,2011,1,17,30);
            fin=[pout fnm{1} '/' fnm{1} '_' int2str(length(ix)) 'x' ...
                 int2str(length(jx)) 'x' int2str(length(kx)) '.' dy];
            if exist(fin)
                fld=fld+readbin(fin,[length(ix) length(jx) length(kx)]);
                n=n+1;
            else
                break
            end
        end
        if n==24
            writebin(fout,fld/n);
        else
            break
        end
    end
end