% extract time series of T averaged between
% in following region
nx=2160;
region_name='Umanak';
minlat=70.2425;
maxlat=71.6693;
minlon=-55.5457;
maxlon=-50.7999;
mindpt=200;
maxdpt=500;

% extract vertical indices for desired region
gdir=['~dmenemen/llc_' int2str(nx) '/grid/'];
rc=-readbin([gdir 'RC.data'],[90 1]);
drf=readbin([gdir 'DRF.data'],[90 1]);
kx=find(rc>mindpt&rc<maxdpt);
dpt=rc(kx);
thk=drf(kx);

% extract horizontal indices
prec='real*4';
fnam=[gdir 'Depth.data'];
[fld fc ix jx]=quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);

% get grid area and mask
RAC={};
fin=[gdir 'RAC.data'];
for f=1:length(fc)
    RAC{f}=read_llc_fkij(fin,nx,fc(f),1,ix{fc(f)},jx{fc(f)});
end
hFacC={};
fin=[gdir 'hFacC.data'];
vol=zeros(length(kx),1);
for f=1:length(fc)
    hFacC{f}=read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx{fc(f)});
    for k=1:length(kx)
                        
        
        
        
        hFacC{f}(:,:,k)=hFacC{f}(:,:,k).*RAC{f}*thk(k);
        vol(k)=vol(k)+sum(sum(hFacC{f}(:,:,k)));
    end
end

% get volume-averaged potential temperature
pin=['~dmenemen/llc_' int2str(nx) '/MITgcm/run'];
ts=92160:80:1586400;
Theta=zeros(length(ts),length(kx));
for t=1:length(ts), mydisp(ts(t))
    if (ts(t)<1198080)
        fin=[pin '_day49_624/Theta.' myint2str(ts(t),10) '.data'];
    else 
        fin=[pin '/Theta.' myint2str(ts(t),10) '.data'];
    end        
    for f=1:length(fc)
        tmp=read_llc_fkij(fin,nx,fc(f),kx,ix{fc(f)},jx{fc(f)});
        for k=1:length(kx)
            Theta(t,k)=Theta(t,k)+sum(sum(tmp(:,:,k).*hFacC{f}(:,:,k)))/vol(k);
        end
    end
end

% save time series
tme=datenum(ts2dte(ts,45,2011,1,17));
eval(['cd ~dmenemen/llc_' int2str(nx) '/regions/' region_name])
save Theta tme thk dpt vol Theta
