% extract time series at surface and 168 m depth
% at 38.95 N, 185.08 E

% vertical levels (.5 and 171.8 m depth)
kx=[1 27];

% horizontal indices (-174.9271E and 38.9444N)
minlat=38.94;
maxlat=38.96;
minlon=185.06;
maxlon=185.10;
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);

% get and save regional fields
pin='~dmenemen/llc_2160/MITgcm/run';
ts=92160:80:1586400;
U=zeros(length(ts),2);
V=zeros(length(ts),2);
for fnm={'U','V'}
    for t=1:length(ts), mydisp(t)
        if (ts(t)<1198080)
            fin=[pin '_day49_624/' fnm{1} '.' myint2str(ts(t),10) '.data'];
        else 
            fin=[pin '/' fnm{1} '.' myint2str(ts(t),10) '.data'];
        end        
        eval([fnm{1} '(t,:)=read_llc_fkij(fin,nx,fc,kx,ix,jx);'])
    end
end

% get coordinates
dates = datenum(ts2dte(ts(1),45,2011,1,17)):(1/24): ...
      datenum(ts2dte(ts(end),45,2011,1,17));
Ulat = read_llc_fkij([gdir 'YC.data'],nx,fc,1,ix,jx);
Ulon = read_llc_fkij([gdir 'XG.data'],nx,fc,1,ix,jx);
Vlat = read_llc_fkij([gdir 'YG.data'],nx,fc,1,ix,jx);
Vlon = read_llc_fkij([gdir 'XC.data'],nx,fc,1,ix,jx);
for k=1:length(kx)
 depths(k) = readbin([gdir 'RC.data'],1,1,'real*4',kx(k)-1);
end

% plot and save velocities
cd ~dmenemen/llc_2160/regions/39N185E
clf
subplot(211)
plot(dates,U)
title(['U (m/s) at ' num2str(Ulon) 'E, ' num2str(Ulat) 'N'])
datetick('x',12)
legend([num2str(depths(1)) ' m'],[num2str(depths(2)) ' m'],'Location','Best')
subplot(212)
plot(dates,V)
title(['V (m/s) at ' num2str(Vlon) 'E, ' num2str(Vlat) 'N'])
datetick('x',12)
legend([num2str(depths(1)) ' m'],[num2str(depths(2)) ' m'],'Location','Best')
print -dpsc llc2160_39N185E
save llc2160_39N185E U* V* d*