cd ~dmenemen/llc_2160/regions/latlon/drift
load global_v3/Release.mat
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={'XG','YG'}                        % load model grid info
    fnm=[pin 'grid/' fld{1} suf];
    eval([fld{1} '=readbin(fnm,[nx ny]);'])
end
xg=XG(:,500);                              % longitude of grid corners
yg=YG(6000,:);                             % latitude of grid corners
xg(find(xg<xg(1)))=xg(find(xg<xg(1)))+360;
xg=[xg; 360];
clear *G c*
r=1
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
[n idx]=histc(Lon(:),xg);
[n idy]=histc(Lat(:),yg);
idx=reshape(idx,size(Lon));            % longitude index (relative to xg)
idy=reshape(idy,size(Lat));            % latitude index  (relative to yg)
IT=sub2ind([nx ny],idx,idy);
Release.SpeciesIndices=Release.Idx;
Release.Idx=IT;
save global/Release Release DriftLength