% This script reads in the ctd cnl files
clear
%dir    = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/ctd_raw/converted/';
dir    = '/Volumes/tgt-data/ctd/converted_ts2/';
outdir = '~/research/cruise/2019_Thompson_TN368_0703-0719/data/ctd_processed/';
fnpre  = 'tn368';


nfile  = [1:47 49:99]; 
% Converted with version 9 xmlcon; 
append = '';

ndata  = 19; % number of columns in data

%Nd = 200; % use 200 points for the vertical resolution at each depth
maxdpth = 650; % max depth for stations in a sigle transect has to be the same
idpth   = 0:-1:-maxdpth;

% lat and lon for Stations, Now they are dierctly read from files
% lon = [-70.823 -70.823];
% lat = [41.391  41.391];
lon = zeros(size(nfile));
lat = lon;

for ifile = 1:length(nfile)
    % simpler way!
     filename = [dir fnpre '-' num2str(nfile(ifile),'%03d') '_' append '1mBins.cnv'];
  try  
    fid = fopen(filename);
  catch  % in case the .cnv file doesn't exist
    continue  % skip and move to the next file
  end
  % skip the first 10 lines, figure out from the data file
    for skipi = 1:9
         tline = fgets(fid);
    end
   % then read in the latitude and longitude
    [lat1 count] = fscanf(fid,'* NMEA Latitude = %d %f N\n',[1,2]);
    [lon1 count] = fscanf(fid,'* NMEA Longitude = %d %f W\n',[1,2]);
   % now combine two
   lat(ifile) = lat1(1) + lat1(2)/60;
   lon(ifile) = -(lon1(1) + lon1(2)/60);
   % then read in the time
    [timestr count] = fscanf(fid,'* NMEA UTC (Time) = %c%c%c %d %d  %d:%d:%d\n',[1,8]);
    switch char(timestr(1:3))
        case 'Jan'
            month = 1;
        case 'Feb' 
            month = 2;
        case 'Mar'
            month = 3;
        case 'Apr'
            month = 4;
        case 'May'
            month = 5;
        case 'Jun'
            month = 6;
        case 'Jul'
            month = 7;
        case 'Aug'
            month = 8;
        case 'Sep'
            month = 9;
        case 'Oct'
            month = 10;
        case 'Nov'
            month = 11;
        case 'Dec'
            month = 12;
        otherwise
            disp(['Wrong month info in the cnv file: ' filename]);
            error('Please check');
    end
    time(ifile) = datenum([timestr(5) month timestr(4) timestr(6) timestr(7) timestr(8)]);
    % skip and read how many data are written
    for skipi = 1:13
         tline = fgets(fid);
    end
   [StaA{ifile},count] = fscanf(fid,'** A%s\n',1); % station ID
   fclose(fid);
end

dpthfile = '~/research/roms/bathymetry/ngdc2.mat';
bdata = load(dpthfile);
wdpth = griddata(bdata.lon,bdata.lat,bdata.h,lon,lat);

for ifile = 1:length(nfile)
    % simpler way!
     filename = [dir fnpre '-' num2str(nfile(ifile),'%03d') '_' append '1mBins.cnv'];
  try  
    fid = fopen(filename);
  catch  % in case the .cnv file doesn't exist
    continue  % skip and move to the next file
  end
   
  % this line automatically determine how many lines to ignore
  idx = 0;
  while ~feof(fid) %Check every line in the file
    templine = fgetl(fid); %Read one line...
    if (strcmp(templine(1),'*') || strcmp(templine(1),'#') )
        idx = idx+1;
    else
        break; % no need to go on if useful data is found
    end
  end
  fclose(fid);
  % print out to check if the skip lines and ndata are OK, checked!
 
   fid = fopen(filename);
   % read all data into 1D array, because we know the format 
    yy = textscan(fid, '%f','Delimiter','() ',...
      'MultipleDelimsAsOne',true, 'HeaderLines', idx);
    fclose(fid);
    data = yy{:};
    data(data == -9.99e-29) = nan; % replace the filled value with NaNs
    %data = reshape(yy,[],ndata);
    pres = data(2:ndata:end);
    temp = data(3:ndata:end); %(data(2:ndata:end)+data(3:ndata:end))./2;
    salt = data(5:ndata:end); %(data(12:ndata:end)+data(13:ndata:end))./2;
    cond = data(4:ndata:end); %(data(4:ndata:end)+data(5:ndata:end))./2;    % unit: S/m
    flor = data(8:ndata:end);
    par  = data(11:ndata:end);
    turb = data(12:ndata:end);
    oxyg = data(9:ndata:end);
    ntrt = data(13:ndata:end);
    cdom = data(7:ndata:end);
    
    % remove the surface 3 m
    nnind = find(pres < 3);
    if ~isempty(nnind)
        pres(nnind) = [];
        temp(nnind) = [];
        salt(nnind) = [];
        cond(nnind) = [];
        flor(nnind) = [];
        turb(nnind) = [];
        par(nnind)  = [];
        oxyg(nnind) = [];  
        ntrt(nnind) = [];  
        cdom(nnind) = [];  
    end
    
    % remove the reverse motion
    dpres = pres(2:end)-pres(1:end-1);
    nnind = find(dpres <= 0);
    while ~isempty(nnind)
        pres(nnind+1) = [];
        temp(nnind+1) = [];
        salt(nnind+1) = [];
        cond(nnind+1) = [];
        flor(nnind+1) = [];
        turb(nnind+1) = [];
        par(nnind+1)  = [];
        oxyg(nnind+1) = [];
        ntrt(nnind+1) = [];
        cdom(nnind+1) = [];
        
        dpres = pres(2:end)-pres(1:end-1);
        nnind = find(dpres <= 0);
    end
    
    % remove the NaNs before the interpolation
    tpres = pres; 
    nnind = find(isnan(temp));
    temp(nnind)  = [];
    tpres(nnind) = [];
    itemp = interp1(tpres,temp,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(salt));
    salt(nnind)  = [];
    tpres(nnind) = [];
    isalt = interp1(tpres,salt,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(cond));
    cond(nnind)  = [];
    tpres(nnind) = [];
    icond = interp1(tpres,cond,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(flor));
    flor(nnind)  = [];
    tpres(nnind) = [];
    iflor = interp1(tpres,flor,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(turb));
    turb(nnind)  = [];
    tpres(nnind) = [];
    iturb = interp1(tpres,turb,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(par));
    par(nnind)  = [];
    tpres(nnind) = [];
    ipar  = interp1(tpres,par,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(oxyg));
    oxyg(nnind)  = [];
    tpres(nnind) = [];
    ioxyg = interp1(tpres,oxyg,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(ntrt));
    ntrt(nnind)  = [];
    tpres(nnind) = [];
    intrt = interp1(tpres,ntrt,-idpth,'linear');

    tpres = pres; 
    nnind = find(isnan(cdom));
    cdom(nnind)  = [];
    tpres(nnind) = [];
    icdom = interp1(tpres,cdom,-idpth,'linear');

    ipden = sw_pden(isalt,itemp,abs(idpth),1);

    ilon   = lon(ifile);
    ilat   = lat(ifile);
    iwdpth = wdpth(ifile);
    istaA  = StaA{ifile};

    omatfn = [outdir fnpre num2str(nfile(ifile),'%03d') '.mat'];
    save(omatfn,'idpth','itemp','isalt','icond','ipden','iflor','iturb','ipar', ...
         'ioxyg','intrt','icdom','filename','ilon','ilat','iwdpth','istaA');
    disp(['Save file: ' omatfn]);
    disp('');
end