%% Read Argo files from the website http://www.coriolis.eu.org/cdc/argo.htm
% These data are in netCDF format, from which this script
% extracts the position, depth, time and cycle number 
% for each useable Argo probe contained in the file. 
% The data is filtered for quality and some initial treatment performed 
% before the final data is then saved in the file probeData.mat

% Jon Yearsley & Julia Sigwart (Aug 2010)
% jon.yearsley@ucd.ie
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all

% The main data file containing the argo probe data in NETCDF format
fileName1 = 'CO_PF_20090618_120705.nc'; % Replace this with your file name

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% No more parameters to change below this line

% Read in the data for the probes from the NETCDF file
ncid = netcdf.open(fileName1,'nc_nowrite'); % Open NETDCF file

varid = netcdf.inqVarID(ncid,'PLATFORM_NUMBER');  % Probe ID
tmp = netcdf.getVar(ncid,varid);
% Convert tmp to a numerical probe ID
[n,m]=size(tmp);
floatID = zeros(m,1);
for i=1:m,
  floatID(i) = str2double(tmp(:,i));
end

varid = netcdf.inqVarID(ncid,'CYCLE_NUMBER'); % Cycle numer for a probe
cycle = netcdf.getVar(ncid,varid);

varid = netcdf.inqVarID(ncid,'LATITUDE'); % Latitude
latitude = netcdf.getVar(ncid,varid);

varid = netcdf.inqVarID(ncid,'LONGITUDE'); % Longitude
longitude = netcdf.getVar(ncid,varid);

varid = netcdf.inqVarID(ncid,'POSITION_QC');  % Position quality control
posquality = netcdf.getVar(ncid,varid);

varid = netcdf.inqVarID(ncid,'PRES');  % Pressure
pressure = netcdf.getVar(ncid,varid);
pressure(pressure==99999) = nan;
maxpressure = max(pressure,[],1)';

% Time at which lat and long recorded... may be later than JULD
% which is time of ascent
varid = netcdf.inqVarID(ncid,'JULD_LOCATION');
jday = netcdf.getVar(ncid,varid);

% Calculate date in years and then months (discarding the year info)
year = jday / 365.24 + 1950;
month = (year - floor(year))*12;

varid = netcdf.inqVarID(ncid,'JULD_QC'); % JULD quality control
dayqual = netcdf.getVar(ncid,varid);

netcdf.close(ncid)

floatList = unique(floatID); % A list of the different float ID's
nFloat = length(floatList);

%% Data subset

% Function to calculate geodetic distance between points
R_earth = 6357; % Radius of Earth in km
distFun = inline('R_av*acos(cos(lat(:,1)).*cos(lat(:,2)).*cos(lon(:,1)-lon(:,2)) + sin(lat(:,1)).*sin(lat(:,2)))','lat','lon','R_av');

float = zeros(length(latitude),1);
pres =  zeros(length(latitude),1);
latStart = zeros(length(latitude),1);
longStart = zeros(length(latitude),1);
dateStart = zeros(length(latitude),1);
latEnd = zeros(length(latitude),1);
longEnd = zeros(length(latitude),1);
dateEnd = zeros(length(latitude),1);
quality = zeros(length(latitude),1);
tdif = zeros(length(latitude),1);
dist = zeros(length(latitude),1);
useData = zeros(length(latitude),1);

% Remove obvious errors and unwanted points in the data
% Look through the data and pull out cycles that are less than 18 days
% long, less than 3000dbar max pressure, and less than 500km drift
% distances... also exclude points with a velocity greater than 30km/day
count = 1;
for f=1:nFloat
   indFloat = find(floatID==floatList(f));
   for t=1:(length(indFloat)-1),
       difft = jday(indFloat(t+1)) - jday(indFloat(t));
       if difft>4 && posquality(indFloat(t+1))=='1' && posquality(indFloat(t))=='1'
           latStart(count) = latitude(indFloat(t)); 
           longStart(count) = longitude(indFloat(t)); 
           dateStart(count) = jday(indFloat(t)); 
           latEnd(count) = latitude(indFloat(t+1)); 
           longEnd(count) = longitude(indFloat(t+1)); 
           dateEnd(count) = jday(indFloat(t+1)); 

           tdif(count) = difft;
           quality(count) = 1;
           float(count) = f;
           pres(count) = maxpressure(indFloat(t));
           
           dist(count) = distFun([latStart(count), latEnd(count)]*pi/180,[longStart(count), longEnd(count)]*pi/180,R_earth);
           if difft<18 && dist(count)<500 && pres(count)<3000 && dist(count)/difft < 30; 
               % Only keep data which is one cycle and for which the distance travelled is less than 500km
               useData(count) = 1;
           end
           count = count + 1;
       end
       
   end
end

% Starting lat, long and date for each probe cycle
latStart = latStart(1:(count-1));
longStart = longStart(1:(count-1));
dateStart = dateStart(1:(count-1));

% Finishing lat, long and date for each probe cycle
latEnd = latEnd(1:(count-1));
longEnd = longEnd(1:(count-1));
dateEnd = dateEnd(1:(count-1));

float = float(1:(count-1));
quality = quality(1:(count-1));
tdif = tdif(1:(count-1));
dist = dist(1:(count-1));



%% Calculate some final quantities and save the probe data 
% in the file 'probeData.mat'

pos = zeros(nnz(useData),2);
fl = zeros(nnz(useData),1);
depth = zeros(nnz(useData),1);
velocity = zeros(nnz(useData),1);
year = zeros(nnz(useData),1);
month = zeros(nnz(useData),1);
timeStep = zeros(nnz(useData),1);
distStep = zeros(nnz(useData),1);
longlatStart = zeros(nnz(useData),2);
longlatEnd = zeros(nnz(useData),2);

ind = find(useData);
for i=1:length(ind)
    fl(i) = float(ind(i));
    depth(i) = pres(ind(i));
    
    % Calculate true data from Julian data
    tmp = dateStart(ind(i))/365.24 + 1950;
    year(i) = floor(tmp); % Year of probe data
    month(i) = ceil((tmp - year(i))*12); % Month of probe data
    timeStep(i) = tdif(ind(i)); % Number of days between probe diving and probe ascending
    distStep(i) = dist(ind(i)); % Great circle distance travelled by each probe during 1 cycle
    longlatStart(i,:) = [longStart(ind(i)) latStart(ind(i))]; % Long and lat of probe dive
    longlatEnd(i,:) = [longEnd(ind(i)) latEnd(ind(i))]; % Long and lat of probe resurface
        
    % The longitude and latitude for the midpoint of each probe cycle
    pos(i,1) = (longEnd(ind(i)) + longStart(ind(i)))/2; 
    pos(i,2) = (latEnd(ind(i)) + latStart(ind(i)))/2;
    
    % An estimate of the average velocity of the probe (km/day) per cycle
    velocity(i) = distStep(i) / timeStep(i);
end

% Save the data
save 'probeData.mat' longlatStart longlatEnd velocity pos year month timeStep distStep depth



