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

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

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Basic MITGCM information (GCM)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
baseDir = fullfile( thisDir, 'mitGCMFiles/mediterranean_sea/');
GCM_datafile=fullfile(baseDir,'gridFiles/data');
GCM_UVWDir = fullfile( thisDir, 'mitGCMFiles/mediterranean_sea/cooling/');

grid=mitgcm_loadgrid(GCM_datafile);

Lon = grid.X;
Lat = grid.Y;
Depth = grid.Z;

% Set the initial tracer position
TLon_orig= 18;
TLat_orig = 36;
TDepth_orig =0;

% longitude in (degree)
delTLon   = ones(1,2)*0.01;
% latitude in (degree)
delTLat   = ones(1,2)*0.005;
% depth in (m)
delTDepth   = [10, 20, 20];

TLon=cumsum([TLon_orig delTLon]);
TLat=cumsum([TLat_orig delTLat]);
TDepth=cumsum([TDepth_orig delTDepth]);
[TLon,TLat,TDepth]=meshgrid(TLon,TLat,TDepth);
cc = [  TLon(:), TLat(:), TDepth(:) ];
nT=size(cc,1);

% 4- set the iteration number in (second)
yyyy=0; 
mm=0;
dd=30;
HH=0;
MM=0;
SS=0;

ts = yyyy*12*30*24*60*60+mm*30*24*60*60+dd*24*60*60+HH*60*60+MM*60+SS;

records=scanformdsfiles([GCM_UVWDir,'U/','U']);
tt=cumsum(ones(1,length(records)))*ts;

% due to memory limitations, we can hold a certain number of files in the
% memory. idx1 and idx2 are the start/end-points files
nrecords=length(records);
n1=1:nrecords-1;
n2=n1+1;

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

