cd ~dmenemen/llc_2160/regions/latlon/drift
load global_v3/Release.mat
dt=1/24;                                   % model output time step (days)
eps=dt/100;                                % epsilon used for roundoff errors
DT=dt*24*60*60;                            % model output time step (s)

% extract and save trajectories
for t=Release.Tme(1):dt:(Release.Tme(end)+DriftLength)
    disp(datestr(t))

    % determine active drifters at this instant in time
    iDrift=find(Release.Tme<=t+eps & Release.Tme>=(t-DriftLength-eps));

    if length(iDrift)>0
        % read location of drifters at this instant in time
        suf2=['_' int2str(length(Release.Lon)) 'x' int2str(length(iDrift)) '_' datestr(t,30)];
        flon=['global_v3/DriftLon' suf2];
        flat=['global_v3/DriftLat' suf2];
        DriftLon=readbin(flon,[length(Release.Lon),length(iDrift)]);
        DriftLat=readbin(flat,[length(Release.Lon),length(iDrift)]);
        for i=1:length(iDrift)
            suf2=['_n' int2str(length(Release.Lon)) '_dt' int2str(DT) '_' datestr(Release.Tme(iDrift(i)),30)];
            flon=['global/DriftLon' suf2];
            flat=['global/DriftLat' suf2];
            writebin(flon,DriftLon(:,i),1,'real*4',round((t-Release.Tme(iDrift(i)))/dt));
            pause(.002)
            writebin(flat,DriftLat(:,i),1,'real*4',round((t-Release.Tme(iDrift(i)))/dt));
            pause(.002)
        end
    end
end


%% fill-in incomplete trajectories
% 
%tmp={'20-Jun-2011 18:00:00',
%     '20-Mar-2012 18:00:00'
%    }
%for i=1:length(tmp)
%    startDte(i)=datenum(tmp{i});
%end
%for s=1:length(startDte)
%    dte=startDte(s):dt:(startDte(s)+DriftLength);
%    Lat=zeros(length(Release.WetLon),length(dte));
%    Lon=zeros(length(Release.WetLon),length(dte));
%    for t=1:length(dte), mydisp(t)
%        iDrift=find(Release.Tme<=dte(t)+eps & Release.Tme>=(dte(t)-DriftLength-eps));
%        suf=['_' int2str(length(Release.WetLon)) 'x' int2str(length(iDrift)) '_' datestr(dte(t),30)];
%        if exist(['global_v3/DriftLon' suf])
%            it=find(Release.Tme(iDrift)==startDte(s));
%            Lon(:,t)=readbin(['global_v3/DriftLon' suf],length(Release.WetLon),1,'real*4',it-1);
%            Lat(:,t)=readbin(['global_v3/DriftLat' suf],length(Release.WetLon),1,'real*4',it-1);
%        else
%            error('choubichou')
%        end
%    end
%    writebin(['DriftLon_n4331_dt3600_' datestr(startDte(s),30)],Lon);
%    writebin(['DriftLat_n4331_dt3600_' datestr(startDte(s),30)],Lat);
%end