% {{{ 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

% get species indices for each location
for i=1:length(Release.Lon)
    ix=find(tmp(:,2)==Release.Lat(i)&tmp(:,3)==Release.Lon(i));
    idx=sort(tmp(ix,1));
    idx(find(diff(idx)==0))=[];
    Release.Idx{i}=idx;
end

% for each species, get index of locations
for i=1:max(tmp(:,1))
    Release.SpeciesLocation{i}=[];
    for j=1:length(Release.Lon)
        if any(Release.Idx{j}==i)
            Release.SpeciesLocation{i}=[Release.SpeciesLocation{i} j];
        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
clear I* J L* X* Y* a* dte f* h* i r*
eval(['mkdir ' pin 'drift/' nme '_v3'])
suf=['_' int2str(nx) 'x' int2str(ny) 'x1.'];
DriftLon=Release.WetLon*ones(1,length(Release.Tme));
DriftLat=Release.WetLat*ones(1,length(Release.Tme));
eval(['save ' pin 'drift/' nme '_v3/Release Release DriftLength'])

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 & Release.Tme>=(t-DriftLength));

    if length(iDrift)>0

        % 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

        % write location of drifters at this instant in time
        suf2=['_' int2str(size(DriftLon,1)) 'x' int2str(length(iDrift)) '_' datestr(t,30)];
        flon=[pin 'drift/' nme '_v3/DriftLon' suf2];
        flat=[pin 'drift/' nme '_v3/DriftLat' suf2];
        writebin(flon,DriftLon(:,iDrift));
        writebin(flat,DriftLat(:,iDrift));
        
        % if possible, try to reduce region needed for interpolation
        ix=find((xc+dx)>=min(min(DriftLon(:,iDrift))) & (xg-dx)<=max(max(DriftLon(:,iDrift))));
        iy=find((yc+dx)>=min(min(DriftLat(:,iDrift))) & (yg-dx)<=max(max(DriftLat(:,iDrift))));
        
        % wrap-around indices, if needed
        WrapAroundFlag=0;
        if any(any(DriftLon(:,iDrift)>xg(end)))
            ix=[ix; 1; 2];
            WrapAroundFlag=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 WrapAroundFlag
            lonc(end-1)=lonc(end-1)+360;
            long(end-1)=long(end-1)+360;
            lonc(end  )=lonc(end  )+360;
            long(end  )=long(end  )+360;
        end
        u=interp2(yc(iy),long,U,DriftLat(:,iDrift),DriftLon(:,iDrift));
        v=interp2(yg(iy),lonc,V,DriftLat(:,iDrift),DriftLon(:,iDrift));

        if any(isnan(u)), error('danger-danger'), end
        if any(isnan(v)), error('danger-danger'), end
        
        % integrate drifters
        DriftLon(:,iDrift)=DriftLon(:,iDrift)+u*DT*m2deg./cos(pi*DriftLat(:,iDrift)/180);
        DriftLat(:,iDrift)=DriftLat(:,iDrift)+v*DT*m2deg;
        
    end
end
% }}}