% {{{ 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
pin='~dmenemen/llc_2160/regions/latlon/';  % model output files location
suf=['_' int2str(nx) 'x' int2str(ny)];
for fld={'XC','YC','XG','YG','hFacC'}     % load model grid indices and mask
  fnm=[pin '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.
% Choose only one of the options below:
% {{{ 1. Release at edges of domain to test wrap-around
nme='edges';
DriftLength=1;                             % drifters integration length (days)
ix=find(YC==mmin(YC) | YC==mmax(YC) | XC==mmin(xc) | (XC+360)==mmax(xc));
tmp=ones(length(ix),1)*datenum(2012,1,1);
DriftInitDte=tmp(:);                       % drifters release times
DriftLon=zeros(length(DriftInitDte),DriftLength*24+1,'single');
DriftLon(:,1)=XC(ix);                      % drifters release longitudes
DriftLat=zeros(length(DriftInitDte),DriftLength*24+1,'single');
DriftLat(:,1)=YC(ix);                      % drifters release latitudes
% }}}
% {{{ 2. Release at east edge of domain to test wrap-around
nme='east_edge';
DriftLength=1;                             % drifters integration length (days)
ix=find((XC+360)==mmax(xc));
tmp=ones(length(ix),1)*datenum(2012,1,1);
DriftInitDte=tmp(:);                       % drifters release times
DriftLon=zeros(length(DriftInitDte),DriftLength*24+1,'single');
DriftLon(:,1)=XC(ix);                      % drifters release longitudes
DriftLat=zeros(length(DriftInitDte),DriftLength*24+1,'single');
DriftLat(:,1)=YC(ix);                      % drifters release latitudes
% }}}
% {{{ 3. Release at west edge of domain to test wrap-around
nme='west_edge';
DriftLength=10;                            % drifters integration length (days)
ix=find(XC==mmin(xc));
tmp=ones(length(ix),1)*datenum(2012,1,1);
DriftInitDte=tmp(:);                       % drifters release times
DriftLon=zeros(length(DriftInitDte),DriftLength*24+1,'single');
DriftLon(:,1)=XC(ix);                      % drifters release longitudes
DriftLat=zeros(length(DriftInitDte),DriftLength*24+1,'single');
DriftLat(:,1)=YC(ix);                      % drifters release latitudes
% }}}
% {{{ 4. Tropical West Africa, first-of-month release

nme='west_africa';
DriftLength=90;                            % drifters integration length (days)
eval(['load ' pin 'grid/coastline'])       % load model coastline indices
ix=find(YC>=-30 & YC<30 & XC>=-18 & XC<18 & coastline);
tmp=ones(length(ix),1)*datenum(2012,1:12,1);
DriftInitDte=tmp(:);                       % drifters release times
DriftLon=zeros(length(DriftInitDte),DriftLength*24+1,'single');
tmp=XC(ix)*ones(1,12);
DriftLon(:,1)=tmp(:);                      % drifters release longitudes
DriftLat=zeros(length(DriftInitDte),DriftLength*24+1,'single');
tmp=YC(ix)*ones(1,12);
DriftLat(:,1)=tmp(:);                      % drifters release latitudes

% }}}
% {{{ 5. Eleven release locations for Magdalene
nme='magdalene';
DriftLength=90;                            % drifters integration length (days)
% {{{ Desired release locations and times
Release.Latitude = [4.771911
                    4.492892
                    4.380511
                    3.964747
                    4.058339
                    4.089933
                    4.114333
                    4.077972
                    3.755764
                    3.082686
                    2.335969];
Release.Longitude = [8.798011
                    8.930833
                    8.896236
                    9.266128
                    9.382239
                    9.722011
                    9.749464
                    9.704531
                    9.700072
                    9.974861
                    9.855525];
Release.Time = ...
[ (datenum(2012,2,1):(1/24):datenum(2012,5,1))'
  (datenum(2012,9,1):(1/24):datenum(2012,11,1))'];
% }}}
% {{{ 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.'];
for t=min(DriftInitDte):dt:(max(DriftInitDte)+DriftLength-dt)
    disp(datestr(t))
    it=find(t>=DriftInitDte & t<(DriftInitDte+DriftLength));
    if length(it)>0
        
        % 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=[pin 'U/U' suf datestr(t,30)];
        tmp=readbin(fnu,[nx ny]);
        U=tmp(ix,iy);                          % zonal velocity (m/s)
        fnv=[pin '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
end
% }}}
eval(['save ' nme ' DriftInitDte DriftLat DriftLon'])