% {{{ 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));
Release.WetLon=XC(ix);                     % drifters release longitudes
Release.WetLat=YC(ix);                     % drifters release latitudes
Release.Tme=datenum(2012,1,1);             % drifters release times
% }}}
% {{{ 2. Eleven release locations for Magdalene
nme='magdalene';
DriftLength=90;                            % drifters integration length (days)
% {{{ Desired release locations and times
Release.Lat = [4.771911
               4.492892
               4.380511
               3.964747
               4.058339
               4.089933
               4.114333
               4.077972
               3.755764
               3.082686
               2.335969];
Release.Lon = [8.798011
               8.930833
               8.896236
               9.266128
               9.382239
               9.722011
               9.749464
               9.704531
               9.700072
               9.974861
               9.855525];
Release.Tme = ...
[ (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.Lon)
    tmp = (YC-Release.Lat(i)).^2  + ...
          ((XC-Release.Lon(i))*cos(pi*Release.Lat(i)/180)).^2;
    tmp(find(LandMask==0))=1e6;
    [I J]=min(tmp(:));
    Release.WetLat(i,1)=YC(J);
    Release.WetLon(i,1)=XC(J);
    LandMask(J)=0;
end
disp([Release.Lat Release.Lon Release.WetLat Release.WetLon])
ix = find( xc>(min(Release.Lon)-2*dx) & ...
           xc<(max(Release.Lon)+2*dx));
iy = find( yc>(min(Release.Lat)-2*dx) & ...
           yc<(max(Release.Lat)+2*dx));
clf
mypcolor(xc(ix),yc(iy),hFacC(ix,iy)');
hold on
plot(Release.Lon,Release.Lat,'k*')
plot(Release.WetLon,Release.WetLat,'ko')
close all
% }}}
% }}}
% {{{ 3. Global biogeographical mangrove range
nme='global';
DriftLength=365;                           % drifters integration length (days)
fin=[pin 'drift/JPL_NASA_MRDH_VLIZ_for_modelling_speciesnumericID.txt'];
eval(['tmp=load(''' fin ''');'])           % load release locations
Release.Lat=tmp(:,2);
Release.Lon=tmp(:,3);
Release.Tme=datenum(2011,4,1):dt:datenum(2012,4,1);

% remove duplicate locations
for i=1:length(Release.Lon)
    if i>=length(Release.Lon), break, end
    ix=find(Release.Lon==Release.Lon(i)&Release.Lat==Release.Lat(i));
    if length(ix)>1
        Release.Lon(ix(2:end))=[];
        Release.Lat(ix(2:end))=[];
    end
end

% find closest wet point to release location
eval(['load ' pin 'grid/coastline'])       % load model coastline indices
ix=find(coastline);
for i=1:length(Release.Lon)
    tmp = (YC(ix)-Release.Lat(i)).^2  + ...
          ((XC(ix)-Release.Lon(i))*cos(pi*Release.Lat(i)/180)).^2;
    [I J]=min(tmp(:));
    Release.WetLat(i,1)=YC(ix(J));
    Release.WetLon(i,1)=XC(ix(J));
    ix(J)=1;
end

clf reset
orient landscape
wysiwyg
plot(XC(ix),YC(ix),'r.','markersize',24)
hold on
plot(Release.WetLon,Release.WetLat,'k.','markersize',24)
plot(Release.Lon,Release.Lat,'.','markersize',24)
axis([-180 180 -40 35])
legend('Coastline Points','Release Locations','Mangrove Atlas',0)
eval(['print -dpsc ' pin 'figs/ReleaseLocations'])
tmp=double([Release.Lon Release.Lat Release.WetLon Release.WetLat]);
eval(['save -ascii ' pin 'figs/ReleaseLocations.txt tmp'])
close all
% }}}

% Then run following code:
% {{{ Compute and save hourly drifter locations
eval(['mkdir ' pin 'drift/' nme])
suf=['_' int2str(nx) 'x' int2str(ny) 'x1.'];
DriftLon=Release.WetLon;
DriftLat=Release.WetLat;
minDriftLon=min(DriftLon);
minDriftLat=min(DriftLat);
maxDriftLon=max(DriftLon);
maxDriftLat=max(DriftLat);

for t=Release.Tme(1):dt:(Release.Tme(end)+DriftLength)
    disp(datestr(t))
    
    % if possible, try to reduce region needed for interpolation
    ix=find((xc+dx)>=minDriftLon & (xg-dx)<=maxDriftLon);
    iy=find((yc+dx)>=minDriftLat & (yg-dx)<=maxDriftLat);

    % wrap-around indices, if needed
    if any(DriftLon>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)
    
    % reset minimum and maximum drifter longitudes and latitudes
    minDriftLon=xc(1)+360;
    minDriftLat=yg(end);
    maxDriftLon=xc(1);
    maxDriftLat=yc(1);    
    
    for DriftInitDte=max(t-DriftLength,Release.Tme(1)):dt:t

        % load DriftLon and DriftLat
        suf2=['_n' int2str(length(DriftLon)) '_dt' int2str(DT) '_' datestr(DriftInitDte,30)];
        flon=[pin 'drift/' nme '/DriftLon' suf2];
        flat=[pin 'drift/' nme '/DriftLat' suf2];
        if DriftInitDte==t
            DriftLon=Release.WetLon;
            DriftLat=Release.WetLat;
        else
            DriftLon=readbin(flon,length(DriftLon),1,'real*4',round((t-DriftInitDte)/dt));
            DriftLat=readbin(flat,length(DriftLat),1,'real*4',round((t-DriftInitDte)/dt));
        end
    
        % limit latitudes
        DriftLat(find(DriftLat<yc(1)))=yc(1);
        DriftLat(find(DriftLat>yg(end)))=yg(end);
        
        % unwrap longitudes
        in=find(DriftLon<xc(1));
        if length(in)>0
            DriftLon(in)=DriftLon(in)+360;
        end
        in=find(DriftLon>(xc(1)+360));
        if length(in)>0
            DriftLon(in)=DriftLon(in)-360;
        end

        % save initial latitudes and longitudes to output files
        if DriftInitDte==t
            writebin(flon,DriftLon);
            writebin(flat,DriftLat);
        end

        % bilinear interpolation of velocity to tracer locations
        lonc=xc(ix);
        long=xg(ix);
        if any(DriftLon>xg(end))
            lonc(end)=lonc(end)+360;
            long(end)=long(end)+360;
        end        
        u=interp2(yc(iy),long,U,DriftLat,DriftLon);
        v=interp2(yg(iy),lonc,V,DriftLat,DriftLon);

        % integrate drifters
        DriftLon=DriftLon+u*DT*m2deg./cos(pi*DriftLat/180);
        DriftLat=DriftLat+v*DT*m2deg;

        % write new drifter locations to file
        writebin(flon,DriftLon,1,'real*4',round((t-DriftInitDte)/dt)+1);
        writebin(flat,DriftLat,1,'real*4',round((t-DriftInitDte)/dt)+1);

        % update minimum and maximum drifter longitudes and latitudes
        minDriftLon=min([DriftLon(:); minDriftLon]);
        minDriftLat=min([DriftLat(:); minDriftLat]);
        maxDriftLon=max([DriftLon(:); maxDriftLon]);
        maxDriftLat=max([DriftLat(:); maxDriftLat]);

    end
end
% }}}