cd ~dmenemen/llc_4320/regions/BottomPressure

% load bottom pressure data
load Stations
nx=4320;
prec='real*4';
gdir='~dmenemen/llc_4320/grid/';
fnam=[gdir 'YC.data'];
pin='~dmenemen/llc_4320/MITgcm/run/';
ts=10368:144:279360;
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),25,2011,9,10));
end

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

    % extract indices in the model
    fac=cos(min(56,abs(lat(s)))*pi/180)/48;
    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 StationsObsPlusModel 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