cd /nobackupp8/cgshuman/BottomPressure

% load bottom pressure data
load Stations
nx=2160;
prec='real*4';
gdir='~dmenemen/llc_2160/grid/';
fnam=[gdir 'YC.data'];
pin='~dmenemen/llc_2160/MITgcm/run_day49_624/';
ts=dte2ts(ts2dte(10368,25,2011,9,10),45,2011,1,17):80:dte2ts(ts2dte(279360,25,2011,9,10),45,2011,1,17);
fnm='PhiBot';
mppr=[]; % model bottom pressure

% convert model time step to a matlab date
mdte=0*ts;
for t=1:length(ts)
	mdte(t)=datenum(ts2dte(ts(t),45,2011,1,17));
end

for s=1:length(stn), disp(stn(s))

    % extract indices in the model
    fac=cos(min(56,abs(lat(s)))*pi/180)/24;
    minlat=lat(s)-fac;
    maxlat=lat(s)+fac;
    minlon=lon(s)-fac;
    maxlon=lon(s)+fac;
    [YC fc ix jx] = ...
        quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);
    if length(ix) == 0
        error('increase fac')
    elseif length(ix)>1 | length(jx)>1
        XC=read_llc_fkij([gdir 'XC.data'],nx,fc,1,ix,jx);
        [JX IX]=meshgrid(jx,ix);
        [Y I]=min(sqrt((XC(:)-lon(s)).^2+(YC(:)-lat(s)).^2));
        ix=IX(I);
        jx=JX(I);
    end

    % extract model-equivalent time series
    for t=1:length(ts), mydisp(t)
        fin=[pin fnm '.' myint2str(ts(t),10) '.data'];
        mppr{s}(t)=read_llc_fkij(fin,nx,fc,1,ix,jx);
    end
    in=find(mppr{s}==0);
    mppr{s}(in)=nan;
end

save StationsObsPlusModel2160 lat lon dte mdte ppr mppr stn

% plot results
%for s=1:length(stn)
%    clf
%    plot(dte{s},ppr{s}-mmean(ppr{s}),mdte,(mppr{s}-mmean(mppr{s}))/9.81)
%    legend('data','model')
%    ndaytick
%    title(['Station ' int2str(stn(s))])
%    pause
%end