%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
clc

thisDir = fileparts( mfilename('fullpath') );

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basic info on the Total files.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
conf = MABSB_conf;

% time period
startdate=datenum(2013,04,27,18,0,0);
finaldate = datenum(2013,04,27,23,0,0);
dtime = startdate:1/24:finaldate;

% Generate all filenames at once using repmat appropriately
fnames=datenum_to_directory_filename(conf.WRFOUT.BaseDir,dtime,...
                                     conf.WRFOUT.FilePrefix,conf.WRFOUT.FileSuffix,conf.WRFOUT.MonthFlag);

%
% load RUWRF Files


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do some particle tracking on the grid.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp( 'Doing some particle tracking on a rectangular Lon,Lat grid.' );

% Domain: MAB: d01: SEA_BREEZE (MABSB)
% generate the drifter release points
lonlims = [-79.52 -68.48]; % (degree)
latlims = [35.37 43.52]; % (degree)
heightlims=[11 3000]; % (meter)

% % Domain: MAB: d02: SEA_BREEZE (MABSB)
% % generate the drifter release points
% lonlims = [-75.98 -72.95]; % (degree)
% latlims = [38.53 40.97]; % (degree)
% heightlims=[11 3000]; % (meter)

dlonlat=[10.0,10.0];
% dlonlat=[1.0, 1.0];
% dlonlat=[0.5,0.5];
% dlonlat=[0.10, 0.10];
du='km';
%
[lon_MABSB_tracer,lat_MABSB_tracer] = LonLat_grid([min(lonlims),min(latlims)],[max(lonlims),max(latlims)],dlonlat,du );
lon_MABSB_tracer=lon_MABSB_tracer(:);
lat_MABSB_tracer=lat_MABSB_tracer(:);
%
ln = lon_MABSB_tracer;
lt = lat_MABSB_tracer;
LL = [ ln(:), lt(:) ];


% Get some good option values for tracking - otherwise output is junk.
abs_tol = 1.0e-3; % Not sure about this
rel_tol = 1.0e-3; % Not sure about this
maxstep = 1/24; % 1 hour
options = odeset('RelTol',rel_tol,'AbsTol',abs_tol,'MaxStep',maxstep);

tspan = startdate:1/24/6:finaldate; % every hour

% Get some good option values for tracking - otherwise output is junk.
abs_tol = 1.0e-3; % Not sure about this
rel_tol = 1.0e-3; % Not sure about this
maxstep = 10*60*60; %0.2 hour in seconds
options = odeset('RelTol',rel_tol,'AbsTol',abs_tol,'MaxStep',maxstep);
odesolver='ode45'; % can change solver if want
tic
[lons,lats,depths,ts]=particle_track_ode_grid3d_LonLat(Lon,Lat, Depth, ...
                                                       U, V, W, TT, tspan, cc, options, odesolver );
toc

% update cc configuration
cc=[lons(end,:)',lats(end,:)',depths(end,:)'];
    
    % Awww
    T3 = TRAJ3struct( size(lats) );
    
    T3.Type = 'TRAJ3';
    T3.TrajectoryDomain='MEDITERRANEAN';
    T3.CreationInfo='mitGCM cooling';
    
    T3.Lon = lons';
    T3.Lat = lats';
    T3.Depth = depths';
    T3.TimeStamp = ts(:)';
    
    % initial positions
    T3.InitialLonLatDepth = [T3.Lon(:,1), T3.Lat(:,1), T3.Depth(:,1)];
    
    % final positions
    T3.FinalLonLatDepth = [T3.Lon(:,end), T3.Lat(:,end), T3.Depth(:,end)];
    
    % Duration of tracks
    T3.TrajectoryDuration = T3.TimeStamp(end)' - T3.TimeStamp(1);
    
    % Some metadata
    

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for iter=1:nrecords-1
    disp(sprintf('iteration: %05d',iter))
    % Read MITgcm U,V,W  meta/data files
    U=rdmds([GCM_UVWDir,'U/','U'],records(n1(iter):n2(iter))); U=permute(U,[2,1,3,4]);
    V=rdmds([GCM_UVWDir,'V/','V'],records(n1(iter):n2(iter))); V=permute(V,[2,1,3,4]);
    W=rdmds([GCM_UVWDir,'W/','W'],records(n1(iter):n2(iter))); W=permute(W,[2,1,3,4]);
    
    % vertical top boundary condition
    Wbc = W(:,:,1,:);
    Wbc(Wbc<0)=0;
    W(:,:,1,:)=Wbc;
    TT= tt(n1(iter):n2(iter));
    t0=tt(n1(iter));
    tf=tt(n2(iter));
    tdump=24*3; tdump=tdump*3600;
    tspan=t0:tdump:tf;
    % Get some good option values for tracking - otherwise output is junk.
    abs_tol = 1.0e-3; % Not sure about this
    rel_tol = 1.0e-3; % Not sure about this
    maxstep = 10*60*60; %0.2 hour in seconds
    options = odeset('RelTol',rel_tol,'AbsTol',abs_tol,'MaxStep',maxstep);
    odesolver='ode45'; % can change solver if want
    tic
    [lons,lats,depths,ts]=particle_track_ode_grid3d_LonLat(Lon,Lat, Depth, ...
        U, V, W, TT, tspan, cc, options, odesolver );
    toc
    
    % update cc configuration
    cc=[lons(end,:)',lats(end,:)',depths(end,:)'];
    
     % Awww
    T3 = TRAJ3struct( size(lats) );
    
    T3.Type = 'TRAJ3';
    T3.TrajectoryDomain='MEDITERRANEAN';
    T3.CreationInfo='mitGCM cooling';
    
    T3.Lon = lons';
    T3.Lat = lats';
    T3.Depth = depths';
    T3.TimeStamp = ts(:)';
    
    % initial positions
    T3.InitialLonLatDepth = [T3.Lon(:,1), T3.Lat(:,1), T3.Depth(:,1)];
    
    % final positions
    T3.FinalLonLatDepth = [T3.Lon(:,end), T3.Lat(:,end), T3.Depth(:,end)];
    
    % Duration of tracks
    T3.TrajectoryDuration = T3.TimeStamp(end)' - T3.TimeStamp(1);
    
    % Some metadata 
    
    
    % save the particle lon/lat/depth
    
    outputDir=thisDir;
    savefile=[outputDir,'xyz.',sprintf('%05d',iter),'.mat'];
    varList = 'lons lats depths';
    varList = strcat(varList,' ts');
    eval(['save ',savefile,' ',varList])
    
%     for i=1:36;
%         h=scatter(lons(:,i),lats(:,i),50,depths(:,i));colorbar;
%         title(sprintf('iteration: %d05',iter))
%         hold on;
%     end
%     hold off;
    
end

