cd ~dmenemen/llc_4320/regions/BottomPressure

% load bottom pressure data
load Station_52403.txt
dte=datenum(Station_52403(:,1),Station_52403(:,2),Station_52403(:,3), ...
            Station_52403(:,4),Station_52403(:,5),Station_52403(:,6));
ppr=Station_52403(:,8);
lat=4.502;
lon=145.592;
clear S*

% extract indices in the model
minlat=lat-1/96*cos(lat*pi/180);
maxlat=lat+1/96*cos(lat*pi/180);
minlon=lon-1/96*cos(lat*pi/180);
maxlon=lon+1/96*cos(lat*pi/180);
nx=4320;
prec='real*4';
gdir='/nobackupp8/dmenemen/llc/llc_4320/grid/';
fnam=[gdir 'YC.data'];
[fld fc ix jx] = ...
    quikread_llc(fnam,nx,1,prec,gdir,minlat,maxlat,minlon,maxlon);

% extract model-equivalent time series
pin='/nobackupp8/dmenemen/llc/llc_4320/MITgcm/run/';
ts=10368:144:279360;
mppr=zeros(length(ts),1);
fnm='PhiBot';
for t=1:length(ts), mydisp(t)
    fin=[pin fnm '.' myint2str(ts(t),10) '.data'];
    mppr(t)=read_llc_fkij(fin,nx,fc,1,ix,jx);
end
iz=find(mppr==0);
mppr(iz)=[];
ts(iz)=[];

% 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

% compare time series
plot(dte,ppr-mean(ppr),mdte,(mppr-mean(mppr))/9.81)
legend('data','model')

% save station data
clear a* f* g* jx i* ma* mi* nx pin pr* t*
save Station_52403

%%%%%%%%%%%%%%%%%%
cd ~dmenemen/llc_4320/regions/BottomPressure
load Station_52403

% compare time series
plot(dte,ppr-mean(ppr),mdte,(mppr-mean(mppr))/9.81)
legend('data','model')