cd ~/llc_4320/regions/BottomPressure

% load bottom pressure data
load Stations
clear dte ppr
nx=2160;
prec='real*4';
gdir='~dmenemen/llc_2160/grid/';
fnam=[gdir 'YC.data'];
pin='~dmenemen/llc_2160/MITgcm/run';
ts=92160:80:1586400;
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)
        if (ts(t)<1198080)
            fin=[pin '_day49_624/' fnm '.' myint2str(ts(t),10) '.data'];
        else 
            fin=[pin '/' fnm '.' myint2str(ts(t),10) '.data'];
        end        
        mppr{s}(t)=read_llc_fkij(fin,nx,fc,1,ix,jx);
    end
    in=find(mppr{s}==0);
    mppr{s}(in)=nan;
end

% remove first NaN
mdte=mdte(2:end);
for i=1:length(stn)
    mppr{i}=mppr{i}(2:end);
end

% check that all time series have same NaNs
in=find(~isnan(mppr{1}));
for i=2:length(stn)
    in2=find(~isnan(mppr{i}));
    if length(in)~=length(in2) | any(in-in2)
        error('in is not equal to in2')
    end
end

% fill in NaNs using interp1
for i=1:length(stn)
    mppr{i}=interp1(mdte(in),mppr{i}(in),mdte);
end

save StationsModel2160long lat lon mdte mppr stn