% This file is an example how to read TS time series

% Load grid information
cd ~dmenemen/llc_4320/regions/SPURS1/grid
nx=144;
ny=115;
nz=89;
suf=[int2str(nx) 'x' int2str(ny)];
Depth=readbin(['Depth_' suf],[nx ny]);
XC   =readbin([   'XC_' suf],[nx ny]);
YC   =readbin([   'YC_' suf],[nx ny]);
load ~dmenemen/llc_4320/grid/thk90

% Integration period
Time=datenum('13-Sep-2011'):(1/24):datenum('15-Nov-2012');
nt=length(Time);

% Create list of file names and indices
FileNames={};
Indices=XC*0;
n=1;
for x=1:24:nx
    mx=min(nx,x+23);
    for y=1:23:ny
        my=min(ny,y+22);
        FileNames{n}=['.x' num2str(x) '-' num2str(mx) ...
                    '.y' num2str(y) '-' num2str(my) ...
                    '.nz' num2str(nz) '.nt' num2str(nt)];
        Indices(x:mx,y:my)=n;
        n=n+1;
    end
end

% Extract time series at location (I,J)
cd ~dmenemen/llc_4320/regions/SPURS1/tiles
I=70;
J=60;
fnm=FileNames{Indices(I,J)};
c=textscan(fnm,'%*2c %d %*c %d %*c %d %*c %d %*2c %d %*2c %d');
ix=c{1}:c{2};
jx=c{3}:c{4};
lon=XC(ix,jx);
lat=YC(ix,jx);
salt =single(readbin([ 'Salt/Salt'  fnm],[length(ix) length(jx) nz nt]));
theta=single(readbin(['Theta/Theta' fnm],[length(ix) length(jx) nz nt]));
i=find(lon(:,1)==XC(I,J));
j=find(lat(1,:)==YC(I,J));

% Plot top 300 m
iz=find(dpt90<300);
clf reset
colormap jet
subplot(211)
mypcolor(Time,-dpt90(iz),squeeze(theta(i,j,iz,:)))
datetick('x',4,'keeplimits')
title(['Theta at ' int2str(lon(i,j)) 'E, ' int2str(lat(i,j)) 'N'])
ylabel('Depth (m)')
colorbar;
subplot(212)
mypcolor(Time,-dpt90(iz),squeeze(salt(i,j,iz,:)))
datetick('x',4,'keeplimits')
title(['Salt at ' int2str(lon(i,j)) 'E, ' int2str(lat(i,j)) 'N'])
xlabel('Month of 2011/2012')
ylabel('Depth (m)')
colorbar;