root/gliderproc/trunk/gliderOptode_Generate_L1_Data.m

%
%  gliderOptode_Generate_L1_Data.m
%
%  Purpose: Generate Level 1 data mat files that include optode DO data.
%
%  NOTE: At top of file, set glider index and deployment number
%
%  Requires:
%  MATLAB folder - contains util, matutil, seawater, plots, strfun, opnml
%
%
%  Author:  William Stark
%           Marine Sciences Department
%           UNC-Chapel Hill
%
%  Created: 21 May 2012
%
%//////////////////////////////////////////////////////////////////////////

clear all;


%!!!!!!  SET PRIOR TO RUNNING CODE  !!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%
% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ...
gliderIndex = 1;

% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ...
deploymentNumber = 1;

%
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



% add paths for required files...
addpath('MATLAB/util/');
addpath('MATLAB/matutil/');
addpath('MATLAB/seawater/');
addpath('MATLAB/plots/');
addpath('MATLAB/strfun/');
addpath('MATLAB/opnml/');
addpath('MATLAB/opnml/FEM/');



% glider name string...
if (gliderIndex==1)
    strGliderName = 'Pelagia';
else
    strGliderName = 'Ramses';
end

% populate arrays for the deployment start and end dates...
% ex. strStart(2, 3) is start date for Ramses, Deployment 3
strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'};
strEnd   = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'};

% deployment number string...
strDeploymentNumber = num2str(deploymentNumber);

% deployment start date string...
strStartDate = strStart(gliderIndex, deploymentNumber);

% deployment end date string...
strEndDate = strEnd(gliderIndex, deploymentNumber);

% define the path to the glider ascii files...
%datadir = strcat('/Users/haloboy/Documents/MASC/MATLAB/CTD_data_correction/GLIDER_CTD_DATA_LEVEL0/',...
datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, '_Deployment', strDeploymentNumber, '/');



% define default bounds for use in plots...
% switch gliderIndex
%     case 1  % Pelagia
%         switch deploymentNumber
%             case 1  % Deployment 1
%                 tempBounds =  [17.0 24.0];
%                 salinBounds = [36.0 36.4];
%                 densBounds =  [1025.0 1026.6];
%                 chlorBounds = [0.0 4.0];
%                 
%             case 2  % Deployment 2
%                 tempBounds =  [17.0 24.0];
%                 salinBounds = [36.0 36.5];
%                 densBounds =  [1024.5 1026.8];
%                 chlorBounds = [0.0 4.0];
%                 
%             case 3  % Deployment 3
%                 tempBounds =  [17.0 24.0];
%                 salinBounds = [35.9 36.7];
%                 densBounds =  [1024.4 1026.4];
%                 chlorBounds = [0.0 4.0];
%         end
%     case 2  % Ramses
%         switch deploymentNumber
%             case 1  % Deployment 1
%                 tempBounds =  [8.0 23.0];
%                 salinBounds = [35.0 36.4];
%                 densBounds =  [1024.5 1027.5];
%                 chlorBounds = [0.0 4.0];
%             
%             case 2  % Deployment 2
%                 tempBounds =  [9.0 25.0];
%                 salinBounds = [35.2 36.6];
%                 densBounds =  [1024.0 1027.5];
%                 chlorBounds = [0.0 4.0];
%             
%             case 3  % Deployment 3
%                 tempBounds =  [10.0 24.5];
%                 salinBounds = [35.3 36.7];
%                 densBounds =  [1024.4 1027.4];
%                 chlorBounds = [0.0 4.0];
%         end
% end


%##########################################################################################





%*** READ IN EBD DATA *****************************************************
% declare variables for storing data...
oxy_dphase=[];
ptime_oxy=[];

% try to load all *.ebdasc files at once...
[files, Dstruct] = wilddir(datadir, '.ebdasc');
nfile = size(files, 1);

for i=1:nfile-1
    % protect against empty ebd file
    if(Dstruct(i).bytes>0)
        data = read_gliderasc2([datadir, files(i,:)]);
        %data = read_gliderasc3([datadir, files(i,:)]);

        % if the number of values (in data.data) is less than the number of 
        % vars (in data.vars), this means that the data were not completely read
        % in.  To correct this, pad data.data with NaNs until its length
        % equals that of data.vars...
        if (length(data.data) < length(data.vars))
            data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post');
        end

        % populate variables with data...
        if(~isempty(data.data))
            oxy_dphase = [oxy_dphase;data.data(:,strmatch('sci_oxy3835_wphase_dphase', data.vars, 'exact'))];   % oxygen (dphase)
            ptime_oxy = [ptime_oxy;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))];       % present time
            %should read in other oxy values for comparison with
            %re-calculated values
        end

        data = [];
    end  
end
%**************************************************************************


% load glider data from existing mat file...
strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, '_Deployment', strDeploymentNumber, '_DataL1.mat');
load(strMatFilePath);


% first, apply the sort() function to ptime_oxy to make sure it increases monotonically...
[Y,I] = sort(ptime_oxy);
ptime_oxy = Y;
oxy_dphase = oxy_dphase(I);


% remove nans from EBD data...
i = find(~isnan(oxy_dphase));
oxy_dphase = oxy_dphase(i); ptime_oxy = ptime_oxy(i);

% interpolate CTD to times of oxy_dphase 
tempi = interp1(ptime,temp,ptime_oxy);
salini = interp1(ptime,salinCorrected,ptime_oxy);
depthi = interp1(ptime,depth,ptime_oxy);

% first, implement temperature-dependent correction to DO concentration
% that utilizes dphase as the input (manual page 30), coefficients from cal
% sheets (a 5x4 matrix of values) and will be glider (and foil) dependent.

% for ramses, 2 June 2010 calibration sheet

C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2];
C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3];
C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5];
C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7];
C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10];

C0 = C0coef(1)+C0coef(2)*tempi+C0coef(3)*tempi.^2+C0coef(4)*tempi.^3;
C1 = C1coef(1)+C1coef(2)*tempi+C1coef(3)*tempi.^2+C1coef(4)*tempi.^3;
C2 = C2coef(1)+C2coef(2)*tempi+C2coef(3)*tempi.^2+C2coef(4)*tempi.^3;
C3 = C3coef(1)+C3coef(2)*tempi+C3coef(3)*tempi.^2+C3coef(4)*tempi.^3;
C4 = C4coef(1)+C4coef(2)*tempi+C4coef(3)*tempi.^2+C4coef(4)*tempi.^3;

o2 = C0+C1.*oxy_dphase+C2.*oxy_dphase.^2+C3.*oxy_dphase.^3+C4.*oxy_dphase.^4;

% second, implement the salinity correction to DO concentration (page 31)
% ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER
% SCALING IS REQUIRED

B0=-6.24097e-3;
B1=-6.93498e-3;
B2=-6.90358e-3;
B3=-4.29155e-3;
C0=-3.11680e-7;

Ts = log((298.15 - tempi) ./ (273.15 + tempi));

o2_salcorr = o2*exp(salini*(B0+B1*Ts+B2*Ts.*Ts+B3*Ts.*Ts.*Ts)+C0*salini.^2);

% third, implement the pressure correction to DO concentration (page 32)

o2_salprescorr = os_salcorr*(1 + 0.04*depthi/1000);

% use polynomial to calculate DO satuations using the measured temp and
% sal (manual page 30)

A0 = 2.00856;
A1 = 3.22400;
A2 = 3.99063;
A3 = 4.80299;
A4 = 0.0978188;
A5 = 1.71069;
B0 = -00624097;
B1 = -00693498;
B2 = -00690358;
B3 = -00429155;
C0 = -000000311680;

% need interpolated temp, salinity, pressure at times of DO obs

rslt = A0 + A1*Ts + A2*Ts.^2 + A3*Ts.^3 + A4*Ts.^4 + A5*Ts.^5 +...
    salini.*(B0 + B1*Ts + B2*Ts.^2 + B3*Ts.^3) + C0*salini.^2;

O2sol = exp(rslt);

o2_sat = (o2_salprescorr * 2.2414) ./ O2sol;


% convert ptime into datenum style...
ptime_oxy_datenum = ptime_oxy/3600/24+datenum(1970, 1, 1, 0, 0, 0);
% ptimehrly = fix(ptime_datenum*24)/24;
% ptimedaily = fix(ptime_datenum);
% ptimedaily2 = unique(ptimedaily);
% ptimedaily2 = ptimedaily2(1:2:end);




lb = datenum('10-Feb-2012 12:22:00');
ub = datenum('10-Feb-2012 12:23:00');
%lb = datenum('22-Mar-2012 12:01:00');
%ub = datenum('23-Mar-2012 12:02:00');

% figure('Position', [500,500,1000,1000]);
% subplot(311);
% ind = find(ptime_oxy_datenum >= lb & ptime_oxy_datenum <= ub);
% plot(ptime_oxy_datenum(ind), oxy_dphase(ind), 'b.');
% datetick;
% xlabel('time', 'fontsize', 12);
% ylabel('DO  (dphase)', 'fontsize', 12);
% subplot(312);
% ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub);
% plot(ptime_datenum(ind2), temp(ind2), 'r.');
% datetick;
% xlabel('time', 'fontsize', 12);
% ylabel('temp', 'fontsize', 12);
% subplot(313);
% ind = find(ptime_oxy_datenum >= lb & ptime_oxy_datenum <= ub);
% ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub);
% plot(ptime_oxy_datenum(ind), ones(length(ind), 1), 'b.', ptime_datenum(ind2), ones(length(ind2), 1), 'r.');
% datetick;


figure('Position', [500,500,1000,1000]);
subplot(211);
ind = find(ptime_datenum >= lb & ptime_datenum <= ub);
plot(ptime_datenum(ind), oxy_dphase(ind), 'b.', ptime_datenum(ind), oxy_dphase_adjusted(ind), 'r.');
datetick;
xlabel('time', 'fontsize', 12);
ylabel('DO  (dphase)', 'fontsize', 12);
subplot(212);
ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub);
plot(ptime_datenum(ind2), temp(ind2), 'r.');
datetick;
xlabel('time', 'fontsize', 12);
ylabel('temp', 'fontsize', 12);








% interpolate DBD data to align with EBD data...
% horizontalVelocity = interp1(ptime1_dbd, horizontalVelocity, ptime);
% depth = interp1(ptime1_dbd, depth, ptime);
% pitch = interp1(ptime1_dbd, pitch, ptime);
% avgDepthRate = interp1(ptime1_dbd, avgDepthRate, ptime);
% glideAngle = interp1(ptime1_dbd, glideAngle, ptime);
% lat = interp1(ptime1_dbd, lat, ptime);
% lon = interp1(ptime1_dbd, lon, ptime);












% create configuration struct...
% units = struct('chlor', 'micrograms liter-1',...
%                'dens', 'kg m-3',...
%                'densCorrected', 'kg m-3',...
%                'depth', 'm',...
%                'gpsLat', 'decimal degrees',...
%                'gpsLon', 'decimal degrees',...
%                'pres', 'decibars',...
%                'ptime', 'seconds since 0000-01-01T00:00',...
%                'ptime_datenum', 'seconds since 0000-01-01T00:00',...
%                'salin', 'psu',...
%                'salinCorrected', 'psu',...
%                'temp', 'deg C');

% variable_description = struct('chlor', 'chlorophyll measured by glider',...
%                               'chlorBounds', 'default chlorophyll bounds for plots',...
%                               'dens', 'density measured by glider',...
%                               'densBounds', 'default density bounds for plots',...
%                               'densCorrected', 'density corrected for thermal lag',...
%                               'depth', 'depth calculated as function of pressure and position latitude',...
%                               'gpsLat', 'position latitude measured by glider GPS',...
%                               'gpsLon', 'position longitude measured by glider GPS',...
%                               'pres', 'pressure measured by glider',...
%                               'ptime', 'time vector reported by glider',...
%                               'ptime_datenum', 'Serial Date Number string',...
%                               'salin', 'salinity measured by glider',...
%                               'salinBounds', 'default salinity bounds for plots',...
%                               'salinCorrected', 'salinity corrected for thermal lag',...
%                               'temp', 'temperature measured by glider',...
%                               'tempBounds', 'default temperature bounds for plots');

% correction_parameters = struct('alpha_offset', correctionParams(1),...
%                                'alpha_slope', correctionParams(2),...
%                                'tau_offset', correctionParams(3),...
%                                'tau_slope', correctionParams(4));
           
% config = struct('glider_name', strGliderName,...
%                 'deployment_number', strDeploymentNumber,...
%                 'start_date', strStartDate,...
%                 'end_date', strEndDate,...
%                 'thermal_lag_correction_parameters', correction_parameters,...
%                 'var_descriptions', variable_description,...
%                 'var_units', units);




% set Level 1 data mat file name...
%strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_WithDO_LEVEL1.mat');



% save glider/deployment data to mat file...
% save(strMatFileName,...
%      'config',...
%      'chlor',...
%      'chlorBounds',...
%      'dens',...
%      'densBounds',...
%      'densCorrected',...
%      'depth',...
%      'gpsLat',...
%      'gpsLon',...
%      'pres',...
%      'ptime',...
%      'ptime_datenum',...
%      'salin',...
%      'salinBounds',...
%      'salinCorrected',...
%      'temp',...
%      'tempBounds');