% {{{ Initialize a few things
nx=8640; ny=4320;                          % model grid dimensions
rSphere=6370e3;                            % radius of Earth (m)
m2deg=360/2/pi/rSphere;                    % meters to degrees at Equator
dx=1/24;                                   % grid spacing at Equator (deg)
dt=1/24;                                   % model output time step (days)
DT=dt*24*60*60;                            % model output time step (s)
dte=datenum(2011,3,6):dt:datenum(2013,4,22,6,0,0); % model output times
pin1='~dmenemen/llc_2160/regions/latlon/'; % llc2160 input file location
pin2='~dmenemen/llc_2160/regions/latlon/drift/transcoast/nguu/'; % output location
suf=['_' int2str(nx) 'x' int2str(ny)];
for fld={'XC','YC','XG','YG','hFacC'}     % load model grid indices and mask
  fnm=[pin1 'grid/' fld{1} suf];
  eval([fld{1} '=single(readbin(fnm,[nx ny]));'])
end
xc=XC(:,500);                              % longitude of grid centers
xg=XG(:,500);                              % longitude of grid corners
xc(find(xc<xc(1)))=xc(find(xc<xc(1)))+360;
xg(find(xg<xg(1)))=xg(find(xg<xg(1)))+360;
yc=YC(6000,:);                             % latitude of grid centers
yg=YG(6000,:);                             % latitude of grid corners
% }}}

% Specify drifters integration length and release times and locations.
% {{{ Release locations transcoast Josphat Nguu
nme='transcoast_nguu';
DriftLength=90;                           % drifters integration length (days)
fin=[pin2 'transcoastnguu_input.txt'];
eval(['tmp=load(''' fin ''');'])          % load release locations
Release.Latitude=tmp(:,1);
Release.Longitude=tmp(:,2);
Release.Time=[(datenum(2012,3,1):dt:datenum(2012,6,30))'];
% }}}

% {{{ Find closest wet point to release locations
LandMask=hFacC;
for i=1:length(Release.Longitude)
    tmp = (YC-Release.Latitude(i)).^2  + ...
          ((XC-Release.Longitude(i))*cos(pi*Release.Latitude(i)/180)).^2;
    tmp(find(LandMask==0))=1e6;
    [I J]=min(tmp(:));
    Release.WetLat(i)=YC(J);
    Release.WetLon(i)=XC(J);
    LandMask(J)=0;
end
disp([Release.Latitude Release.Longitude Release.WetLat' Release.WetLon'])
ix = find( xc>(min(Release.Longitude)-2*dx) & ...
           xc<(max(Release.Longitude)+2*dx));
iy = find( yc>(min(Release.Latitude)-2*dx) & ...
           yc<(max(Release.Latitude)+2*dx));
clf
mypcolor(xc(ix),yc(iy),hFacC(ix,iy)');
hold on
plot(Release.Longitude,Release.Latitude,'k*')
plot(Release.WetLon,Release.WetLat,'ko')
close all
% }}}

% Drifters release times
tmp=ones(length(Release.Latitude),1)*Release.Time';
DriftInitDte=tmp(:);

% Drifter release longitudes
DriftLon=zeros(length(DriftInitDte),DriftLength*24+1,'single');
tmp=Release.WetLon'*ones(1,length(Release.Time));
DriftLon(:,1)=tmp(:);

% Drifter release latitudes
DriftLat=zeros(length(DriftInitDte),DriftLength*24+1,'single');
tmp=Release.WetLat'*ones(1,length(Release.Time));
DriftLat(:,1)=tmp(:);

% Then run following code:
% {{{ Compute and save hourly drifter locations
DriftIdx=ones(length(DriftInitDte),1);
suf=['_' int2str(nx) 'x' int2str(ny) 'x1.'];
NN=0;
for t=min(DriftInitDte):dt:(max(DriftInitDte)+DriftLength-dt)
    disp(datestr(t))
    it=find(t>=DriftInitDte & t<(DriftInitDte+DriftLength));
    if length(it)>0
        NN=NN+1;
        
        % find index of current drifter positions
        Idx=sub2ind(size(DriftLon),it,DriftIdx(it));
    
        % limit latitudes
        tlat=DriftLat(Idx);
        tlat(find(tlat<yc(1)))=yc(1);
        tlat(find(tlat>yg(end)))=yg(end);
        DriftLat(Idx)=tlat;
        
        % unwrap longitudes
        tlon=DriftLon(Idx);
        in=find(tlon<xc(1));
        if length(in)>0
            tlon(in)=tlon(in)+360;
        end
        in=find(tlon>(xc(1)+360));
        if length(in)>0
            tlon(in)=tlon(in)-360;
        end
        DriftLon(Idx)=tlon;
        
        % if possible, try to reduce region needed for interpolation
        ix=find((xc+dx)>=min(DriftLon(Idx)) & (xg-dx)<=max(DriftLon(Idx)));
        iy=find((yc+dx)>=min(DriftLat(Idx)) & (yg-dx)<=max(DriftLat(Idx)));

        % wrap-around indices, if needed
        if any(DriftLon(Idx)>xg(end)), ix=[ix; 1]; end
        
        % read in U and V fields
        fnu=[pin1 'U/U' suf datestr(t,30)];
        tmp=readbin(fnu,[nx ny]);
        U=tmp(ix,iy);                          % zonal velocity (m/s)
        fnv=[pin1 'V/V' suf datestr(t,30)];
        tmp=readbin(fnv,[nx ny]);
        V=tmp(ix,iy);                          % meridional velocity (m/s)
        
        % bilinear interpolation of velocity to tracer locations
        lonc=xc(ix);
        long=xg(ix);
        if any(DriftLon(Idx)>xg(end))
            lonc(end)=lonc(end)+360;
            long(end)=long(end)+360;
        end        
        u=interp2(yc(iy),long,U,DriftLat(Idx),DriftLon(Idx));
        v=interp2(yg(iy),lonc,V,DriftLat(Idx),DriftLon(Idx));
        
        % find index of next drifter positions
        DriftIdx(it)=DriftIdx(it)+1;
        Idx2=sub2ind(size(DriftLon),it,DriftIdx(it));
        
        % integrate drifters
        DriftLon(Idx2)=DriftLon(Idx)+u*DT*m2deg./cos(pi*DriftLat(Idx)/180);
        DriftLat(Idx2)=DriftLat(Idx)+v*DT*m2deg;
    
    end
    if NN==100
        eval(['save ' pin2 nme ' DriftInitDte DriftLat DriftLon'])
        NN=0;
    end        
end
% }}}
eval(['save ' pin2 nme ' DriftInitDte DriftLat DriftLon'])