cd ~dmenemen/llc_2160/regions/latlon/drift
load global_v3/Release.mat
N=length(Release.Lon);
pin='~dmenemen/llc_2160/regions/latlon/'; % model output files location
nx=8640; ny=4320;                         % model grid dimensions
suf=['_' int2str(nx) 'x' int2str(ny)];
for fld={'XC','YC','hFacC'}               % load model grid info
  fnm=[pin 'grid/' fld{1} suf];
  eval([fld{1} '=readbin(fnm,[nx ny]);'])
end
XC(find(XC<XC(1)))=XC(find(XC<XC(1)))+360;
xc=XC(:,500);                             % longitude of grid centers
yc=YC(6000,:);                            % latitude of grid corners
land=abs(hFacC-1);                        % array of ones for land and zeros for water

t=1;                                      % trajectory index
r=1;                                      % release time index
eval(['load CAP/Capture' myint2str(r,4)])
flat=['global/DriftLat_n4331_dt3600_' datestr(Release.Tme(r),30)];
flon=['global/DriftLon_n4331_dt3600_' datestr(Release.Tme(r),30)];
D=dir(flon);
Lon=readbin(flon,[N D.bytes/4/N]);        % propagule trajectories longitude
Lat=readbin(flat,[N D.bytes/4/N]);        % propagule trajectories latitude
ix=find(xc>=min(Lon(t,:)-.2)&xc<=max(Lon(t,:)+.2));
iy=find(yc>=min(Lat(t,:)-.2)&yc<=max(Lat(t,:)+.2));
clf
shadeland(xc(ix),yc(iy),land(ix,iy)',.8*[1 1 1]);
hold on
plot(Lon(t,:),Lat(t,:))
plot(Lon(t,1),Lat(t,1),'g*','linewidth',3,'markersize',16)
plot(XC(CaptureIdx{t}),YC(CaptureIdx{t}),'ro','linewidth',2)